dune-istl
2.2.0
|
00001 #ifndef DUNE_AMG_KAMG_HH 00002 #define DUNE_AMG_KAMG_HH 00003 00004 #include<dune/istl/preconditioners.hh> 00005 #include"amg.hh" 00006 00007 namespace Dune 00008 { 00009 namespace Amg 00010 { 00011 00026 template<class AMG> 00027 class KAmgTwoGrid 00028 : public Preconditioner<typename AMG::Domain,typename AMG::Range> 00029 { 00031 typedef typename AMG::Domain Domain; 00033 typedef typename AMG::Range Range; 00034 public: 00035 00036 enum { 00038 category = AMG::category 00039 }; 00040 00048 KAmgTwoGrid(AMG& amg, InverseOperator<Domain,Range>* coarseSolver) 00049 : amg_(amg), coarseSolver_(coarseSolver) 00050 {} 00051 00053 void pre(typename AMG::Domain& x, typename AMG::Range& b) 00054 {} 00055 00057 void post(typename AMG::Domain& x) 00058 {} 00059 00061 void apply(typename AMG::Domain& v, const typename AMG::Range& d) 00062 { 00063 *amg_.rhs = d; 00064 *amg_.lhs = v; 00065 // Copy data 00066 *amg_.update=0; 00067 00068 amg_.presmooth(); 00069 bool processFineLevel = 00070 amg_.moveToCoarseLevel(); 00071 00072 if(processFineLevel){ 00073 typename AMG::Range b=*amg_.rhs; 00074 typename AMG::Domain x=*amg_.update; 00075 InverseOperatorResult res; 00076 coarseSolver_->apply(x, b, res); 00077 *amg_.update=x; 00078 } 00079 00080 amg_.moveToFineLevel(processFineLevel); 00081 00082 amg_.postsmooth(); 00083 v=*amg_.update; 00084 } 00085 00090 InverseOperator<Domain,Range>* coarseSolver() 00091 { 00092 return coarseSolver_; 00093 } 00094 00096 ~KAmgTwoGrid() 00097 {} 00098 00099 private: 00101 AMG& amg_; 00103 InverseOperator<Domain,Range>* coarseSolver_; 00104 }; 00105 00106 00107 00118 template<class M, class X, class S, class PI=SequentialInformation, 00119 class K=BiCGSTABSolver<X>, class A=std::allocator<X> > 00120 class KAMG : public Preconditioner<X,X> 00121 { 00122 public: 00124 typedef AMG<M,X,S,PI,A> Amg; 00126 typedef K KrylovSolver; 00128 typedef typename Amg::OperatorHierarchy OperatorHierarchy; 00130 typedef typename Amg::CoarseSolver CoarseSolver; 00132 typedef typename Amg::ParallelInformation ParallelInformation; 00134 typedef typename Amg::SmootherArgs SmootherArgs; 00136 typedef typename Amg::Operator Operator; 00138 typedef typename Amg::Domain Domain; 00140 typedef typename Amg::Range Range; 00142 typedef typename Amg::ParallelInformationHierarchy ParallelInformationHierarchy; 00144 typedef typename Amg::ScalarProduct ScalarProduct; 00145 00146 enum { 00148 category = Amg::category 00149 }; 00163 KAMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, 00164 const SmootherArgs& smootherArgs, std::size_t gamma, 00165 std::size_t preSmoothingSteps =1, std::size_t postSmoothingSteps = 1, 00166 std::size_t maxLevelKrylovSteps = 3 , double minDefectReduction =1e-1); 00167 00184 template<class C> 00185 KAMG(const Operator& fineOperator, const C& criterion, 00186 const SmootherArgs& smootherArgs, std::size_t gamma=1, 00187 std::size_t preSmoothingSteps=1, std::size_t postSmoothingSteps=1, 00188 std::size_t maxLevelKrylovSteps=3, double minDefectReduction=1e-1, 00189 const ParallelInformation& pinfo=ParallelInformation()); 00190 00192 void pre(Domain& x, Range& b); 00194 void post(Domain& x); 00196 void apply(Domain& x, const Range& b); 00197 00198 std::size_t maxlevels(); 00199 00200 private: 00202 Amg amg; 00203 00205 std::size_t maxLevelKrylovSteps; 00206 00208 double levelDefectReduction; 00209 00211 std::vector<typename Amg::ScalarProduct*> scalarproducts; 00212 00214 std::vector<KAmgTwoGrid<Amg>*> ksolvers; 00215 }; 00216 00217 template<class M, class X, class S, class P, class K, class A> 00218 KAMG<M,X,S,P,K,A>::KAMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, 00219 const SmootherArgs& smootherArgs, 00220 std::size_t gamma, std::size_t preSmoothingSteps, 00221 std::size_t postSmoothingSteps, 00222 std::size_t ksteps, double reduction) 00223 : amg(matrices, coarseSolver, smootherArgs, gamma, preSmoothingSteps, 00224 postSmoothingSteps), maxLevelKrylovSteps(ksteps), levelDefectReduction(reduction) 00225 {} 00226 00227 template<class M, class X, class S, class P, class K, class A> 00228 template<class C> 00229 KAMG<M,X,S,P,K,A>::KAMG(const Operator& fineOperator, const C& criterion, 00230 const SmootherArgs& smootherArgs, std::size_t gamma, 00231 std::size_t preSmoothingSteps, std::size_t postSmoothingSteps, 00232 std::size_t ksteps, double reduction, 00233 const ParallelInformation& pinfo) 00234 : amg(fineOperator, criterion, smootherArgs, gamma, preSmoothingSteps, 00235 postSmoothingSteps, false, pinfo), maxLevelKrylovSteps(ksteps), levelDefectReduction(reduction) 00236 {} 00237 00238 00239 template<class M, class X, class S, class P, class K, class A> 00240 void KAMG<M,X,S,P,K,A>::pre(Domain& x, Range& b) 00241 { 00242 amg.pre(x,b); 00243 scalarproducts.reserve(amg.levels()); 00244 ksolvers.reserve(amg.levels()); 00245 00246 typename OperatorHierarchy::ParallelMatrixHierarchy::Iterator 00247 matrix = amg.matrices_->matrices().coarsest(); 00248 typename ParallelInformationHierarchy::Iterator 00249 pinfo = amg.matrices_->parallelInformation().coarsest(); 00250 bool hasCoarsest=(amg.levels()==amg.maxlevels()); 00251 KrylovSolver* ks=0; 00252 00253 if(hasCoarsest){ 00254 if(matrix==amg.matrices_->matrices().finest()) 00255 return; 00256 --matrix; 00257 --pinfo; 00258 ksolvers.push_back(new KAmgTwoGrid<Amg>(amg, amg.solver_)); 00259 }else 00260 ksolvers.push_back(new KAmgTwoGrid<Amg>(amg, ks)); 00261 00262 std::ostringstream s; 00263 00264 if(matrix!=amg.matrices_->matrices().finest()) 00265 while(true){ 00266 scalarproducts.push_back(Amg::ScalarProductChooser::construct(*pinfo)); 00267 ks = new KrylovSolver(*matrix, *(scalarproducts.back()), 00268 *(ksolvers.back()), levelDefectReduction, 00269 maxLevelKrylovSteps, 0); 00270 ksolvers.push_back(new KAmgTwoGrid<Amg>(amg, ks)); 00271 --matrix; 00272 --pinfo; 00273 if(matrix==amg.matrices_->matrices().finest()) 00274 break; 00275 } 00276 } 00277 00278 00279 template<class M, class X, class S, class P, class K, class A> 00280 void KAMG<M,X,S,P,K,A>::post(Domain& x) 00281 { 00282 typedef typename std::vector<KAmgTwoGrid<Amg>*>::reverse_iterator KIterator; 00283 typedef typename std::vector<ScalarProduct*>::iterator SIterator; 00284 for(KIterator kiter = ksolvers.rbegin(); kiter != ksolvers.rend(); 00285 ++kiter) 00286 { 00287 if((*kiter)->coarseSolver()!=amg.solver_) 00288 delete (*kiter)->coarseSolver(); 00289 delete *kiter; 00290 } 00291 for(SIterator siter = scalarproducts.begin(); siter!=scalarproducts.end(); 00292 ++siter) 00293 delete *siter; 00294 amg.post(x); 00295 00296 } 00297 00298 template<class M, class X, class S, class P, class K, class A> 00299 void KAMG<M,X,S,P,K,A>::apply(Domain& x, const Range& b) 00300 { 00301 amg.initIteratorsWithFineLevel(); 00302 if(ksolvers.size()==0) 00303 { 00304 Range tb=b; 00305 InverseOperatorResult res; 00306 amg.solver_->apply(x,tb,res); 00307 }else 00308 ksolvers.back()->apply(x,b); 00309 } 00310 00311 template<class M, class X, class S, class P, class K, class A> 00312 std::size_t KAMG<M,X,S,P,K,A>::maxlevels() 00313 { 00314 return amg.maxlevels(); 00315 } 00316 00318 } // Amg 00319 } // Dune 00320 00321 #endif