dune-istl
2.2.0
|
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set ts=4 sw=2 et sts=2: 00003 // $Id: owneroverlapcopy.hh 1504 2011-10-31 09:26:48Z rebecca $ 00004 #ifndef DUNE_OWNEROVERLAPCOPY_HH 00005 #define DUNE_OWNEROVERLAPCOPY_HH 00006 00007 #include<new> 00008 #include<iostream> 00009 #include<vector> 00010 #include<list> 00011 #include<map> 00012 #include<set> 00013 00014 #include"cmath" 00015 00016 // MPI header 00017 #if HAVE_MPI 00018 #include<mpi.h> 00019 #endif 00020 00021 00022 #include<dune/common/tuples.hh> 00023 #include<dune/common/enumset.hh> 00024 00025 #if HAVE_MPI 00026 #include <dune/common/parallel/indexset.hh> 00027 #include <dune/common/parallel/communicator.hh> 00028 #include <dune/common/parallel/remoteindices.hh> 00029 #include<dune/common/mpicollectivecommunication.hh> 00030 #endif 00031 00032 #include"solvercategory.hh" 00033 #include"istlexception.hh" 00034 #include<dune/common/collectivecommunication.hh> 00035 00036 template<int dim, template<class,class> class Comm> 00037 void testRedistributed(int s); 00038 00039 00040 namespace Dune { 00041 00057 struct OwnerOverlapCopyAttributeSet 00058 { 00059 enum AttributeSet { 00060 owner=1, overlap=2, copy=3 }; 00061 }; 00062 00074 template <class G, class L> 00075 class IndexInfoFromGrid 00076 { 00077 public: 00079 typedef G GlobalIdType; 00080 00082 typedef L LocalIdType; 00083 00090 typedef tuple<GlobalIdType,LocalIdType,int> IndexTripel; 00097 typedef tuple<int,GlobalIdType,int> RemoteIndexTripel; 00098 00104 void addLocalIndex (const IndexTripel& x) 00105 { 00106 if (get<2>(x)!=OwnerOverlapCopyAttributeSet::owner && 00107 get<2>(x)!=OwnerOverlapCopyAttributeSet::overlap && 00108 get<2>(x)!=OwnerOverlapCopyAttributeSet::copy) 00109 DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set"); 00110 localindices.insert(x); 00111 } 00112 00118 void addRemoteIndex (const RemoteIndexTripel& x) 00119 { 00120 if (get<2>(x)!=OwnerOverlapCopyAttributeSet::owner && 00121 get<2>(x)!=OwnerOverlapCopyAttributeSet::overlap && 00122 get<2>(x)!=OwnerOverlapCopyAttributeSet::copy) 00123 DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set"); 00124 remoteindices.insert(x); 00125 } 00126 00131 const std::set<IndexTripel>& localIndices () const 00132 { 00133 return localindices; 00134 } 00135 00140 const std::set<RemoteIndexTripel>& remoteIndices () const 00141 { 00142 return remoteindices; 00143 } 00144 00148 void clear () 00149 { 00150 localindices.clear(); 00151 remoteindices.clear(); 00152 } 00153 00154 private: 00156 std::set<IndexTripel> localindices; 00158 std::set<RemoteIndexTripel> remoteindices; 00159 }; 00160 00161 00162 #if HAVE_MPI 00163 00170 template <class GlobalIdType, class LocalIdType=int> 00171 class OwnerOverlapCopyCommunication 00172 { 00173 template<typename M, typename G, typename L> 00174 friend void loadMatrixMarket(M&, 00175 const std::string&, 00176 OwnerOverlapCopyCommunication<G,L>&, 00177 bool); 00178 // used types 00179 typedef typename IndexInfoFromGrid<GlobalIdType,LocalIdType>::IndexTripel IndexTripel; 00180 typedef typename IndexInfoFromGrid<GlobalIdType,LocalIdType>::RemoteIndexTripel RemoteIndexTripel; 00181 typedef typename std::set<IndexTripel>::const_iterator localindex_iterator; 00182 typedef typename std::set<RemoteIndexTripel>::const_iterator remoteindex_iterator; 00183 typedef typename OwnerOverlapCopyAttributeSet::AttributeSet AttributeSet; 00184 typedef Dune::ParallelLocalIndex<AttributeSet> LI; 00185 public: 00186 typedef Dune::ParallelIndexSet<GlobalIdType,LI,512> PIS; 00187 typedef Dune::RemoteIndices<PIS> RI; 00188 typedef Dune::RemoteIndexListModifier<PIS,typename RI::Allocator,false> RILM; 00189 typedef typename RI::RemoteIndex RX; 00190 typedef Dune::BufferedCommunicator BC; 00191 typedef Dune::Interface IF; 00192 typedef EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner> OwnerSet; 00193 typedef EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy> CopySet; 00194 typedef Combine<EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner>,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::overlap>,AttributeSet> OwnerOverlapSet; 00195 typedef Dune::AllSet<AttributeSet> AllSet; 00196 protected: 00197 00198 00200 template<typename T> 00201 struct CopyGatherScatter 00202 { 00203 typedef typename CommPolicy<T>::IndexedType V; 00204 00205 static V gather(const T& a, std::size_t i) 00206 { 00207 return a[i]; 00208 } 00209 00210 static void scatter(T& a, V v, std::size_t i) 00211 { 00212 a[i] = v; 00213 } 00214 }; 00215 template<typename T> 00216 struct AddGatherScatter 00217 { 00218 typedef typename CommPolicy<T>::IndexedType V; 00219 00220 static V gather(const T& a, std::size_t i) 00221 { 00222 return a[i]; 00223 } 00224 00225 static void scatter(T& a, V v, std::size_t i) 00226 { 00227 a[i] += v; 00228 } 00229 }; 00230 00231 void buildOwnerOverlapToAllInterface () const 00232 { 00233 if (OwnerOverlapToAllInterfaceBuilt) 00234 OwnerOverlapToAllInterface.free(); 00235 typedef Combine<EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner>,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::overlap>,AttributeSet> OwnerOverlapSet; 00236 typedef Combine<OwnerOverlapSet,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy>,AttributeSet> AllSet; 00237 OwnerOverlapSet sourceFlags; 00238 AllSet destFlags; 00239 OwnerOverlapToAllInterface.build(ri,sourceFlags,destFlags); 00240 OwnerOverlapToAllInterfaceBuilt = true; 00241 } 00242 00243 void buildOwnerToAllInterface () const 00244 { 00245 if (OwnerToAllInterfaceBuilt) 00246 OwnerToAllInterface.free(); 00247 OwnerSet sourceFlags; 00248 AllSet destFlags; 00249 OwnerToAllInterface.build(ri,sourceFlags,destFlags); 00250 OwnerToAllInterfaceBuilt = true; 00251 } 00252 00253 void buildOwnerCopyToAllInterface () const 00254 { 00255 if (OwnerCopyToAllInterfaceBuilt) 00256 OwnerCopyToAllInterface.free(); 00257 typedef Combine<EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner>,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy>,AttributeSet> OwnerCopySet; 00258 typedef Combine<OwnerCopySet,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::overlap>,AttributeSet> AllSet; 00259 OwnerCopySet sourceFlags; 00260 AllSet destFlags; 00261 OwnerCopyToAllInterface.build(ri,sourceFlags,destFlags); 00262 OwnerCopyToAllInterfaceBuilt = true; 00263 } 00264 00265 void buildOwnerCopyToOwnerCopyInterface () const 00266 { 00267 if (OwnerCopyToOwnerCopyInterfaceBuilt) 00268 OwnerCopyToOwnerCopyInterface.free(); 00269 typedef Combine<EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner>,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy>,AttributeSet> OwnerCopySet; 00270 OwnerCopySet sourceFlags; 00271 OwnerCopySet destFlags; 00272 OwnerCopyToOwnerCopyInterface.build(ri,sourceFlags,destFlags); 00273 OwnerCopyToOwnerCopyInterfaceBuilt = true; 00274 } 00275 00276 void buildCopyToAllInterface () const 00277 { 00278 if (CopyToAllInterfaceBuilt) 00279 CopyToAllInterface.free(); 00280 CopySet sourceFlags; 00281 AllSet destFlags; 00282 CopyToAllInterface.build(ri,sourceFlags,destFlags); 00283 CopyToAllInterfaceBuilt = true; 00284 } 00285 00286 public: 00287 00292 void setSolverCategory (SolverCategory set) { 00293 category = set; 00294 } 00295 00301 const SolverCategory::Category getSolverCategory () const { 00302 return category; 00303 } 00304 00305 const CollectiveCommunication<MPI_Comm>& communicator() const 00306 { 00307 return cc; 00308 } 00309 00316 template<class T> 00317 void copyOwnerToAll (const T& source, T& dest) const 00318 { 00319 if (!OwnerToAllInterfaceBuilt) 00320 buildOwnerToAllInterface (); 00321 BC communicator; 00322 communicator.template build<T>(OwnerToAllInterface); 00323 communicator.template forward<CopyGatherScatter<T> >(source,dest); 00324 communicator.free(); 00325 } 00326 00333 template<class T> 00334 void copyCopyToAll (const T& source, T& dest) const 00335 { 00336 if (!CopyToAllInterfaceBuilt) 00337 buildCopyToAllInterface (); 00338 BC communicator; 00339 communicator.template build<T>(CopyToAllInterface); 00340 communicator.template forward<CopyGatherScatter<T> >(source,dest); 00341 communicator.free(); 00342 } 00343 00350 template<class T> 00351 void addOwnerOverlapToAll (const T& source, T& dest) const 00352 { 00353 if (!OwnerOverlapToAllInterfaceBuilt) 00354 buildOwnerOverlapToAllInterface (); 00355 BC communicator; 00356 communicator.template build<T>(OwnerOverlapToAllInterface); 00357 communicator.template forward<AddGatherScatter<T> >(source,dest); 00358 communicator.free(); 00359 } 00360 00367 template<class T> 00368 void addOwnerCopyToAll (const T& source, T& dest) const 00369 { 00370 if (!OwnerCopyToAllInterfaceBuilt) 00371 buildOwnerCopyToAllInterface (); 00372 BC communicator; 00373 communicator.template build<T>(OwnerCopyToAllInterface); 00374 communicator.template forward<AddGatherScatter<T> >(source,dest); 00375 communicator.free(); 00376 } 00377 00384 template<class T> 00385 void addOwnerCopyToOwnerCopy (const T& source, T& dest) const 00386 { 00387 if (!OwnerCopyToOwnerCopyInterfaceBuilt) 00388 buildOwnerCopyToOwnerCopyInterface (); 00389 BC communicator; 00390 communicator.template build<T>(OwnerCopyToOwnerCopyInterface); 00391 communicator.template forward<AddGatherScatter<T> >(source,dest); 00392 communicator.free(); 00393 } 00394 00395 00403 template<class T1, class T2> 00404 void dot (const T1& x, const T1& y, T2& result) const 00405 { 00406 // set up mask vector 00407 if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size())) 00408 { 00409 mask.resize(x.size()); 00410 for (typename std::vector<double>::size_type i=0; i<mask.size(); i++) 00411 mask[i] = 1; 00412 for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i) 00413 if (i->local().attribute()!=OwnerOverlapCopyAttributeSet::owner) 00414 mask[i->local().local()] = 0; 00415 } 00416 result = 0; 00417 00418 for (typename T1::size_type i=0; i<x.size(); i++) 00419 result += x[i]*(y[i])*mask[i]; 00420 result = cc.sum(result); 00421 return; 00422 } 00423 00430 template<class T1> 00431 double norm (const T1& x) const 00432 { 00433 // set up mask vector 00434 if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size())) 00435 { 00436 mask.resize(x.size()); 00437 for (typename std::vector<double>::size_type i=0; i<mask.size(); i++) 00438 mask[i] = 1; 00439 for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i) 00440 if (i->local().attribute()!=OwnerOverlapCopyAttributeSet::owner) 00441 mask[i->local().local()] = 0; 00442 } 00443 double result = 0; 00444 for (typename T1::size_type i=0; i<x.size(); i++) 00445 result += x[i].two_norm2()*mask[i]; 00446 return sqrt(cc.sum(result)); 00447 } 00448 00449 typedef Dune::EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy> CopyFlags; 00450 00452 typedef Dune::ParallelIndexSet<GlobalIdType,LI,512> ParallelIndexSet; 00453 00455 typedef Dune::RemoteIndices<PIS> RemoteIndices; 00456 00459 typedef Dune::GlobalLookupIndexSet<ParallelIndexSet> GlobalLookupIndexSet; 00460 00465 const ParallelIndexSet& indexSet() const 00466 { 00467 return pis; 00468 } 00469 00474 const RemoteIndices& remoteIndices() const 00475 { 00476 return ri; 00477 } 00478 00483 ParallelIndexSet& indexSet() 00484 { 00485 return pis; 00486 } 00487 00488 00493 RemoteIndices& remoteIndices() 00494 { 00495 return ri; 00496 } 00497 00498 void buildGlobalLookup() 00499 { 00500 if(globalLookup_){ 00501 if(pis.seqNo()==oldseqNo) 00502 // Nothing changed! 00503 return; 00504 delete globalLookup_; 00505 } 00506 00507 globalLookup_ = new GlobalLookupIndexSet(pis); 00508 oldseqNo = pis.seqNo(); 00509 } 00510 00511 void buildGlobalLookup(std::size_t size) 00512 { 00513 if(globalLookup_){ 00514 if(pis.seqNo()==oldseqNo) 00515 // Nothing changed! 00516 return; 00517 delete globalLookup_; 00518 } 00519 globalLookup_ = new GlobalLookupIndexSet(pis, size); 00520 oldseqNo = pis.seqNo(); 00521 } 00522 00523 void freeGlobalLookup() 00524 { 00525 delete globalLookup_; 00526 globalLookup_=0; 00527 } 00528 00529 const GlobalLookupIndexSet& globalLookup() const 00530 { 00531 assert(globalLookup_ != 0); 00532 return *globalLookup_; 00533 } 00534 00540 template<class T1> 00541 void project (T1& x) const 00542 { 00543 for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i) 00544 if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy) 00545 x[i->local().local()] = 0; 00546 } 00547 00557 OwnerOverlapCopyCommunication (MPI_Comm comm_, 00558 SolverCategory::Category cat_ = SolverCategory::overlapping, 00559 bool freecomm_ = false) 00560 : comm(comm_), cc(comm_), pis(), ri(pis,pis,comm_), 00561 OwnerToAllInterfaceBuilt(false), OwnerOverlapToAllInterfaceBuilt(false), 00562 OwnerCopyToAllInterfaceBuilt(false), OwnerCopyToOwnerCopyInterfaceBuilt(false), 00563 CopyToAllInterfaceBuilt(false), globalLookup_(0), category(cat_), 00564 freecomm(freecomm_) 00565 {} 00566 00575 OwnerOverlapCopyCommunication (SolverCategory::Category cat_ = SolverCategory::overlapping) 00576 : comm(MPI_COMM_WORLD), cc(MPI_COMM_WORLD), pis(), ri(pis,pis,MPI_COMM_WORLD), 00577 OwnerToAllInterfaceBuilt(false), OwnerOverlapToAllInterfaceBuilt(false), 00578 OwnerCopyToAllInterfaceBuilt(false), OwnerCopyToOwnerCopyInterfaceBuilt(false), 00579 CopyToAllInterfaceBuilt(false), globalLookup_(0), category(cat_), freecomm(false) 00580 {} 00581 00589 OwnerOverlapCopyCommunication (const IndexInfoFromGrid<GlobalIdType, LocalIdType>& indexinfo, 00590 MPI_Comm comm_, 00591 SolverCategory::Category cat_ = SolverCategory::overlapping, 00592 bool freecomm_ = false) 00593 : comm(comm_), cc(comm_), OwnerToAllInterfaceBuilt(false), 00594 OwnerOverlapToAllInterfaceBuilt(false), OwnerCopyToAllInterfaceBuilt(false), 00595 OwnerCopyToOwnerCopyInterfaceBuilt(false), CopyToAllInterfaceBuilt(false), 00596 globalLookup_(0), category(cat_), freecomm(freecomm_) 00597 { 00598 // set up an ISTL index set 00599 pis.beginResize(); 00600 for (localindex_iterator i=indexinfo.localIndices().begin(); i!=indexinfo.localIndices().end(); ++i) 00601 { 00602 if (get<2>(*i)==OwnerOverlapCopyAttributeSet::owner) 00603 pis.add(get<0>(*i),LI(get<1>(*i),OwnerOverlapCopyAttributeSet::owner,true)); 00604 if (get<2>(*i)==OwnerOverlapCopyAttributeSet::overlap) 00605 pis.add(get<0>(*i),LI(get<1>(*i),OwnerOverlapCopyAttributeSet::overlap,true)); 00606 if (get<2>(*i)==OwnerOverlapCopyAttributeSet::copy) 00607 pis.add(get<0>(*i),LI(get<1>(*i),OwnerOverlapCopyAttributeSet::copy,true)); 00608 // std::cout << cc.rank() << ": adding index " << get<0>(*i) << " " << get<1>(*i) << " " << get<2>(*i) << std::endl; 00609 } 00610 pis.endResize(); 00611 00612 // build remote indices WITHOUT communication 00613 // std::cout << cc.rank() << ": build remote indices" << std::endl; 00614 ri.setIndexSets(pis,pis,cc); 00615 if (indexinfo.remoteIndices().size()>0) 00616 { 00617 remoteindex_iterator i=indexinfo.remoteIndices().begin(); 00618 int p = get<0>(*i); 00619 RILM modifier = ri.template getModifier<false,true>(p); 00620 typename PIS::const_iterator pi=pis.begin(); 00621 for ( ; i!=indexinfo.remoteIndices().end(); ++i) 00622 { 00623 // handle processor change 00624 if (p!=get<0>(*i)) 00625 { 00626 p = get<0>(*i); 00627 modifier = ri.template getModifier<false,true>(p); 00628 pi=pis.begin(); 00629 } 00630 00631 // position to correct entry in parallel index set 00632 while (pi->global()!=get<1>(*i) && pi!=pis.end()) 00633 ++pi; 00634 if (pi==pis.end()) 00635 DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set"); 00636 00637 // insert entry 00638 // std::cout << cc.rank() << ": adding remote index " << get<0>(*i) << " " << get<1>(*i) << " " << get<2>(*i) << std::endl; 00639 if (get<2>(*i)==OwnerOverlapCopyAttributeSet::owner) 00640 modifier.insert(RX(OwnerOverlapCopyAttributeSet::owner,&(*pi))); 00641 if (get<2>(*i)==OwnerOverlapCopyAttributeSet::overlap) 00642 modifier.insert(RX(OwnerOverlapCopyAttributeSet::overlap,&(*pi))); 00643 if (get<2>(*i)==OwnerOverlapCopyAttributeSet::copy) 00644 modifier.insert(RX(OwnerOverlapCopyAttributeSet::copy,&(*pi))); 00645 } 00646 }else{ 00647 // Force remote indices to be synced! 00648 ri.template getModifier<false,true>(0); 00649 } 00650 } 00651 00652 // destructor: free memory in some objects 00653 ~OwnerOverlapCopyCommunication () 00654 { 00655 ri.free(); 00656 if (OwnerToAllInterfaceBuilt) OwnerToAllInterface.free(); 00657 if (OwnerOverlapToAllInterfaceBuilt) OwnerOverlapToAllInterface.free(); 00658 if (OwnerCopyToAllInterfaceBuilt) OwnerCopyToAllInterface.free(); 00659 if (OwnerCopyToOwnerCopyInterfaceBuilt) OwnerCopyToOwnerCopyInterface.free(); 00660 if (CopyToAllInterfaceBuilt) CopyToAllInterface.free(); 00661 if (globalLookup_) delete globalLookup_; 00662 if (freecomm==true) 00663 if(comm!=MPI_COMM_NULL) 00664 MPI_Comm_free(&comm); 00665 } 00666 00667 private: 00668 OwnerOverlapCopyCommunication (const OwnerOverlapCopyCommunication&) 00669 {} 00670 MPI_Comm comm; 00671 CollectiveCommunication<MPI_Comm> cc; 00672 PIS pis; 00673 RI ri; 00674 mutable IF OwnerToAllInterface; 00675 mutable bool OwnerToAllInterfaceBuilt; 00676 mutable IF OwnerOverlapToAllInterface; 00677 mutable bool OwnerOverlapToAllInterfaceBuilt; 00678 mutable IF OwnerCopyToAllInterface; 00679 mutable bool OwnerCopyToAllInterfaceBuilt; 00680 mutable IF OwnerCopyToOwnerCopyInterface; 00681 mutable bool OwnerCopyToOwnerCopyInterfaceBuilt; 00682 mutable IF CopyToAllInterface; 00683 mutable bool CopyToAllInterfaceBuilt; 00684 mutable std::vector<double> mask; 00685 int oldseqNo; 00686 GlobalLookupIndexSet* globalLookup_; 00687 SolverCategory::Category category; 00688 bool freecomm; 00689 }; 00690 00691 #endif 00692 00693 00696 } // end namespace 00697 00698 #endif