linalg/apcg.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) 2012-2018 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef NMIII_APCG_HH
#define NMIII_APCG_HH
 
#include <numeric>
#include <string>
 
#include <boost/circular_buffer.hpp>
#include <boost/timer/timer.hpp>
 
#include "dune/common/timer.hh"
#include "dune/istl/istlexception.hh"
#include "dune/istl/operators.hh"
#include "dune/istl/preconditioners.hh"
#include "dune/istl/scalarproducts.hh"
#include "dune/istl/solvers.hh"
 
#include "linalg/symmetricOperators.hh"
#include "utilities/detailed_exception.hh"
#include "utilities/geometric_sequence.hh"
#include "utilities/scalar.hh"
 
namespace Kaskade
{
  template <class Domain, class Range>
  bool requiresInitializedInput(Dune::Preconditioner<Domain,Range> const* p)
  {
    if (dynamic_cast<Dune::Richardson<Domain,Range> const*>(p))
      return false;
    if (auto sp = dynamic_cast<SymmetricPreconditioner<Domain,Range> const*>(p))
      return sp->requiresInitializedInput();
    return true;
  }
 
  // ----------------------------------------------------------------------------------------------
  
  template <class R>
  class PCGTerminationCriterion 
  {
  public:
    typedef R Real;
    
    virtual void clear() = 0;
    
    virtual PCGTerminationCriterion<R>& tolerance(Real tol) = 0;
    
    virtual void step(Real gamma2) = 0;
    
    virtual void residual(Real sigma) = 0;
    
    virtual operator bool() const = 0;
  };
  
  template <class R>
  class PCGCountingTerminationCriterion: public PCGTerminationCriterion<R>
  {
  public:
    typedef R Real;
    
    PCGCountingTerminationCriterion(int n_)
    : n(n_), count(0)
    {}
    
    virtual void clear()
    {
      count = 0;
    }
    
    virtual PCGTerminationCriterion<R>& tolerance(Real ) 
    {
      return *this;
    }
    
    virtual void step(Real ) 
    {
      ++count;
    }
    
    virtual void residual(Real ) 
    {}
    
    virtual operator bool() const 
    {
      return count >= n;
    }
    
  private: 
    int n;
    int count;
  };
 
  template <class R>
  class PCGEnergyErrorTerminationCriterion: public PCGTerminationCriterion<R>
  {
  public:
    typedef R Real;
    
    PCGEnergyErrorTerminationCriterion(Real atol, int maxit);
    
    virtual void clear();
    
    virtual PCGEnergyErrorTerminationCriterion<Real>& tolerance(Real atol);
    
     PCGEnergyErrorTerminationCriterion<Real>& lookahead(int lah);
    
    virtual void step(Real gamma2);
    
    virtual void residual(Real sigma) {}
    
    virtual operator bool() const;
    
    Real error() const;
    
    void setMaxit(int m)
    {
      maxit = m;
    }
    
    std::vector<Real> const& gamma2() const {return gammas2; }
    
  private:
    Real tol;
    int maxit;
    // squared gammas
    std::vector<Real> gammas2;
    int lookah;
  };
  
 
  template<class X, class Xstar>
  class Pcg : public Dune::InverseOperator<X,Xstar> 
  {
  public:
    typedef X domain_type;
    typedef Xstar range_type;
    typedef typename X::field_type field_type;
    
    typedef typename ScalarTraits<field_type>::Real Real;
 
    Pcg(SymmetricLinearOperator<X,Xstar>& op_, SymmetricPreconditioner<X,Xstar>& prec_, 
        PCGTerminationCriterion<Real>& terminate_, int verbose_=0)
    : proxyOp(op_,zeroDp), op(op_),
      proxyPrec(prec_,zeroDp), prec(prec_), 
      terminate(terminate_), verbose(verbose_), requiresInit(requiresInitializedInput(&prec))
    { }
 
    Pcg(Dune::LinearOperator<X,Xstar>& op_, Dune::Preconditioner<X,Xstar>& prec_,
        DualPairing<X,Xstar> const& dp_, PCGTerminationCriterion<Real>& terminate_, int verbose_=0)
    : proxyOp(op_,dp_), op(proxyOp),
      proxyPrec(prec_,dp_), prec(proxyPrec), 
      terminate(terminate_), verbose(verbose_), requiresInit(requiresInitializedInput(&prec))
    {    }
    
    virtual void apply (X& u, Xstar& b, Dune::InverseOperatorResult& res) 
    {
      res.clear();                // clear solver statistics
      terminate.clear();          // clear termination criterion
      Dune::Timer watch;          // start a timer
      
      boost::timer::cpu_timer mvTimer;
      mvTimer.stop();
      boost::timer::cpu_timer apTimer;
      apTimer.stop();
      boost::timer::cpu_timer prTimer;
      prTimer.stop();
    
      Xstar r(b);
      op.applyscaleadd(-1,u,r);  // r = b-Au
 
      // Keeping rq and q in unique_ptr allows efficient computation of q = rq + beta*q below
      // by swapping the pointers.
      // More elegant would be a swap method directly on X, but this is hard to implement
      // as there is no swap either for boost::fusion vectors nor for Dune::BlockVectors.
      // Even more elegant would be a loop fusion expression framework for BLAS level 1 operations.
      std::unique_ptr<X> rq(new X(u)); *rq = 0;
      prec.apply(*rq,r); // rq = B^{-1} r
      
      std::unique_ptr<X> q(new X(*rq));
 
      Xstar Aq(b);
    
 
      // some local variables
      field_type alpha,beta,sigma,gamma;
      sigma = op.dp(*rq,r); // preconditioned residual norm squard
      terminate.residual(ScalarTraits<field_type>::real(sigma));
      field_type const sigma0 = std::abs(sigma);
      
      // Check for trivial case
      if (sigma0 == 0)
        return;
 
      if (verbose>0) 
      {            // printing
        std::cout << "=== Kaskade7 PCG ===\n";
        if (verbose>1)
          std::cout << "Iter   precond res.       ratio   avg. rate  step ener  \n"
                    << "     0" << std::setw(13) << sigma0 << std::endl;
      }
 
      // the PCG loop
      int i=1; 
      for ( ; true; i++ ) 
      {
        // minimize in given search direction p
        mvTimer.resume();
        field_type qAq = op.applyDp(*q,Aq);             // compute Aq = A*q and qAq = <q,A*q>
        if (std::isnan(qAq)) std:: cerr << "qAq is nan\n";        
        mvTimer.stop();
        
        if (std::abs(qAq) < 1e-28*sigma) 
          break; // oops.
          
        if (qAq <= 0)
          throw NonpositiveMatrixException("encountered nonpositive energy product " + std::to_string(qAq)
                                           + " vs sigma " + std::to_string(sigma) + " in PCG.",
                                           __FILE__,__LINE__);
          
        alpha = sigma/qAq;
        apTimer.resume();
        u.axpy(alpha,*q);
        apTimer.stop();
        gamma = sigma*alpha;
        terminate.step(ScalarTraits<field_type>::real(gamma));
        resDebug.push_back(std::sqrt(sigma));
 
        // convergence test
        if (terminate)
          break;
 
        apTimer.resume();
        r.axpy(-alpha,Aq); // r = r - alpha*A*q
        apTimer.stop();
        
        prTimer.resume();
        if (requiresInit)
          *rq = 0;
        field_type sigmaNew = prec.applyDp(*rq,r); // compute rq = B^{-1}r and <rq,r>
        prTimer.stop();
 
        // determine new search direction
        terminate.residual(ScalarTraits<field_type>::real(sigmaNew));
        if (verbose>1)             // print
        {
          std::cout << std::setw(6) << i 
                    << std::setw(13) << std::scientific << std::setprecision(5) << std::abs(sigmaNew) 
                    << std::setw(12) << std::fixed      << std::setprecision(3) << std::abs(sigmaNew/sigma) 
                    << std::setw(12) << std::fixed      << std::setprecision(3) << std::pow(std::abs(sigmaNew)/sigma0,1.0/i) 
                    << std::setw(12) << std::scientific << std::setprecision(5) << gamma;
          if (verbose>2)
          {
            Xstar Au = u; Au = 0;
            op.apply(u,Au);
            auto e = op.dp(u,Au);
            Au *= 0.5;
            Au -= b;
            std::cout << " " << op.dp(u,Au) << " " << e << " " << u.two_norm();
          }
          std::cout << std::endl;
        }
 
        // convergence check to prevent division by zero if we obtained an exact solution
        // (which may happen for low dimensional systems)
        if (std::abs(sigmaNew) < 1e-28*sigma0)
          break;
        
        beta = sigmaNew/sigma;
        sigma = sigmaNew;
        apTimer.resume();
        rq->axpy(beta,*q); // compute q = rq + beta*q
        std::swap(rq,q);
        apTimer.stop();
      }
 
      if (verbose>0)                // printing for non verbose
        std::cout << std::setw(6) << i << std::setw(12) << std::abs(sigma) << std::endl;
 
      res.iterations = i;               // fill statistics
      res.reduction = std::sqrt(std::abs(sigma)/sigma0);
      res.conv_rate  = pow(res.reduction,1.0/i);
      res.elapsed = watch.elapsed();
 
      if (verbose>0)                 // final print 
      {
        std::cout << "=== rate=" << res.conv_rate
                  << ", T=" << res.elapsed
                  << ", TIT=" << res.elapsed/i
                  << ", IT=" << i << std::endl;
        
        std::cout << "Matrix-Vector time: " << mvTimer.format() 
                  << "vector ops        : " << apTimer.format() 
                  << "preconditioner    : " << prTimer.format() << '\n';
      } 
    }
 
    virtual void apply (X& x, Xstar& b, double tol, Dune::InverseOperatorResult& res) 
    {
      terminate.tolerance(tol);
      (*this).apply(x,b,res);
    }
 
    virtual Dune::SolverCategory::Category category() const override
    {
      return Dune::SolverCategory::sequential;
    }
 
  private:
    // The solver may be constructed with either symmetric operator/preconditioner interface provided,
    // or with Dune interfaces provided. The implementation relies on the symmetric interface
    // in order to exploit performance gains with simultaneous evaluation of matrix-vector products
    // and dual pairings.
    // Therefore we hold "working" references to symmetric operator/preconditioners. Those reference
    // either the provided ones, or wrapper objects relaying the evaluation to provided Dune interfaces.
    // For the second case, we hold the (small) wrapper objects ourselves. In the first case, those
    // wrapper objects have to be initialized correctly (even though they are never accessed). For this 
    // we need a dual pairing, hence the dummy zeroDp.
    ZeroDualPairing<X,Xstar>                zeroDp;  // just a dummy
    
    SymmetricLinearOperatorWrapper<X,Xstar> proxyOp; // proxy in case a Dune::LinearOperator is provided
    SymmetricLinearOperator<X,Xstar>& op;  // working reference to operator
 
    SymmetricPreconditionerWrapper<X,Xstar>  proxyPrec; // proxy in case a Dune::Preconditioner is provided
    SymmetricPreconditioner<X,Xstar>&        prec;      // working reference to preconditioner
    
    PCGTerminationCriterion<Real>& terminate;
    int verbose;
    bool requiresInit;
    
    std::vector<field_type> resDebug;
  };
} // namespace Kaskade
#endif