dune-istl  2.2.0
amg.hh
Go to the documentation of this file.
00001 // $Id: amg.hh 1547 2012-04-19 17:19:40Z mblatt $
00002 #ifndef DUNE_AMG_AMG_HH
00003 #define DUNE_AMG_AMG_HH
00004 
00005 #include<memory>
00006 #include<dune/common/exceptions.hh>
00007 #include<dune/istl/paamg/smoother.hh>
00008 #include<dune/istl/paamg/transfer.hh>
00009 #include<dune/istl/paamg/hierarchy.hh>
00010 #include<dune/istl/solvers.hh>
00011 #include<dune/istl/scalarproducts.hh>
00012 #include<dune/istl/superlu.hh>
00013 #include<dune/istl/solvertype.hh>
00014 #include<dune/common/typetraits.hh>
00015 #include<dune/common/exceptions.hh>
00016 
00017 namespace Dune
00018 {
00019   namespace Amg
00020   {
00038     template<class M, class X, class S, class P, class K, class A>
00039     class KAMG;
00040     
00041     template<class T>
00042     class KAmgTwoGrid;
00043     
00051     template<class M, class X, class S, class PI=SequentialInformation,
00052              class A=std::allocator<X> >
00053     class AMG : public Preconditioner<X,X>
00054     {
00055       template<class M1, class X1, class S1, class P1, class K1, class A1>
00056       friend class KAMG;
00057       
00058       friend class KAmgTwoGrid<AMG>;
00059       
00060     public:
00062       typedef M Operator;
00069       typedef PI ParallelInformation;
00071       typedef MatrixHierarchy<M, ParallelInformation, A> OperatorHierarchy;
00073       typedef typename OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy;
00074 
00076       typedef X Domain;
00078       typedef X Range;
00080       typedef InverseOperator<X,X> CoarseSolver;
00086       typedef S Smoother;
00087   
00089       typedef typename SmootherTraits<Smoother>::Arguments SmootherArgs;
00090       
00091       enum {
00093         category = S::category
00094       };
00095 
00111       AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, 
00112           const SmootherArgs& smootherArgs, std::size_t gamma,
00113           std::size_t preSmoothingSteps,
00114           std::size_t postSmoothingSteps, 
00115           bool additive=false) DUNE_DEPRECATED;
00116 
00126       AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, 
00127           const SmootherArgs& smootherArgs, const Parameters& parms);
00128 
00147       template<class C>
00148       AMG(const Operator& fineOperator, const C& criterion,
00149           const SmootherArgs& smootherArgs, std::size_t gamma,
00150           std::size_t preSmoothingSteps, 
00151           std::size_t postSmoothingSteps,
00152           bool additive=false, 
00153           const ParallelInformation& pinfo=ParallelInformation()) DUNE_DEPRECATED;
00154 
00166       template<class C>
00167       AMG(const Operator& fineOperator, const C& criterion,
00168           const SmootherArgs& smootherArgs,
00169           const ParallelInformation& pinfo=ParallelInformation());
00170 
00171       ~AMG();
00172 
00174       void pre(Domain& x, Range& b);
00175 
00177       void apply(Domain& v, const Range& d);
00178       
00180       void post(Domain& x);
00181 
00186       template<class A1>
00187       void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
00188       
00189       std::size_t levels();
00190       
00191       std::size_t maxlevels();
00192 
00201       void recalculateHierarchy()
00202       {
00203         matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
00204       }
00205 
00210       bool usesDirectCoarseLevelSolver() const;
00211       
00212     private:
00214       void mgc();
00215 
00216       typename Hierarchy<Smoother,A>::Iterator smoother;
00217       typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix;
00218       typename ParallelInformationHierarchy::Iterator pinfo;
00219       typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
00220       typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
00221       typename Hierarchy<Domain,A>::Iterator lhs;
00222       typename Hierarchy<Domain,A>::Iterator update;
00223       typename Hierarchy<Range,A>::Iterator rhs;
00224       
00225       void additiveMgc();
00226 
00228       void presmooth();
00229 
00231       void postsmooth();
00232       
00236       void moveToFineLevel(bool processedFineLevel);
00237 
00239       bool moveToCoarseLevel();
00240 
00242       void initIteratorsWithFineLevel();
00243 
00245       OperatorHierarchy* matrices_;
00247       SmootherArgs smootherArgs_;
00249       Hierarchy<Smoother,A> smoothers_;
00251       CoarseSolver* solver_;
00253       Hierarchy<Range,A>* rhs_;
00255       Hierarchy<Domain,A>* lhs_;
00257       Hierarchy<Domain,A>* update_;
00259       typedef Dune::ScalarProductChooser<X,PI,M::category> ScalarProductChooser;
00261       typedef typename ScalarProductChooser::ScalarProduct ScalarProduct;
00263       ScalarProduct* scalarProduct_;
00265       std::size_t gamma_;
00267       std::size_t preSteps_;
00269       std::size_t postSteps_;
00270       std::size_t level;
00271       bool buildHierarchy_;
00272       bool additive;
00273       bool coarsesolverconverged;
00274       Smoother *coarseSmoother_;
00276       std::size_t verbosity_;
00277     };
00278 
00279     template<class M, class X, class S, class PI, class A>
00280     AMG<M,X,S,PI,A>::AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, 
00281                         const SmootherArgs& smootherArgs,
00282                         std::size_t gamma, std::size_t preSmoothingSteps,
00283                         std::size_t postSmoothingSteps, bool additive_)
00284       : matrices_(&matrices), smootherArgs_(smootherArgs),
00285         smoothers_(), solver_(&coarseSolver), scalarProduct_(0),
00286         gamma_(gamma), preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(false),
00287         additive(additive_), coarsesolverconverged(true),
00288         coarseSmoother_(), verbosity_(2)
00289     {
00290       assert(matrices_->isBuilt());
00291       
00292       // build the necessary smoother hierarchies
00293       matrices_->coarsenSmoother(smoothers_, smootherArgs_);
00294     }
00295 
00296     template<class M, class X, class S, class PI, class A>
00297     AMG<M,X,S,PI,A>::AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, 
00298                          const SmootherArgs& smootherArgs,
00299                          const Parameters& parms)
00300       : matrices_(&matrices), smootherArgs_(smootherArgs),
00301         smoothers_(), solver_(&coarseSolver), scalarProduct_(0),
00302         gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()), 
00303         postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
00304         additive(parms.getAdditive()), coarsesolverconverged(true),
00305         coarseSmoother_(), verbosity_(parms.debugLevel())
00306     {
00307       assert(matrices_->isBuilt());
00308       
00309       // build the necessary smoother hierarchies
00310       matrices_->coarsenSmoother(smoothers_, smootherArgs_);
00311     }
00312 
00313     template<class M, class X, class S, class PI, class A>
00314     template<class C>
00315     AMG<M,X,S,PI,A>::AMG(const Operator& matrix,
00316                         const C& criterion,
00317                         const SmootherArgs& smootherArgs,
00318                         std::size_t gamma, std::size_t preSmoothingSteps,
00319                         std::size_t postSmoothingSteps,
00320                         bool additive_,
00321                         const PI& pinfo)
00322       : smootherArgs_(smootherArgs),
00323         smoothers_(), solver_(), scalarProduct_(0), gamma_(gamma),
00324         preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(true),
00325         additive(additive_), coarsesolverconverged(true),
00326         coarseSmoother_(), verbosity_(criterion.debugLevel())
00327     {
00328       dune_static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
00329                          "Matrix and Solver must match in terms of category!");
00330       // TODO: reestablish compile time checks.
00331       //dune_static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
00332       //                         "Matrix and Solver must match in terms of category!");
00333       Timer watch;
00334       matrices_ = new OperatorHierarchy(const_cast<Operator&>(matrix), pinfo);
00335             
00336       matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
00337       
00338       // build the necessary smoother hierarchies
00339       matrices_->coarsenSmoother(smoothers_, smootherArgs_);
00340 
00341       if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
00342         std::cout<<"Building Hierarchy of "<<matrices_->maxlevels()<<" levels took "<<watch.elapsed()<<" seconds."<<std::endl;
00343     }
00344 
00345     template<class M, class X, class S, class PI, class A>
00346     template<class C>
00347     AMG<M,X,S,PI,A>::AMG(const Operator& matrix,
00348                         const C& criterion,
00349                          const SmootherArgs& smootherArgs,
00350                         const PI& pinfo)
00351       : smootherArgs_(smootherArgs),
00352         smoothers_(), solver_(), scalarProduct_(0), 
00353         gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()), 
00354         postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
00355         additive(criterion.getAdditive()), coarsesolverconverged(true),
00356         coarseSmoother_(), verbosity_(criterion.debugLevel())
00357     {
00358       dune_static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
00359                          "Matrix and Solver must match in terms of category!");
00360       // TODO: reestablish compile time checks.
00361       //dune_static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
00362       //                         "Matrix and Solver must match in terms of category!");
00363       Timer watch;
00364       matrices_ = new OperatorHierarchy(const_cast<Operator&>(matrix), pinfo);
00365             
00366       matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
00367       
00368       // build the necessary smoother hierarchies
00369       matrices_->coarsenSmoother(smoothers_, smootherArgs_);
00370 
00371       if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
00372         std::cout<<"Building Hierarchy of "<<matrices_->maxlevels()<<" levels took "<<watch.elapsed()<<" seconds."<<std::endl;
00373     }
00374     
00375     template<class M, class X, class S, class PI, class A>
00376     AMG<M,X,S,PI,A>::~AMG()
00377     {
00378       if(buildHierarchy_){
00379         delete matrices_;
00380       }
00381       if(scalarProduct_)
00382         delete scalarProduct_;
00383     }
00384 
00385     
00386     template<class M, class X, class S, class PI, class A>
00387     void AMG<M,X,S,PI,A>::pre(Domain& x, Range& b)
00388     {
00389       if(smoothers_.levels()>0)
00390         smoothers_.finest()->pre(x,b);
00391       else
00392         // No smoother to make x consistent! Do it by hand
00393         matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
00394       Range* copy = new Range(b);
00395       rhs_ = new Hierarchy<Range,A>(*copy);
00396       Domain* dcopy = new Domain(x);
00397       lhs_ = new Hierarchy<Domain,A>(*dcopy);
00398       dcopy = new Domain(x);
00399       update_ = new Hierarchy<Domain,A>(*dcopy);
00400       matrices_->coarsenVector(*rhs_);
00401       matrices_->coarsenVector(*lhs_);
00402       matrices_->coarsenVector(*update_);
00403       
00404       // Preprocess all smoothers
00405       typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
00406       typedef typename Hierarchy<Range,A>::Iterator RIterator;
00407       typedef typename Hierarchy<Domain,A>::Iterator DIterator;
00408       Iterator coarsest = smoothers_.coarsest();
00409       Iterator smoother = smoothers_.finest();
00410       RIterator rhs = rhs_->finest();
00411       DIterator lhs = lhs_->finest();
00412       if(smoothers_.levels()>0){
00413           
00414       assert(lhs_->levels()==rhs_->levels());
00415       assert(smoothers_.levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
00416       assert(smoothers_.levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
00417       
00418       if(smoother!=coarsest)
00419         for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
00420           smoother->pre(*lhs,*rhs);
00421       smoother->pre(*lhs,*rhs);
00422       }
00423       
00424       
00425       // The preconditioner might change x and b. So we have to 
00426       // copy the changes to the original vectors.
00427       x = *lhs_->finest();
00428       b = *rhs_->finest();
00429       
00430       if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()){
00431         // We have the carsest level. Create the coarse Solver
00432         SmootherArgs sargs(smootherArgs_);
00433         sargs.iterations = 1;
00434         
00435         typename ConstructionTraits<Smoother>::Arguments cargs;
00436         cargs.setArgs(sargs);
00437         if(matrices_->redistributeInformation().back().isSetup()){
00438           // Solve on the redistributed partitioning     
00439           cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
00440           cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
00441           
00442           coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
00443           scalarProduct_ = ScalarProductChooser::construct(matrices_->parallelInformation().coarsest().getRedistributed());
00444         }else{    
00445           cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
00446           cargs.setComm(*matrices_->parallelInformation().coarsest());
00447           
00448           coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
00449           scalarProduct_ = ScalarProductChooser::construct(*matrices_->parallelInformation().coarsest());
00450         }
00451 #if HAVE_SUPERLU
00452       // Use superlu if we are purely sequential or with only one processor on the coarsest level.
00453         if(is_same<ParallelInformation,SequentialInformation>::value // sequential mode 
00454            || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
00455            || (matrices_->parallelInformation().coarsest().isRedistributed() 
00456                && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
00457                && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)){ // redistribute and 1 proc
00458           if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
00459           std::cout<<"Using superlu"<<std::endl;
00460           if(matrices_->parallelInformation().coarsest().isRedistributed())
00461             {
00462               if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
00463                 // We are still participating on this level
00464                 solver_  = new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat());
00465               else
00466                 solver_ = 0;
00467             }else
00468               solver_  = new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat());
00469         }else
00470 #endif
00471           {
00472             if(matrices_->parallelInformation().coarsest().isRedistributed())
00473               {
00474                 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
00475                   // We are still participating on this level
00476                   solver_ = new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()), 
00477                                                   *scalarProduct_, 
00478                                                   *coarseSmoother_, 1E-2, 10000, 0);
00479                 else
00480                   solver_ = 0;
00481               }else
00482               solver_ = new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()), 
00483                                               *scalarProduct_, 
00484                                               *coarseSmoother_, 1E-2, 1000, 0);
00485           }
00486       }
00487     }
00488     template<class M, class X, class S, class PI, class A>
00489     std::size_t AMG<M,X,S,PI,A>::levels()
00490     {
00491       return matrices_->levels();
00492     }
00493     template<class M, class X, class S, class PI, class A>
00494     std::size_t AMG<M,X,S,PI,A>::maxlevels()
00495     {
00496       return matrices_->maxlevels();
00497     }
00498 
00500     template<class M, class X, class S, class PI, class A>
00501     void AMG<M,X,S,PI,A>::apply(Domain& v, const Range& d)
00502     {
00503       if(additive){
00504         *(rhs_->finest())=d;
00505         additiveMgc();
00506         v=*lhs_->finest();
00507       }else{
00508         // Init all iterators for the current level
00509         initIteratorsWithFineLevel();
00510 
00511         
00512         *lhs = v;
00513         *rhs = d;
00514         *update=0;
00515         level=0;
00516                   
00517         mgc();
00518         
00519         if(postSteps_==0||matrices_->maxlevels()==1)
00520           pinfo->copyOwnerToAll(*update, *update);
00521         
00522         v=*update;
00523       }
00524       
00525     }
00526 
00527     template<class M, class X, class S, class PI, class A>
00528     void AMG<M,X,S,PI,A>::initIteratorsWithFineLevel()
00529     {
00530       smoother = smoothers_.finest();
00531       matrix = matrices_->matrices().finest();
00532       pinfo = matrices_->parallelInformation().finest();
00533       redist = 
00534         matrices_->redistributeInformation().begin();
00535       aggregates = matrices_->aggregatesMaps().begin();
00536       lhs = lhs_->finest();
00537       update = update_->finest();
00538       rhs = rhs_->finest();
00539     }
00540     
00541     template<class M, class X, class S, class PI, class A>
00542     bool AMG<M,X,S,PI,A>
00543     ::moveToCoarseLevel()
00544     {
00545       
00546       bool processNextLevel=true;
00547       
00548       if(redist->isSetup()){
00549         redist->redistribute(static_cast<const Range&>(*rhs), rhs.getRedistributed());
00550         processNextLevel =rhs.getRedistributed().size()>0;
00551         if(processNextLevel){           
00552           //restrict defect to coarse level right hand side.
00553           typename Hierarchy<Range,A>::Iterator fineRhs = rhs++;
00554           ++pinfo;
00555           Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
00556             ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(fineRhs.getRedistributed()), *pinfo);
00557         }
00558       }else{            
00559         //restrict defect to coarse level right hand side.
00560         typename Hierarchy<Range,A>::Iterator fineRhs = rhs++;
00561           ++pinfo;
00562           Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
00563 	    ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
00564       }
00565       
00566       if(processNextLevel){
00567         // prepare coarse system
00568         ++lhs;
00569         ++update;
00570         ++matrix;
00571         ++level;
00572         ++redist;
00573 
00574         if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()){
00575           // next level is not the globally coarsest one
00576           ++smoother;
00577           ++aggregates;
00578         }
00579         // prepare the update on the next level
00580         *update=0;
00581       }
00582       return processNextLevel;
00583     }
00584     
00585     template<class M, class X, class S, class PI, class A>
00586     void AMG<M,X,S,PI,A>
00587     ::moveToFineLevel(bool processNextLevel)
00588     {
00589       if(processNextLevel){
00590         if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()){
00591           // previous level is not the globally coarsest one
00592             --smoother;
00593             --aggregates;
00594         }
00595         --redist;
00596         --level;
00597         //prolongate and add the correction (update is in coarse left hand side)
00598         --matrix;
00599         
00600         //typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
00601         --lhs;  
00602         --pinfo;
00603       }
00604       if(redist->isSetup()){
00605         // Need to redistribute during prolongate
00606         lhs.getRedistributed()=0;
00607         Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
00608           ::prolongate(*(*aggregates), *update, *lhs, lhs.getRedistributed(), matrices_->getProlongationDampingFactor(), 
00609                        *pinfo, *redist);
00610       }else{
00611         *lhs=0;
00612         Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
00613           ::prolongate(*(*aggregates), *update, *lhs, 
00614                        matrices_->getProlongationDampingFactor(), *pinfo);
00615       }
00616       
00617       
00618       if(processNextLevel){
00619         --update;
00620         --rhs;
00621       }
00622       
00623       *update += *lhs;
00624     }
00625     
00626     
00627     template<class M, class X, class S, class PI, class A>
00628     void AMG<M,X,S,PI,A>
00629     ::presmooth()
00630     {
00631       
00632       for(std::size_t i=0; i < preSteps_; ++i){
00633             *lhs=0;
00634             SmootherApplier<S>::preSmooth(*smoother, *lhs, *rhs);
00635             // Accumulate update
00636             *update += *lhs;
00637             
00638             // update defect
00639             matrix->applyscaleadd(-1,static_cast<const Domain&>(*lhs), *rhs);
00640             pinfo->project(*rhs);
00641           }
00642     }
00643     
00644      template<class M, class X, class S, class PI, class A>
00645     void AMG<M,X,S,PI,A>
00646      ::postsmooth()
00647     { 
00648       
00649         for(std::size_t i=0; i < postSteps_; ++i){
00650           // update defect
00651           matrix->applyscaleadd(-1,static_cast<const Domain&>(*lhs), *rhs);
00652           *lhs=0;
00653           pinfo->project(*rhs);
00654           SmootherApplier<S>::postSmooth(*smoother, *lhs, *rhs);
00655           // Accumulate update
00656           *update += *lhs;
00657         }
00658     }
00659     
00660     
00661     template<class M, class X, class S, class PI, class A>
00662     bool AMG<M,X,S,PI,A>::usesDirectCoarseLevelSolver() const
00663     {
00664       return IsDirectSolver< CoarseSolver>::value;
00665     }
00666     
00667     template<class M, class X, class S, class PI, class A>
00668     void AMG<M,X,S,PI,A>::mgc(){
00669       if(matrix == matrices_->matrices().coarsest() && levels()==maxlevels()){
00670         // Solve directly
00671         InverseOperatorResult res;
00672         res.converged=true; // If we do not compute this flag will not get updated
00673         if(redist->isSetup()){
00674             redist->redistribute(*rhs, rhs.getRedistributed());
00675           if(rhs.getRedistributed().size()>0){
00676             // We are still participating in the computation
00677             pinfo.getRedistributed().copyOwnerToAll(rhs.getRedistributed(), rhs.getRedistributed());
00678             solver_->apply(update.getRedistributed(), rhs.getRedistributed(), res);
00679           }
00680           redist->redistributeBackward(*update, update.getRedistributed());
00681           pinfo->copyOwnerToAll(*update, *update);
00682         }else{
00683           pinfo->copyOwnerToAll(*rhs, *rhs);
00684           solver_->apply(*update, *rhs, res);
00685         }
00686 
00687         if (!res.converged)
00688           coarsesolverconverged = false;
00689       }else{
00690         // presmoothing
00691         presmooth();
00692 
00693 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
00694         bool processNextLevel = moveToCoarseLevel();
00695         
00696         if(processNextLevel){
00697           // next level
00698           for(std::size_t i=0; i<gamma_; i++)
00699             mgc();
00700         }
00701         
00702         moveToFineLevel(processNextLevel);
00703 #else
00704         *lhs=0;
00705 #endif  
00706 
00707         if(matrix == matrices_->matrices().finest()){
00708           coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
00709           if(!coarsesolverconverged)
00710             DUNE_THROW(MathError, "Coarse solver did not converge");
00711         }
00712         // postsmoothing
00713         postsmooth();
00714         
00715       }
00716     }
00717 
00718     template<class M, class X, class S, class PI, class A>
00719     void AMG<M,X,S,PI,A>::additiveMgc(){
00720       
00721       // restrict residual to all levels
00722       typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
00723       typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();      
00724       typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
00725       typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
00726       
00727       for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates){
00728         ++pinfo;
00729         Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
00730 	  ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
00731       }
00732       
00733       // pinfo is invalid, set to coarsest level
00734       //pinfo = matrices_->parallelInformation().coarsest
00735       // calculate correction for all levels
00736       lhs = lhs_->finest();
00737       typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_.finest();
00738       
00739       for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother){
00740         // presmoothing
00741         *lhs=0;
00742         smoother->apply(*lhs, *rhs);
00743       }
00744       
00745       // Coarse level solve
00746 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION 
00747       InverseOperatorResult res;
00748       pinfo->copyOwnerToAll(*rhs, *rhs);
00749       solver_->apply(*lhs, *rhs, res);
00750       
00751       if(!res.converged)
00752         DUNE_THROW(MathError, "Coarse solver did not converge");
00753 #else
00754       *lhs=0;
00755 #endif
00756       // Prologate and add up corrections from all levels
00757       --pinfo;
00758       --aggregates;
00759       
00760       for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo){
00761         Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
00762 	  ::prolongate(*(*aggregates), *coarseLhs, *lhs, 1, *pinfo);
00763       }
00764     }
00765 
00766     
00768     template<class M, class X, class S, class PI, class A>
00769     void AMG<M,X,S,PI,A>::post(Domain& x)
00770     {
00771       if(buildHierarchy_){
00772         if(solver_)
00773           delete solver_;
00774         if(coarseSmoother_)
00775           ConstructionTraits<Smoother>::deconstruct(coarseSmoother_);
00776       }
00777       
00778       // Postprocess all smoothers
00779       typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
00780       typedef typename Hierarchy<Range,A>::Iterator RIterator;
00781       typedef typename Hierarchy<Domain,A>::Iterator DIterator;
00782       Iterator coarsest = smoothers_.coarsest();
00783       Iterator smoother = smoothers_.finest();
00784       DIterator lhs = lhs_->finest();
00785       if(smoothers_.levels()>0){
00786         if(smoother != coarsest  || matrices_->levels()<matrices_->maxlevels())
00787           smoother->post(*lhs);
00788         if(smoother!=coarsest)
00789           for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
00790             smoother->post(*lhs);
00791         smoother->post(*lhs);
00792       }
00793 
00794       delete &(*lhs_->finest());
00795       delete lhs_;
00796       delete &(*update_->finest());
00797       delete update_;
00798       delete &(*rhs_->finest());
00799       delete rhs_;
00800     }
00801 
00802     template<class M, class X, class S, class PI, class A>
00803     template<class A1>
00804     void AMG<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
00805     {
00806       matrices_->getCoarsestAggregatesOnFinest(cont);
00807     }
00808     
00809   } // end namespace Amg
00810 }// end namespace Dune
00811   
00812 #endif