fgmres.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/solver/iterative.hpp>
 
// includes, FEAT
#include <vector>
 
namespace FEAT
{
  namespace Solver
  {
    template<
      typename Matrix_,
      typename Filter_>
    class FGMRES :
      public PreconditionedIterativeSolver<typename Matrix_::VectorTypeR>
    {
    public:
      typedef Matrix_ MatrixType;
      typedef Filter_ FilterType;
      typedef typename MatrixType::VectorTypeR VectorType;
      typedef typename MatrixType::DataType DataType;
      typedef PreconditionedIterativeSolver<VectorType> BaseClass;
 
      typedef SolverBase<VectorType> PrecondType;
 
    protected:
      const MatrixType& _system_matrix;
      const FilterType& _system_filter;
      Index _krylov_dim;
      DataType _inner_res_scale;
      std::vector<VectorType> _vec_v, _vec_z;
      std::vector<DataType> _c, _s, _q;
      std::vector<std::vector<DataType>> _h;
 
    public:
      explicit FGMRES(const MatrixType& matrix, const FilterType& filter, Index krylov_dim,
        DataType inner_res_scale = DataType(0), std::shared_ptr<PrecondType> precond = nullptr) :
        BaseClass("FGMRES(" + stringify(krylov_dim) + ")", precond),
        _system_matrix(matrix),
        _system_filter(filter),
        _krylov_dim(krylov_dim)
      {
        // we always need to compute defect norms
        this->_force_def_norm_calc = true;
 
        // set communicator by system matrix
        this->_set_comm_by_matrix(matrix);
        set_inner_res_scale(inner_res_scale);
      }
 
      explicit FGMRES(const String& section_name, const PropertyMap* section,
        const MatrixType& matrix, const FilterType& filter, std::shared_ptr<PrecondType> precond = nullptr) :
        BaseClass("FGMRES", section_name, section, precond),
        _system_matrix(matrix),
        _system_filter(filter),
        _inner_res_scale(DataType(0))
      {
        // we always need to compute defect norms
        this->_force_def_norm_calc = true;
 
        // set communicator by system matrix
        this->_set_comm_by_matrix(matrix);
 
        // Set _inner_res_scale by parameter or use default value
        auto inner_res_scale_p = section->query("inner_res_scale");
        if(inner_res_scale_p.second && (!inner_res_scale_p.first.parse(this->_inner_res_scale) || (this->_inner_res_scale < DataType(0))))
          throw ParseError(section_name + ".inner_res_scale", inner_res_scale_p.first, "a non-negative float");
 
        // Check if we have set _krylov_vim
        auto krylov_dim_p = section->query("krylov_dim");
        if(!krylov_dim_p.second)
          throw ParseError("FGMRES config section is missing the mandatory krylov_dim!");
 
        if(!krylov_dim_p.first.parse(this->_krylov_dim) || (this->_krylov_dim <= Index(0)))
          throw ParseError(section_name + ".krylov_dim", krylov_dim_p.first, "a positive integer");
 
        this->set_plot_name("FGMRES("+stringify(_krylov_dim)+")");
      }
 
      virtual ~FGMRES()
      {
      }
 
      virtual String name() const override
      {
        return "FGMRES";
      }
 
      virtual void force_defect_norm_calc(bool DOXY(force)) override
      {
        // force is already applied in constructor
      }
 
      virtual void init_symbolic() override
      {
        BaseClass::init_symbolic();
 
        _c.reserve(_krylov_dim);
        _s.reserve(_krylov_dim);
        _q.reserve(_krylov_dim + 1);
        _h.resize(_krylov_dim);
 
        for(Index i(0); i < _krylov_dim; ++i)
        {
          _h.at(i).resize(i+1);
        }
 
        _vec_v.push_back(this->_system_matrix.create_vector_r());
        for(Index i(0); i < _krylov_dim; ++i)
        {
          _vec_v.push_back(this->_vec_v.front().clone(LAFEM::CloneMode::Layout));
          _vec_z.push_back(this->_vec_v.front().clone(LAFEM::CloneMode::Layout));
        }
      }
 
      virtual void done_symbolic() override
      {
        _vec_v.clear();
        _vec_z.clear();
        BaseClass::done_symbolic();
      }
 
      virtual void set_krylov_dim(Index krylov_dim)
      {
        XASSERT(krylov_dim > Index(0));
        _krylov_dim = krylov_dim;
      }
 
      virtual void set_inner_res_scale(DataType inner_res_scale)
      {
        XASSERT(inner_res_scale >= DataType(0));
        _inner_res_scale = inner_res_scale;
      }
 
      virtual Status apply(VectorType& vec_sol, const VectorType& vec_rhs) override
      {
        // save input rhs vector as initial defect
        this->_vec_v.at(0).copy(vec_rhs);
 
        // clear solution vector
        vec_sol.format();
 
        // apply
        this->_status = _apply_intern(vec_sol, vec_rhs);
        this->plot_summary();
        return this->_status;
      }
 
      virtual Status correct(VectorType& vec_sol, const VectorType& vec_rhs) override
      {
        // compute initial defect
        this->_system_matrix.apply(this->_vec_v.at(0), vec_sol, vec_rhs, -DataType(1));
        this->_system_filter.filter_def(this->_vec_v.at(0));
 
        // apply
        this->_status = _apply_intern(vec_sol, vec_rhs);
        this->plot_summary();
        return this->_status;
      }
 
    protected:
      virtual Status _apply_intern(VectorType& vec_sol, const VectorType& vec_rhs)
      {
        IterationStats pre_iter(*this);
        Statistics::add_solver_expression(std::make_shared<ExpressionStartSolve>(this->name()));
        const MatrixType& matrix(this->_system_matrix);
        const FilterType& filter(this->_system_filter);
 
        // compute initial defect
        Status status = this->_set_initial_defect(this->_vec_v.at(0), vec_sol);
 
        pre_iter.destroy();
 
        // outer GMRES loop
        while(status == Status::progress)
        {
          IterationStats stat(*this);
 
          _q.clear();
          _s.clear();
          _c.clear();
          _q.push_back(this->_def_cur);
 
          // normalize v[0]
          this->_vec_v.at(0).scale(this->_vec_v.at(0), DataType(1) / _q.back());
 
          // inner GMRES loop
          Index i(0);
          while(i < this->_krylov_dim)
          {
            // apply preconditioner
            if(!this->_apply_precond(this->_vec_z.at(i), this->_vec_v.at(i), filter))
            {
              stat.destroy();
              Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::aborted, this->get_num_iter()));
              return Status::aborted;
            }
            //filter.filter_cor(this->_vec_z.at(i));
 
            // v[i+1] := A*z[i]
            matrix.apply(this->_vec_v.at(i+1), this->_vec_z.at(i));
            filter.filter_def(this->_vec_v.at(i+1));
 
            // Gram-Schmidt process
            for(Index k(0); k <= i; ++k)
            {
              this->_h.at(i).at(k) = this->_vec_v.at(i+1).dot(this->_vec_v.at(k));
              this->_vec_v.at(i+1).axpy(this->_vec_v.at(k), -this->_h.at(i).at(k));
            }
 
            // normalize v[i+1]
            DataType alpha = this->_vec_v.at(i+1).norm2();
            this->_vec_v.at(i+1).scale(this->_vec_v.at(i+1), DataType(1) / alpha);
 
            // apply Givens rotations
            for(Index k(0); k < i; ++k)
            {
              DataType t(this->_h.at(i).at(k));
              this->_h.at(i).at(k  ) = this->_c.at(k) * t + this->_s.at(k) * this->_h.at(i).at(k+1);
              this->_h.at(i).at(k+1) = this->_s.at(k) * t - this->_c.at(k) * this->_h.at(i).at(k+1);
            }
 
            // compute beta
            DataType beta = Math::sqrt(Math::sqr(this->_h.at(i).at(i)) + Math::sqr(alpha));
 
            // compute next plane rotation
            _s.push_back(alpha / beta);
            _c.push_back(this->_h.at(i).at(i) / beta);
 
            this->_h.at(i).at(i) = beta;
            this->_q.push_back(this->_s.back() * this->_q.at(i));
            this->_q.at(i) *= this->_c.back();
 
            // push our new defect
            if(++i < this->_krylov_dim)
            {
              // get the absolute defect
              DataType def_cur = Math::abs(this->_q.back());
 
              // did we diverge?
              if((def_cur > this->_div_abs) || (def_cur > (this->_div_rel * this->_def_init)))
                break;
 
              // increase iteration count
              ++this->_num_iter;
 
              // minimum number of iterations performed?
              if(!(this->_num_iter < this->_min_iter))
              {
                // did we converge?
                if((def_cur <= _inner_res_scale * this->_tol_abs) &&
                  ((def_cur <= _inner_res_scale * (this->_tol_rel * this->_def_init)) || (def_cur <= _inner_res_scale * this->_tol_abs_low) ))
                  break;
 
                // maximum number of iterations performed?
                if(this->_num_iter >= this->_max_iter)
                  break;
              }
 
              // set our pseudo defect
              // compute new defect
              this->_def_cur = def_cur;
 
              // plot?
              if(this->_plot_iter())
              {
                String msg = this->_plot_name
                  +  "* " + stringify(this->_num_iter).pad_front(this->_iter_digits)
                  + " : " + stringify_fp_sci(this->_def_cur);
                this->_print_line(msg);
              }
            }
          }
 
          Index n = Math::min(i, this->_krylov_dim);
 
          // solve H*q = q
          for(Index k(n); k > 0;)
          {
            --k;
            this->_q.at(k) /= this->_h.at(k).at(k);
            for(Index j(k); j > 0;)
            {
              --j;
              this->_q.at(j) -= this->_h.at(k).at(j) * this->_q.at(k);
            }
          }
 
          // update solution
          for(Index k(0); k < n; ++k)
            vec_sol.axpy(this->_vec_z.at(k), this->_q.at(k));
 
          // compute "real" residual
          matrix.apply(this->_vec_v.at(0), vec_sol, vec_rhs, -DataType(1));
          filter.filter_def(this->_vec_v.at(0));
 
          // set the current defect
          status = this->_set_new_defect(this->_vec_v.at(0), vec_sol);
        }
 
        // finished
        Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), status, this->get_num_iter()));
        return status;
      }
    }; // class FGMRES<...>
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<FGMRES<Matrix_, Filter_>> new_fgmres(
      const Matrix_& matrix, const Filter_& filter, Index krylov_dim,
      typename Matrix_::DataType inner_res_scale = typename Matrix_::DataType(0),
      std::shared_ptr<SolverBase<typename Matrix_::VectorTypeL>> precond = nullptr)
    {
      return std::make_shared<FGMRES<Matrix_, Filter_>>(matrix, filter, krylov_dim, inner_res_scale, precond);
    }
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<FGMRES<Matrix_, Filter_>> new_fgmres(
      const String& section_name, const PropertyMap* section,
      const Matrix_& matrix, const Filter_& filter,
      std::shared_ptr<SolverBase<typename Matrix_::VectorTypeL>> precond = nullptr)
    {
      return std::make_shared<FGMRES<Matrix_, Filter_>>(section_name, section, matrix, filter, precond);
    }
  } // namespace Solver
} // namespace FEAT