dune-istl  2.2.0
owneroverlapcopy.hh
Go to the documentation of this file.
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