dune-istl 3.0-git
owneroverlapcopy.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_ISTL_OWNEROVERLAPCOPY_HH
4#define DUNE_ISTL_OWNEROVERLAPCOPY_HH
5
6#include <new>
7#include <iostream>
8#include <vector>
9#include <list>
10#include <map>
11#include <set>
12
13#include "cmath"
14
15// MPI header
16#if HAVE_MPI
17#include <mpi.h>
18#endif
19
20
21#include <dune/common/tuples.hh>
22#include <dune/common/enumset.hh>
23
24#if HAVE_MPI
25#include <dune/common/parallel/indexset.hh>
26#include <dune/common/parallel/communicator.hh>
27#include <dune/common/parallel/remoteindices.hh>
28#include <dune/common/parallel/mpicollectivecommunication.hh>
29#endif
30
31#include "solvercategory.hh"
32#include "istlexception.hh"
33#include <dune/common/parallel/collectivecommunication.hh>
35
36template<int dim, template<class,class> class Comm>
38
39
40namespace Dune {
41
63
75 template <class G, class L>
77 {
78 public:
80 typedef G GlobalIdType;
81
83 typedef L LocalIdType;
84
99
106 {
110 DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set");
111 localindices.insert(x);
112 }
113
120 {
124 DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set");
125 remoteindices.insert(x);
126 }
127
132 const std::set<IndexTripel>& localIndices () const
133 {
134 return localindices;
135 }
136
141 const std::set<RemoteIndexTripel>& remoteIndices () const
142 {
143 return remoteindices;
144 }
145
149 void clear ()
150 {
151 localindices.clear();
152 remoteindices.clear();
153 }
154
155 private:
157 std::set<IndexTripel> localindices;
159 std::set<RemoteIndexTripel> remoteindices;
160 };
161
162
163#if HAVE_MPI
164
171 template <class GlobalIdType, class LocalIdType=int>
173 {
174 template<typename M, typename G, typename L>
175 friend void loadMatrixMarket(M&,
176 const std::string&,
178 bool);
179 // used types
182 typedef typename std::set<IndexTripel>::const_iterator localindex_iterator;
183 typedef typename std::set<RemoteIndexTripel>::const_iterator remoteindex_iterator;
186 public:
190 typedef typename RI::RemoteIndex RX;
192 typedef Dune::Interface IF;
197 protected:
198
199
201 template<typename T>
203 {
205
206 static V gather(const T& a, std::size_t i)
207 {
208 return a[i];
209 }
210
211 static void scatter(T& a, V v, std::size_t i)
212 {
213 a[i] = v;
214 }
215 };
216 template<typename T>
218 {
220
221 static V gather(const T& a, std::size_t i)
222 {
223 return a[i];
224 }
225
226 static void scatter(T& a, V v, std::size_t i)
227 {
228 a[i] += v;
229 }
230 };
231
233 {
234 if (OwnerOverlapToAllInterfaceBuilt)
235 OwnerOverlapToAllInterface.free();
240 OwnerOverlapToAllInterface.build(ri,sourceFlags,destFlags);
241 OwnerOverlapToAllInterfaceBuilt = true;
242 }
243
245 {
246 if (OwnerToAllInterfaceBuilt)
247 OwnerToAllInterface.free();
250 OwnerToAllInterface.build(ri,sourceFlags,destFlags);
251 OwnerToAllInterfaceBuilt = true;
252 }
253
255 {
256 if (OwnerCopyToAllInterfaceBuilt)
257 OwnerCopyToAllInterface.free();
262 OwnerCopyToAllInterface.build(ri,sourceFlags,destFlags);
263 OwnerCopyToAllInterfaceBuilt = true;
264 }
265
267 {
268 if (OwnerCopyToOwnerCopyInterfaceBuilt)
269 OwnerCopyToOwnerCopyInterface.free();
273 OwnerCopyToOwnerCopyInterface.build(ri,sourceFlags,destFlags);
274 OwnerCopyToOwnerCopyInterfaceBuilt = true;
275 }
276
278 {
279 if (CopyToAllInterfaceBuilt)
280 CopyToAllInterface.free();
283 CopyToAllInterface.build(ri,sourceFlags,destFlags);
284 CopyToAllInterfaceBuilt = true;
285 }
286
287 public:
288
293 category = set;
294 }
295
301 return category;
302 }
303
305 {
306 return cc;
307 }
308
315 template<class T>
316 void copyOwnerToAll (const T& source, T& dest) const
317 {
318 if (!OwnerToAllInterfaceBuilt)
321 communicator.template build<T>(OwnerToAllInterface);
323 communicator.free();
324 }
325
332 template<class T>
333 void copyCopyToAll (const T& source, T& dest) const
334 {
335 if (!CopyToAllInterfaceBuilt)
338 communicator.template build<T>(CopyToAllInterface);
340 communicator.free();
341 }
342
349 template<class T>
350 void addOwnerOverlapToAll (const T& source, T& dest) const
351 {
352 if (!OwnerOverlapToAllInterfaceBuilt)
355 communicator.template build<T>(OwnerOverlapToAllInterface);
357 communicator.free();
358 }
359
366 template<class T>
367 void addOwnerCopyToAll (const T& source, T& dest) const
368 {
369 if (!OwnerCopyToAllInterfaceBuilt)
372 communicator.template build<T>(OwnerCopyToAllInterface);
374 communicator.free();
375 }
376
383 template<class T>
384 void addOwnerCopyToOwnerCopy (const T& source, T& dest) const
385 {
386 if (!OwnerCopyToOwnerCopyInterfaceBuilt)
389 communicator.template build<T>(OwnerCopyToOwnerCopyInterface);
391 communicator.free();
392 }
393
394
402 template<class T1, class T2>
403 void dot (const T1& x, const T1& y, T2& result) const
404 {
405 // set up mask vector
406 if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size()))
407 {
408 mask.resize(x.size());
409 for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
410 mask[i] = 1;
411 for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
412 if (i->local().attribute()!=OwnerOverlapCopyAttributeSet::owner)
413 mask[i->local().local()] = 0;
414 }
415 result = T2(0.0);
416
417 for (typename T1::size_type i=0; i<x.size(); i++)
418 result += x[i]*(y[i])*mask[i];
419 result = cc.sum(result);
420 return;
421 }
422
429 template<class T1>
431 {
432 // set up mask vector
433 if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size()))
434 {
435 mask.resize(x.size());
436 for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
437 mask[i] = 1;
438 for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
439 if (i->local().attribute()!=OwnerOverlapCopyAttributeSet::owner)
440 mask[i->local().local()] = 0;
441 }
442 typename T1::field_type result = typename T1::field_type(0.0);
443 for (typename T1::size_type i=0; i<x.size(); i++)
444 result += x[i].two_norm2()*mask[i];
445 return static_cast<double>(sqrt(cc.sum(result)));
446 }
447
449
452
455
459
465 {
466 return pis;
467 }
468
474 {
475 return ri;
476 }
477
483 {
484 return pis;
485 }
486
487
493 {
494 return ri;
495 }
496
498 {
499 if(globalLookup_) {
500 if(pis.seqNo()==oldseqNo)
501 // Nothing changed!
502 return;
503 delete globalLookup_;
504 }
505
506 globalLookup_ = new GlobalLookupIndexSet(pis);
507 oldseqNo = pis.seqNo();
508 }
509
510 void buildGlobalLookup(std::size_t size)
511 {
512 if(globalLookup_) {
513 if(pis.seqNo()==oldseqNo)
514 // Nothing changed!
515 return;
516 delete globalLookup_;
517 }
518 globalLookup_ = new GlobalLookupIndexSet(pis, size);
519 oldseqNo = pis.seqNo();
520 }
521
523 {
524 delete globalLookup_;
525 globalLookup_=0;
526 }
527
529 {
530 assert(globalLookup_ != 0);
531 return *globalLookup_;
532 }
533
539 template<class T1>
540 void project (T1& x) const
541 {
542 for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
543 if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy)
544 x[i->local().local()] = 0;
545 }
546
558 bool freecomm_ = false)
559 : comm(comm_), cc(comm_), pis(), ri(pis,pis,comm_),
560 OwnerToAllInterfaceBuilt(false), OwnerOverlapToAllInterfaceBuilt(false),
561 OwnerCopyToAllInterfaceBuilt(false), OwnerCopyToOwnerCopyInterfaceBuilt(false),
562 CopyToAllInterfaceBuilt(false), globalLookup_(0), category(cat_),
563 freecomm(freecomm_)
564 {}
565
575 : comm(MPI_COMM_WORLD), cc(MPI_COMM_WORLD), pis(), ri(pis,pis,MPI_COMM_WORLD),
576 OwnerToAllInterfaceBuilt(false), OwnerOverlapToAllInterfaceBuilt(false),
577 OwnerCopyToAllInterfaceBuilt(false), OwnerCopyToOwnerCopyInterfaceBuilt(false),
578 CopyToAllInterfaceBuilt(false), globalLookup_(0), category(cat_), freecomm(false)
579 {}
580
589 MPI_Comm comm_,
591 bool freecomm_ = false)
592 : comm(comm_), cc(comm_), OwnerToAllInterfaceBuilt(false),
593 OwnerOverlapToAllInterfaceBuilt(false), OwnerCopyToAllInterfaceBuilt(false),
594 OwnerCopyToOwnerCopyInterfaceBuilt(false), CopyToAllInterfaceBuilt(false),
595 globalLookup_(0), category(cat_), freecomm(freecomm_)
596 {
597 // set up an ISTL index set
598 pis.beginResize();
599 for (localindex_iterator i=indexinfo.localIndices().begin(); i!=indexinfo.localIndices().end(); ++i)
600 {
606 pis.add(get<0>(*i),LI(get<1>(*i),OwnerOverlapCopyAttributeSet::copy,true));
607 // std::cout << cc.rank() << ": adding index " << get<0>(*i) << " " << get<1>(*i) << " " << get<2>(*i) << std::endl;
608 }
609 pis.endResize();
610
611 // build remote indices WITHOUT communication
612 // std::cout << cc.rank() << ": build remote indices" << std::endl;
613 ri.setIndexSets(pis,pis,cc);
614 if (indexinfo.remoteIndices().size()>0)
615 {
616 remoteindex_iterator i=indexinfo.remoteIndices().begin();
617 int p = get<0>(*i);
618 RILM modifier = ri.template getModifier<false,true>(p);
619 typename PIS::const_iterator pi=pis.begin();
620 for ( ; i!=indexinfo.remoteIndices().end(); ++i)
621 {
622 // handle processor change
623 if (p!=get<0>(*i))
624 {
625 p = get<0>(*i);
626 modifier = ri.template getModifier<false,true>(p);
627 pi=pis.begin();
628 }
629
630 // position to correct entry in parallel index set
631 while (pi->global()!=get<1>(*i) && pi!=pis.end())
632 ++pi;
633 if (pi==pis.end())
634 DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set");
635
636 // insert entry
637 // std::cout << cc.rank() << ": adding remote index " << get<0>(*i) << " " << get<1>(*i) << " " << get<2>(*i) << std::endl;
644 }
645 }else{
646 // Force remote indices to be synced!
647 ri.template getModifier<false,true>(0);
648 }
649 }
650
651 // destructor: free memory in some objects
653 {
654 ri.free();
655 if (OwnerToAllInterfaceBuilt) OwnerToAllInterface.free();
656 if (OwnerOverlapToAllInterfaceBuilt) OwnerOverlapToAllInterface.free();
657 if (OwnerCopyToAllInterfaceBuilt) OwnerCopyToAllInterface.free();
658 if (OwnerCopyToOwnerCopyInterfaceBuilt) OwnerCopyToOwnerCopyInterface.free();
659 if (CopyToAllInterfaceBuilt) CopyToAllInterface.free();
660 if (globalLookup_) delete globalLookup_;
661 if (freecomm==true)
662 if(comm!=MPI_COMM_NULL)
663 {
664#ifdef MPI_2
665 // If it is possible to query whether MPI_Finalize
666 // was called, only free the communicator before
667 // calling MPI_Finalize.
668 int wasFinalized = 0;
670 if(!wasFinalized)
671#endif
672 MPI_Comm_free(&comm);
673 }
674 }
675
676 private:
678 {}
679 MPI_Comm comm;
680 CollectiveCommunication<MPI_Comm> cc;
681 PIS pis;
682 RI ri;
683 mutable IF OwnerToAllInterface;
684 mutable bool OwnerToAllInterfaceBuilt;
685 mutable IF OwnerOverlapToAllInterface;
686 mutable bool OwnerOverlapToAllInterfaceBuilt;
687 mutable IF OwnerCopyToAllInterface;
688 mutable bool OwnerCopyToAllInterfaceBuilt;
689 mutable IF OwnerCopyToOwnerCopyInterface;
690 mutable bool OwnerCopyToOwnerCopyInterfaceBuilt;
691 mutable IF CopyToAllInterface;
692 mutable bool CopyToAllInterfaceBuilt;
693 mutable std::vector<double> mask;
694 int oldseqNo;
695 GlobalLookupIndexSet* globalLookup_;
697 bool freecomm;
698 };
699
700#endif
701
702
705} // end namespace
706
707#endif
Provides classes for reading and writing MatrixMarket Files with an extension for parallel matrices.
void testRedistributed(int s)
Definition basearray.hh:19
Statistics about compression achieved in implicit mode.
Definition bcrsmatrix.hh:81
derive error class from the base class in common
Definition istlexception.hh:16
Attribute set for overlapping schwarz.
Definition owneroverlapcopy.hh:58
AttributeSet
Definition owneroverlapcopy.hh:59
@ owner
Definition owneroverlapcopy.hh:60
@ copy
Definition owneroverlapcopy.hh:60
@ overlap
Definition owneroverlapcopy.hh:60
Information about the index distribution.
Definition owneroverlapcopy.hh:77
tuple< GlobalIdType, LocalIdType, int > IndexTripel
A triple describing a local index.
Definition owneroverlapcopy.hh:91
void addRemoteIndex(const RemoteIndexTripel &x)
Add a new remote index triple to the set of remote indices.
Definition owneroverlapcopy.hh:119
G GlobalIdType
The type of the global index.
Definition owneroverlapcopy.hh:80
const std::set< IndexTripel > & localIndices() const
Get the set of indices local to the process.
Definition owneroverlapcopy.hh:132
const std::set< RemoteIndexTripel > & remoteIndices() const
Get the set of remote indices.
Definition owneroverlapcopy.hh:141
L LocalIdType
The type of the local index.
Definition owneroverlapcopy.hh:83
void clear()
Remove all indices from the sets.
Definition owneroverlapcopy.hh:149
tuple< int, GlobalIdType, int > RemoteIndexTripel
A triple describing a remote index.
Definition owneroverlapcopy.hh:98
void addLocalIndex(const IndexTripel &x)
Add a new index triple to the set of local indices.
Definition owneroverlapcopy.hh:105
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition owneroverlapcopy.hh:173
EnumItem< AttributeSet, OwnerOverlapCopyAttributeSet::copy > CopySet
Definition owneroverlapcopy.hh:194
const GlobalLookupIndexSet & globalLookup() const
Definition owneroverlapcopy.hh:528
FieldTraits< typenameT1::field_type >::real_type norm(const T1 &x) const
Compute the global euclidian norm of a vector.
Definition owneroverlapcopy.hh:430
void buildOwnerOverlapToAllInterface() const
Definition owneroverlapcopy.hh:232
Dune::ParallelIndexSet< GlobalIdType, LI, 512 > PIS
Definition owneroverlapcopy.hh:187
void buildOwnerCopyToAllInterface() const
Definition owneroverlapcopy.hh:254
void buildOwnerCopyToOwnerCopyInterface() const
Definition owneroverlapcopy.hh:266
OwnerOverlapCopyCommunication(const IndexInfoFromGrid< GlobalIdType, LocalIdType > &indexinfo, MPI_Comm comm_, SolverCategory::Category cat_=SolverCategory::overlapping, bool freecomm_=false)
Constructor.
Definition owneroverlapcopy.hh:588
void addOwnerCopyToOwnerCopy(const T &source, T &dest) const
Communicate values from owner and copy data points to owner and copy data points and add them to thos...
Definition owneroverlapcopy.hh:384
void buildCopyToAllInterface() const
Definition owneroverlapcopy.hh:277
Dune::EnumItem< AttributeSet, OwnerOverlapCopyAttributeSet::copy > CopyFlags
Definition owneroverlapcopy.hh:448
RemoteIndices & remoteIndices()
Get the underlying remote indices.
Definition owneroverlapcopy.hh:492
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition owneroverlapcopy.hh:464
Dune::RemoteIndices< PIS > RI
Definition owneroverlapcopy.hh:188
void buildGlobalLookup(std::size_t size)
Definition owneroverlapcopy.hh:510
void addOwnerOverlapToAll(const T &source, T &dest) const
Communicate values from owner data points to all other data points and add them to those values.
Definition owneroverlapcopy.hh:350
Dune::RemoteIndices< PIS > RemoteIndices
The type of the remote indices.
Definition owneroverlapcopy.hh:454
void project(T1 &x) const
Set vector to zero at copy dofs.
Definition owneroverlapcopy.hh:540
Dune::AllSet< AttributeSet > AllSet
Definition owneroverlapcopy.hh:196
Combine< EnumItem< AttributeSet, OwnerOverlapCopyAttributeSet::owner >, EnumItem< AttributeSet, OwnerOverlapCopyAttributeSet::overlap >, AttributeSet > OwnerOverlapSet
Definition owneroverlapcopy.hh:195
void copyCopyToAll(const T &source, T &dest) const
Communicate values from copy data points to all other data points.
Definition owneroverlapcopy.hh:333
Dune::GlobalLookupIndexSet< ParallelIndexSet > GlobalLookupIndexSet
The type of the reverse lookup of indices.
Definition owneroverlapcopy.hh:458
Dune::Interface IF
Definition owneroverlapcopy.hh:192
~OwnerOverlapCopyCommunication()
Definition owneroverlapcopy.hh:652
void buildGlobalLookup()
Definition owneroverlapcopy.hh:497
Dune::BufferedCommunicator BC
Definition owneroverlapcopy.hh:191
OwnerOverlapCopyCommunication(MPI_Comm comm_, SolverCategory::Category cat_=SolverCategory::overlapping, bool freecomm_=false)
Construct the communication without any indices.
Definition owneroverlapcopy.hh:556
ParallelIndexSet & indexSet()
Get the underlying parallel index set.
Definition owneroverlapcopy.hh:482
void dot(const T1 &x, const T1 &y, T2 &result) const
Compute a global dot product of two vectors.
Definition owneroverlapcopy.hh:403
OwnerOverlapCopyCommunication(SolverCategory::Category cat_=SolverCategory::overlapping)
Construct the communication without any indices using MPI_COMM_WORLD.
Definition owneroverlapcopy.hh:574
const CollectiveCommunication< MPI_Comm > & communicator() const
Definition owneroverlapcopy.hh:304
EnumItem< AttributeSet, OwnerOverlapCopyAttributeSet::owner > OwnerSet
Definition owneroverlapcopy.hh:193
void setSolverCategory(SolverCategory set)
Set right Solver Category (default is overlapping).
Definition owneroverlapcopy.hh:292
void copyOwnerToAll(const T &source, T &dest) const
Communicate values from owner data points to all other data points.
Definition owneroverlapcopy.hh:316
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition owneroverlapcopy.hh:473
friend void loadMatrixMarket(M &, const std::string &, OwnerOverlapCopyCommunication< G, L > &, bool)
Load a parallel matrix/vector stored in matrix market format.
Definition matrixmarket.hh:1031
RI::RemoteIndex RX
Definition owneroverlapcopy.hh:190
void addOwnerCopyToAll(const T &source, T &dest) const
Communicate values from owner and copy data points to all other data points and add them to those val...
Definition owneroverlapcopy.hh:367
void freeGlobalLookup()
Definition owneroverlapcopy.hh:522
SolverCategory::Category getSolverCategory() const
Get Solver Category.
Definition owneroverlapcopy.hh:300
Dune::RemoteIndexListModifier< PIS, typename RI::Allocator, false > RILM
Definition owneroverlapcopy.hh:189
Dune::ParallelIndexSet< GlobalIdType, LI, 512 > ParallelIndexSet
The type of the parallel index set.
Definition owneroverlapcopy.hh:451
void buildOwnerToAllInterface() const
Definition owneroverlapcopy.hh:244
gather/scatter callback for communcation
Definition owneroverlapcopy.hh:203
static V gather(const T &a, std::size_t i)
Definition owneroverlapcopy.hh:206
static void scatter(T &a, V v, std::size_t i)
Definition owneroverlapcopy.hh:211
CommPolicy< T >::IndexedType V
Definition owneroverlapcopy.hh:204
CommPolicy< T >::IndexedType V
Definition owneroverlapcopy.hh:219
static V gather(const T &a, std::size_t i)
Definition owneroverlapcopy.hh:221
static void scatter(T &a, V v, std::size_t i)
Definition owneroverlapcopy.hh:226
Categories for the solvers.
Definition solvercategory.hh:18
Category
Definition solvercategory.hh:19
@ overlapping
Category for overlapping solvers.
Definition solvercategory.hh:25