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