kamg.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_AMG_KAMG_HH
#define DUNE_AMG_KAMG_HH
 
#include <dune/istl/preconditioners.hh>
#include "amg.hh"
 
namespace Dune
{
  namespace Amg
  {
 
    template<class AMG>
    class KAmgTwoGrid
      : public Preconditioner<typename AMG::Domain,typename AMG::Range>
    {
      typedef typename AMG::Domain Domain;
      typedef typename AMG::Range Range;
    public:
 
      enum {
        category = AMG::category
      };
 
      KAmgTwoGrid(AMG& amg, shared_ptr<InverseOperator<Domain,Range> > coarseSolver)
        : amg_(amg), coarseSolver_(coarseSolver)
      {}
 
      void pre(typename AMG::Domain& x, typename AMG::Range& b)
      {}
 
      void post(typename AMG::Domain& x)
      {}
 
      void apply(typename AMG::Domain& v, const typename AMG::Range& d)
      {
        // Copy data
        *levelContext_->update=0;
        *levelContext_->rhs = d;
        *levelContext_->lhs = v;
 
        presmooth(*levelContext_, amg_.preSteps_);
        bool processFineLevel =
          amg_.moveToCoarseLevel(*levelContext_);
 
        if(processFineLevel) {
          typename AMG::Range b=*levelContext_->rhs;
          typename AMG::Domain x=*levelContext_->update;
          InverseOperatorResult res;
          coarseSolver_->apply(x, b, res);
          *levelContext_->update=x;
        }
 
        amg_.moveToFineLevel(*levelContext_, processFineLevel);
 
        postsmooth(*levelContext_, amg_.postSteps_);
        v=*levelContext_->update;
      }
 
      InverseOperator<Domain,Range>* coarseSolver()
      {
        return coarseSolver_;
      }
 
      void setLevelContext(Dune::shared_ptr<typename AMG::LevelContext> p)
      {
        levelContext_=p;
      }
 
      ~KAmgTwoGrid()
      {}
 
    private:
      AMG& amg_;
      Dune::shared_ptr<InverseOperator<Domain,Range> > coarseSolver_;
      Dune::shared_ptr<typename AMG::LevelContext> levelContext_;
    };
 
 
 
    template<class M, class X, class S, class PI=SequentialInformation,
        class K=GeneralizedPCGSolver<X>, class A=std::allocator<X> >
    class KAMG : public Preconditioner<X,X>
    {
    public:
      typedef AMG<M,X,S,PI,A> Amg;
      typedef K KrylovSolver;
      typedef typename Amg::OperatorHierarchy OperatorHierarchy;
      typedef typename Amg::CoarseSolver CoarseSolver;
      typedef typename Amg::ParallelInformation ParallelInformation;
      typedef typename Amg::SmootherArgs SmootherArgs;
      typedef typename Amg::Operator Operator;
      typedef typename Amg::Domain Domain;
      typedef typename Amg::Range Range;
      typedef typename Amg::ParallelInformationHierarchy ParallelInformationHierarchy;
      typedef typename Amg::ScalarProduct ScalarProduct;
 
      enum {
        category = Amg::category
      };
      KAMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
           const SmootherArgs& smootherArgs, std::size_t gamma,
           std::size_t preSmoothingSteps =1, std::size_t postSmoothingSteps = 1,
           std::size_t maxLevelKrylovSteps = 3 , double minDefectReduction =1e-1);
 
      template<class C>
      KAMG(const Operator& fineOperator, const C& criterion,
           const SmootherArgs& smootherArgs, std::size_t gamma=1,
           std::size_t preSmoothingSteps=1, std::size_t postSmoothingSteps=1,
           std::size_t maxLevelKrylovSteps=3, double minDefectReduction=1e-1,
           const ParallelInformation& pinfo=ParallelInformation());
 
      void pre(Domain& x, Range& b);
      void post(Domain& x);
      void apply(Domain& v, const Range& d);
 
      std::size_t maxlevels();
 
    private:
      Amg amg;
 
      std::size_t maxLevelKrylovSteps;
 
      double levelDefectReduction;
 
      std::vector<shared_ptr<typename Amg::ScalarProduct> > scalarproducts;
 
      std::vector<shared_ptr<KAmgTwoGrid<Amg> > > ksolvers;
    };
 
    template<class M, class X, class S, class P, class K, class A>
    KAMG<M,X,S,P,K,A>::KAMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
                            const SmootherArgs& smootherArgs,
                            std::size_t gamma, std::size_t preSmoothingSteps,
                            std::size_t postSmoothingSteps,
                            std::size_t ksteps, double reduction)
      : amg(matrices, coarseSolver, smootherArgs, gamma, preSmoothingSteps,
            postSmoothingSteps), maxLevelKrylovSteps(ksteps), levelDefectReduction(reduction)
    {}
 
    template<class M, class X, class S, class P, class K, class A>
    template<class C>
    KAMG<M,X,S,P,K,A>::KAMG(const Operator& fineOperator, const C& criterion,
                            const SmootherArgs& smootherArgs, std::size_t gamma,
                            std::size_t preSmoothingSteps, std::size_t postSmoothingSteps,
                            std::size_t ksteps, double reduction,
                            const ParallelInformation& pinfo)
      : amg(fineOperator, criterion, smootherArgs, gamma, preSmoothingSteps,
            postSmoothingSteps, false, pinfo), maxLevelKrylovSteps(ksteps), levelDefectReduction(reduction)
    {}
 
 
    template<class M, class X, class S, class P, class K, class A>
    void KAMG<M,X,S,P,K,A>::pre(Domain& x, Range& b)
    {
      amg.pre(x,b);
      scalarproducts.reserve(amg.levels());
      ksolvers.reserve(amg.levels());
 
      typename OperatorHierarchy::ParallelMatrixHierarchy::Iterator
      matrix = amg.matrices_->matrices().coarsest();
      typename ParallelInformationHierarchy::Iterator
      pinfo = amg.matrices_->parallelInformation().coarsest();
      bool hasCoarsest=(amg.levels()==amg.maxlevels());
 
      if(hasCoarsest) {
        if(matrix==amg.matrices_->matrices().finest())
          return;
        --matrix;
        --pinfo;
        ksolvers.push_back(shared_ptr<KAmgTwoGrid<Amg> >(new KAmgTwoGrid<Amg>(amg, amg.solver_)));
      }else
        ksolvers.push_back(shared_ptr<KAmgTwoGrid<Amg> >(new KAmgTwoGrid<Amg>(amg, shared_ptr<InverseOperator<Domain,Range> >())));
 
      std::ostringstream s;
 
      if(matrix!=amg.matrices_->matrices().finest())
        while(true) {
          scalarproducts.push_back(shared_ptr<typename Amg::ScalarProduct>(Amg::ScalarProductChooser::construct(*pinfo)));
          shared_ptr<InverseOperator<Domain,Range> > ks =
            shared_ptr<InverseOperator<Domain,Range> >(new KrylovSolver(*matrix, *(scalarproducts.back()),
                                                                        *(ksolvers.back()), levelDefectReduction,
                                                                        maxLevelKrylovSteps, 0));
          ksolvers.push_back(shared_ptr<KAmgTwoGrid<Amg> >(new KAmgTwoGrid<Amg>(amg, ks)));
          --matrix;
          --pinfo;
          if(matrix==amg.matrices_->matrices().finest())
            break;
        }
    }
 
 
    template<class M, class X, class S, class P, class K, class A>
    void KAMG<M,X,S,P,K,A>::post(Domain& x)
    {
      amg.post(x);
 
    }
 
    template<class M, class X, class S, class P, class K, class A>
    void KAMG<M,X,S,P,K,A>::apply(Domain& v, const Range& d)
    {
      if(ksolvers.size()==0)
      {
        Range td=d;
        InverseOperatorResult res;
        amg.solver_->apply(v,td,res);
      }else
      {
        typedef typename Amg::LevelContext LevelContext;
        Dune::shared_ptr<LevelContext> levelContext(new LevelContext);
        amg.initIteratorsWithFineLevel(*levelContext);
        typedef typename std::vector<shared_ptr<KAmgTwoGrid<Amg> > >::iterator Iter;
        for(Iter solver=ksolvers.begin(); solver!=ksolvers.end(); ++solver)
          (*solver)->setLevelContext(levelContext);
        ksolvers.back()->apply(v,d);
      }
    }
 
    template<class M, class X, class S, class P, class K, class A>
    std::size_t KAMG<M,X,S,P,K,A>::maxlevels()
    {
      return amg.maxlevels();
    }
 
  } // Amg
} // Dune
 
#endif