dune-istl  2.2.0
transfer.hh
Go to the documentation of this file.
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set ts=8 sw=2 et sts=2:
00003 #ifndef DUNE_AMGTRANSFER_HH
00004 #define DUNE_AMGTRANSFER_HH
00005 
00006 #include<dune/istl/bvector.hh>
00007 #include<dune/istl/matrixredistribute.hh>
00008 #include<dune/istl/paamg/pinfo.hh>
00009 #include<dune/istl/owneroverlapcopy.hh>
00010 #include<dune/istl/paamg/aggregates.hh>
00011 #include<dune/common/exceptions.hh>
00012 
00013 namespace Dune
00014 {
00015   namespace Amg
00016   {
00017     
00028     template<class V1, class V2, class T>
00029     class Transfer
00030     {
00031     
00032     public:
00033       typedef V1 Vertex;
00034       typedef V2 Vector;
00035     
00036       template<typename T1, typename R>
00037       static void prolongate(const AggregatesMap<Vertex>& aggregates, Vector& coarse, Vector& fine, 
00038                              Vector& fineRedist,T1 damp, R& redistributor=R());
00039 
00040       template<typename T1, typename R>
00041       static void prolongate(const AggregatesMap<Vertex>& aggregates, Vector& coarse, Vector& fine, 
00042                              T1 damp);
00043 
00044       static void restrict(const AggregatesMap<Vertex>& aggregates, Vector& coarse, const Vector& fine,
00045                            T& comm);
00046     };
00047 
00048     template<class V,class V1>
00049     class Transfer<V,V1, SequentialInformation>
00050     {
00051     public:
00052       typedef V Vertex;
00053       typedef V1 Vector;
00054       typedef RedistributeInformation<SequentialInformation> Redist;
00055       template<typename T1>
00056       static void prolongate(const AggregatesMap<Vertex>& aggregates, Vector& coarse, Vector& fine, 
00057                              Vector& fineRedist, T1 damp,
00058                              const SequentialInformation& comm=SequentialInformation(),
00059                              const Redist& redist=Redist());
00060       template<typename T1>
00061       static void prolongate(const AggregatesMap<Vertex>& aggregates, Vector& coarse, Vector& fine, 
00062                              T1 damp,
00063                              const SequentialInformation& comm=SequentialInformation());
00064 
00065 
00066       static void restrict(const AggregatesMap<Vertex>& aggregates, Vector& coarse, const Vector& fine,
00067                            const SequentialInformation& comm);
00068     };
00069 
00070 #if HAVE_MPI
00071 
00072     template<class V,class V1, class T1, class T2>
00073     class Transfer<V,V1,OwnerOverlapCopyCommunication<T1,T2> > 
00074     {
00075     public:
00076       typedef V Vertex;
00077       typedef V1 Vector;
00078       typedef RedistributeInformation<OwnerOverlapCopyCommunication<T1,T2> > Redist;
00079       template<typename T3>
00080       static void prolongate(const AggregatesMap<Vertex>& aggregates, Vector& coarse, Vector& fine, 
00081                              Vector& fineRedist, T3 damp, OwnerOverlapCopyCommunication<T1,T2>& comm,
00082                              const Redist& redist);
00083       template<typename T3>
00084       static void prolongate(const AggregatesMap<Vertex>& aggregates, Vector& coarse, Vector& fine, 
00085                              T3 damp, OwnerOverlapCopyCommunication<T1,T2>& comm);
00086 
00087       static void restrict(const AggregatesMap<Vertex>& aggregates, Vector& coarse, const Vector& fine,
00088                            OwnerOverlapCopyCommunication<T1,T2>& comm);
00089     };
00090 
00091 #endif
00092 
00093     template<class V, class V1>
00094     template<typename T>
00095     inline void
00096     Transfer<V,V1,SequentialInformation>::prolongate(const AggregatesMap<Vertex>& aggregates,
00097                                                      Vector& coarse, Vector& fine, Vector& fineRedist,
00098                                                      T damp,
00099                                                      const SequentialInformation& comm,
00100                                                      const Redist& redist)
00101     {
00102       prolongate(aggregates, coarse, fine, damp);
00103     }
00104     template<class V, class V1>
00105      template<typename T>
00106     inline void 
00107     Transfer<V,V1,SequentialInformation>::prolongate(const AggregatesMap<Vertex>& aggregates,
00108                                                      Vector& coarse, Vector& fine,
00109                                                      T damp,
00110                                                      const SequentialInformation& comm)
00111     {
00112       typedef typename Vector::iterator Iterator;
00113 
00114       Iterator end = coarse.end();
00115       Iterator begin= coarse.begin();
00116       for(;begin!=end;++begin)
00117         *begin*=damp;
00118       end=fine.end();
00119       begin=fine.begin();
00120       
00121       for(Iterator block=begin; block != end; ++block){
00122         std::ptrdiff_t index=block-begin;
00123         const Vertex& vertex = aggregates[index];
00124         if(vertex != AggregatesMap<Vertex>::ISOLATED)
00125           *block += coarse[aggregates[index]];
00126       }
00127     }
00128 
00129     template<class V, class V1>
00130     inline void 
00131     Transfer<V,V1,SequentialInformation>::restrict(const AggregatesMap<Vertex>& aggregates,
00132                                                    Vector& coarse, 
00133                                                    const Vector& fine,
00134                                                    const SequentialInformation& comm)
00135     {
00136       // Set coarse vector to zero
00137       coarse=0;
00138       
00139       typedef typename Vector::const_iterator Iterator;
00140       Iterator end = fine.end();
00141       Iterator begin=fine.begin();
00142 
00143       for(Iterator block=begin; block != end; ++block){
00144         const Vertex& vertex = aggregates[block-begin];
00145         if(vertex != AggregatesMap<Vertex>::ISOLATED)
00146           coarse[vertex] += *block;
00147       }
00148     }
00149 
00150 #if HAVE_MPI
00151     template<class V, class V1, class T1, class T2>
00152     template<typename T3>
00153     inline void Transfer<V,V1,OwnerOverlapCopyCommunication<T1,T2> >::prolongate(const AggregatesMap<Vertex>& aggregates,
00154                                                                                 Vector& coarse, Vector& fine, 
00155                                                                                 Vector& fineRedist, T3 damp, 
00156                                                                                 OwnerOverlapCopyCommunication<T1,T2>& comm,
00157                                                                                 const Redist& redist)
00158     {
00159       if(fineRedist.size()>0)
00160         // we operated on the coarse level
00161         Transfer<V,V1,SequentialInformation>::prolongate(aggregates, coarse, fineRedist, damp);
00162       
00163       // TODO This could be accomplished with one communication, too!
00164       redist.redistributeBackward(fine, fineRedist);
00165       comm.copyOwnerToAll(fine,fine);
00166     }
00167     
00168     template<class V, class V1, class T1, class T2>
00169     template<typename T3>
00170     inline void Transfer<V,V1,OwnerOverlapCopyCommunication<T1,T2> >::prolongate(const AggregatesMap<Vertex>& aggregates,
00171                                                                                 Vector& coarse, Vector& fine, 
00172                                                                                 T3 damp, 
00173                                                                                 OwnerOverlapCopyCommunication<T1,T2>& comm)
00174     {
00175       Transfer<V,V1,SequentialInformation>::prolongate(aggregates, coarse, fine, damp);
00176     }
00177     template<class V, class V1, class T1, class T2>
00178     inline void Transfer<V,V1,OwnerOverlapCopyCommunication<T1,T2> >::restrict(const AggregatesMap<Vertex>& aggregates,
00179                                                                                Vector& coarse, const Vector& fine, 
00180                                                                                OwnerOverlapCopyCommunication<T1,T2>& comm)
00181     {
00182       Transfer<V,V1,SequentialInformation>::restrict(aggregates, coarse, fine, SequentialInformation());
00183       // We need this here to avoid it in the smoothers on the coarse level.
00184       // There (in the preconditioner d is const.
00185       comm.project(coarse);
00186     }
00187 #endif
00188 
00189       }// namspace Amg
00190     } // namspace Dune
00191 #endif