dune-istl
2.2.0
|
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