multiplicativeMultigrid.hh

Go to the documentation of this file.

 
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2015-2020 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef MG_MULTIPLICATIVEMULTIGRID_HH
#define MG_MULTIPLICATIVEMULTIGRID_HH
 
#include <type_traits>
 
#include <boost/timer/timer.hpp>
 
#include <dune/istl/preconditioner.hh>
 
#include "fem/spaces.hh"
#include "linalg/conjugation.hh"
#include "linalg/direct.hh"
#include "linalg/domainDecompositionPreconditioner.hh"
#include "linalg/jacobiPreconditioner.hh"
#include "linalg/symmetricOperators.hh"
#include "mg/prolongation.hh"
#include "utilities/memory.hh"
#include "utilities/timing.hh"
 
namespace Kaskade
{
 
 
  // ----------------------------------------------------------------------------------------------------------------------
 
  template <class Entry, class Index, class Smoother, class Prolongation>
  class MultiplicativeMultiGrid: public SymmetricPreconditioner<typename Smoother::domain_type,typename Smoother::range_type>
  {
  public:
    using field_type = typename Smoother::field_type;
    using domain_type = typename Smoother::domain_type;
    using range_type = typename Smoother::range_type;
    using CoarsePreconditioner = Dune::Preconditioner<domain_type,range_type>; // symmetric would be better? but elaborate implementation for direct solvers
 
    MultiplicativeMultiGrid() = default;
 
    template <class MakeSmoother>
    MultiplicativeMultiGrid(MultiGridStack<Prolongation,Entry,Index>&& mgStack_, MakeSmoother const& makeSmoother,
                            std::unique_ptr<CoarsePreconditioner>&& coarsePreconditioner_, int nPre_=3, int nPost_=3)
    : mgStack(std::move(mgStack_)), nPre(nPre_), nPost(nPost_), coarsePreconditioner(std::move(coarsePreconditioner_)),
      linesearch(false), smootherStepSize(1.0)
    {
      // Create the smoothers for all levels.
      Timings& timer = Timings::instance();
      timer.start("smoother construction");
      for (int l=1; l<mgStack.levels(); ++l)
        smoothers.push_back(makeSmoother(mgStack.a(l)));
      timer.stop("smoother construction");
    }
 
    MultiplicativeMultiGrid(MultiplicativeMultiGrid&& other) = default;
 
    MultiplicativeMultiGrid& operator=(MultiplicativeMultiGrid&& other) = default;
 
 
    virtual void apply(domain_type& x, range_type const& r)
    {
      assert(x.two_norm()==0);
      range_type b = r;
      runMG(x,b,mgStack.levels()-1);
    }
 
    virtual field_type applyDp(domain_type& x, range_type const& r)
    {
      apply(x,r);
      return x*r;
    }
 
    virtual bool requiresInitializedInput() const
    {
      return true;
    }
 
    void setSmoothings(int nPre_, int nPost_)
    {
      nPre = nPre_;
      nPost = nPost_;
      assert(nPre>=0 && nPost>=0 && nPre+nPost>0);
    }
 
    void setSmootherStepSize(double w)
    {
      assert(w>0);
      smootherStepSize = w;
    }
 
    double getSmootherStepSize() const { return smootherStepSize; }
 
    void setLinesearch(bool ls)
    {
      linesearch = ls;
    }
 
    CoarsePreconditioner& getCoarsePreconditioner()
    {
      return *coarsePreconditioner;
    }
 
    std::pair<int,int> getSmoothings()
    {
      return std::make_pair(nPre,nPost);
    }
 
  private:
    // precondition: x==0
    void runMG(domain_type& x, range_type& r, int level)
    {
      Timings& timer = Timings::instance();
 
      if (level==0)
      {
        timer.start("coarse grid solver");
        coarsePreconditioner->pre(x,r);
        coarsePreconditioner->apply(x,r);
        coarsePreconditioner->post(x);
        timer.stop("coarse grid solver");
      }
      else
      {
        auto& s = smoothers[level-1];             // Smoother S
        auto const& a = mgStack.a(level);         // Galerkin matrix A
        double w = smootherStepSize;
        domain_type dx(x);                        // temporary vector
        range_type adx(r);
 
        // pre-smoothing
        bool const sRequiresZeroInput = s.requiresInitializedInput();
        s.pre(x,r);
        for (int i=0; i<nPre; ++i)
        {
          if (sRequiresZeroInput) dx = 0;
          timer.start("smoother application");
          s.apply(dx,r);                          // dx = B^{-1} r
          timer.stop("smoother application");
 
          timer.start("matrix-vector");
          if (linesearch)
          {
            double dxadx = a.mv(dx,adx);          // adx = A*dx           residual update direction
            w = dx*r / (dxadx);                   // w = dx*r / dx*A*dx   "optimal" step length
            r.axpy(-w,adx);                       // r = r - w*adx        update residual
          }
          else
            a.usmv(-w,dx,r);                      // r = r - w*A*dx       update residual (without intermediate adx access)
          timer.stop("matrix-vector");
 
          x.axpy(w,dx);                           // x = x + w*dx         update iterate
        }
        s.post(x);
 
        // coarse grid correction
        auto const& p = mgStack.p(level-1);
        timer.start("restriction");
        range_type cr; cr.resize(p.M());
        p.mtv(r,cr);                             // cr = P^T r
        timer.stop("restriction");
        domain_type cx; cx.resize(p.M()); cx = 0;
        runMG(cx,cr,level-1);                    // cx ~ (P^T A P)^{-1} cr
        timer.start("prolongation");
        p.mv(cx,dx);                             // dx = P*cx
        timer.stop("prolongation");
        x += dx;                                 // x = x+dx
        a.usmv(-1,dx,r);                         // r = r-A*dx
 
 
        // post-smoothing
        s.pre(x,r);
        for (int i=0; i<nPost; ++i)
        {
          if (sRequiresZeroInput) dx = 0;
          timer.start("smoother application");
          s.apply(dx,r);                          // dx = B^{-1} r
          timer.stop("smoother application");
 
          if (i+1<nPost || linesearch)            // update residual only if needed lateron
          {
            timer.start("matrix-vector");
            if (linesearch)
            {
              double dxadx = a.mv(dx,adx);        // adx = A*dx
              w = dx*r / (dxadx);                 // w = dx*r / dx*A*dx
              r.axpy(-w,adx);                     // r = r - w*adx
            }
            else
              a.usmv(-w,dx,r);                    // r = r - w*A*dx
            timer.stop("matrix-vector");
          }
          x.axpy(w,dx);                           // x = x + w*dx
        }
        s.post(x);
 
        // optional line search
        if (linesearch)
        {
          a.mv(x,dx);                             // dx = Ax
          double omega = (x*r) / (dx*x);          // omega = r^T x / x^T A x
          x *= 1+omega;                           // x = x + omega * x
        }
      }
    }
 
    MultiGridStack<Prolongation,Entry,Index>  mgStack;
    std::vector<Smoother>                     smoothers;
    int                                       nPre, nPost;
    std::unique_ptr<CoarsePreconditioner>     coarsePreconditioner;
    bool                                      linesearch;
    double                                    smootherStepSize;
  };
 
  // ----------------------------------------------------------------------------------------------------------------------
 
  template <class Entry, class Index, class Prolongation, class MakeSmoother, class CoarsePreconditioner>
  auto makeMultiplicativeMultiGrid(MultiGridStack<Prolongation,Entry,Index>&& mgStack,
                                   MakeSmoother const& makeSmoother,
                                   std::unique_ptr<CoarsePreconditioner>&& coarsePreconditioner,
                                   int nPre=3, int nPost=3)
  {
    using Smoother = std::result_of_t<MakeSmoother(NumaBCRSMatrix<Entry,Index>)>;
    return MultiplicativeMultiGrid<Entry,Index,Smoother,Prolongation>(
                std::move(mgStack),makeSmoother,std::move(coarsePreconditioner),nPre,nPost);
  }
 
  // ----------------------------------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------------------------------
 
  class MakeJacobiSmoother
  {
  public:
    template <typename Matrix>
    auto operator()(Matrix const& a) const
    {
      return JacobiPreconditioner(a);
    }
  };
 
  template <typename Space, typename Scalar=float>
  class MakeAdditiveSchwarzSmoother
  {
  public:
    MakeAdditiveSchwarzSmoother(Space const& space_, double targetCondition_=-1.0)
    : space(space_), targetCondition(targetCondition_)
    {}
 
    template <typename Entry, typename Index>
    auto operator()(NumaBCRSMatrix<Entry,Index> const& a) const
    {
      return PatchDomainDecompositionPreconditioner<Space,Entry::rows,NestedDissectionStorageTag<Scalar,Scalar>,
                                                    Index>(space,a,targetCondition);
    }
 
  private:
    Space const& space;
    double targetCondition;
  };
 
  template <typename Matrix>
  auto makeDirectPreconditioner(Matrix&& A, DirectType directType=DirectType::MUMPS)
  {
    return DirectSolver<typename MatrixTraits<Matrix>::NaturalDomain,typename MatrixTraits<Matrix>::NaturalRange>(A,directType);
  }
 
  template <class Domain, class Range, class SparseIndex>
  class DirectPreconditionerFloatWrapper : public Dune::Preconditioner<Domain,Range>
  {
  public:
    using FloatDomain = Domain;
    using FloatRange = Range;
    static constexpr int domainBlockSize = Domain::block_type::dimension;
    static constexpr int rangeBlockSize = Range::block_type::dimension;
    using FloatMatrix = NumaBCRSMatrix<Dune::FieldMatrix<float,rangeBlockSize,domainBlockSize>,SparseIndex>;
    using DoubleMatrix = NumaBCRSMatrix<Dune::FieldMatrix<double,rangeBlockSize,domainBlockSize>,SparseIndex>;
    using DoubleDomain = typename MatrixTraits<DoubleMatrix>::NaturalDomain;
    using DoubleRange = typename MatrixTraits<DoubleMatrix>::NaturalRange;
 
    DirectPreconditionerFloatWrapper(FloatMatrix const& matA, DirectType directType) : directSolver(DoubleMatrix(matA), directType) {}
 
    virtual void pre (FloatDomain &, FloatRange &) {}
 
    virtual void apply (FloatDomain & x, FloatRange const& b)
    {
      DoubleDomain xd(x.size());
      DoubleRange bd(b.size());
      auto endBd = vectorToSequence(b, &(bd[0][0]));
      if(endBd != &(bd[0][0]) + bd.dim()) {
        throw DirectSolverException("Copying float to double block-vector failed. End Iterators differ.",__FILE__,__LINE__);
      }
      directSolver.apply(xd, bd);
      auto endX = vectorToSequence(xd, &(x[0][0]));
      if(endX != &(x[0][0]) + x.dim()) {
        throw DirectSolverException("Copying float to double block-vector failed. End Iterators differ.",__FILE__,__LINE__);
      }
    }
 
    virtual void post (FloatDomain &) {}
 
  private:
    DirectSolver<DoubleDomain, DoubleRange> directSolver;
  };
 
  template <int domainBlockSize, int rangeBlockSize, class Index>
  auto makeDirectPreconditioner(NumaBCRSMatrix<Dune::FieldMatrix<float,rangeBlockSize,domainBlockSize>,Index> && A, DirectType directType=DirectType::MUMPS)
  {
    using Matrix = NumaBCRSMatrix<Dune::FieldMatrix<float,rangeBlockSize,domainBlockSize>,Index>;
    return DirectPreconditionerFloatWrapper<typename MatrixTraits<Matrix>::NaturalDomain,typename MatrixTraits<Matrix>::NaturalRange,Index>(A, directType);
  }
 
  // ----------------------------------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------------------------------
 
  template <class Entry, class Index, class GridMan>
  auto makeJacobiMultiGrid(NumaBCRSMatrix<Entry,Index> const& A, GridMan const& gridman,
                           int nPre=3, int nPost=3, bool onlyLowerTriangle=false, DirectType directType=DirectType::MUMPS)
  {
    Timings& timer = Timings::instance();
 
    timer.start("MG stack creation");
    auto mgStack = makeGeometricMultiGridStack(gridman,duplicate(A),onlyLowerTriangle);
    timer.stop("MG stack creation");
 
    timer.start("direct solver creation");
    auto coarsePreconditioner = makeDirectPreconditioner(std::move(mgStack.coarseGridMatrix()),directType);
    timer.stop("direct solver creation");
 
    timer.start("MG creation");
    auto mg = makeMultiplicativeMultiGrid(std::move(mgStack),MakeJacobiSmoother(),moveUnique(std::move(coarsePreconditioner)),nPre,nPost);
    timer.stop("MG creation");
    mg.setSmootherStepSize(0.5);
 
    return mg;
  }
 
  template <class Entry, class Index, class Space, class MakeSmoother>
  auto makeMultiplicativePMultiGrid(NumaBCRSMatrix<Entry,Index>&& A, Space const& space,
                                    MakeSmoother const& makeSmoother, int nPre=3, int nPost=3, bool onlyLowerTriangle=false, DirectType directType=DirectType::MUMPS)
  {
    Timings& timer = Timings::instance();
    timer.start("create P multigrid stack");
    auto mgStack = makePMultiGridStack(space,std::move(A),onlyLowerTriangle);
    timer.stop("create P multigrid stack");
 
    timer.start("create MG as coarse solver");
    auto coarsePreconditioner = makeJacobiMultiGrid(std::move(mgStack.coarseGridMatrix()),space.gridManager(),nPre,nPost,onlyLowerTriangle,directType);
    coarsePreconditioner.setSmootherStepSize(0.6);
    timer.stop("create MG as coarse solver");
 
    timer.start("create multigrid");
    auto mg = makeMultiplicativeMultiGrid(std::move(mgStack),makeSmoother,moveUnique(std::move(coarsePreconditioner)),nPre,nPost);
    timer.stop("create multigrid");
 
    return mg;
  }
 
  template <class Entry, class Index, class Space, class CoarseSpace>
  auto makeBlockJacobiPMultiGrid(NumaBCRSMatrix<Entry,Index> A, Space const& space,
                                 NumaBCRSMatrix<Entry,Index> coarseA, CoarseSpace const& coarseSpace,
                                 int nPre=3, int nPost=3, DirectType directType=DirectType::MUMPS)
  {
    Timings& timer = Timings::instance();
    timer.start("create P multigrid stack");
    auto mgStack = makePMultiGridStack(space,std::move(A),coarseSpace,std::move(coarseA));
    timer.stop("create P multigrid stack");
 
    timer.start("create MG as coarse solver");
    auto coarsePreconditioner = makeJacobiMultiGrid(std::move(mgStack.coarseGridMatrix()),space.gridManager(),nPre,nPost,false,directType);
    coarsePreconditioner.setSmootherStepSize(0.6);
    timer.stop("create MG as coarse solver");
 
    timer.start("create multigrid");
    auto mg = makeMultiplicativeMultiGrid(std::move(mgStack),MakeAdditiveSchwarzSmoother<Space>(space),
                                          moveUnique(std::move(coarsePreconditioner)),nPre,nPost);
    mg.setSmootherStepSize(Space::Grid::dimension==3 ? 0.4 : 0.5);
    timer.stop("create multigrid");
 
    return mg;
  }
 
  template <class Entry, class Index, class Space>
  auto makeBlockJacobiPMultiGrid(NumaBCRSMatrix<Entry,Index> A, Space const& space,
                                 int nPre=3, int nPost=3, bool onlyLowerTriangle=false, DirectType directType=DirectType::MUMPS)
  {
    auto mg = makeMultiplicativePMultiGrid(std::move(A),space,MakeAdditiveSchwarzSmoother<Space>(space),
                                           nPre,nPost,onlyLowerTriangle,directType);
    mg.setSmootherStepSize(Space::Grid::dimension==3 ? 0.4 : 0.5);
    return mg;
  }
 
  template <class Entry, class Index, class FineSpace>
  auto makeJacobiPMultiGrid(NumaBCRSMatrix<Entry,Index> A, FineSpace const& space,
                            int nPre=3, int nPost=3, bool onlyLowerTriangle=false, DirectType directType=DirectType::MUMPS)
  {
    auto mg = makeMultiplicativePMultiGrid(std::move(A),space,MakeJacobiSmoother(),nPre,nPost,onlyLowerTriangle,directType);
    mg.setSmootherStepSize(0.5);
    return mg;
  }
 
 
 
 
  // ---------------------------------------------------------------------------------------------------------
 
 
}
 
#endif