iterative.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/base.hpp>
#include <kernel/util/dist.hpp>
#include <kernel/util/property_map.hpp>
#include <kernel/util/statistics.hpp>
 
namespace FEAT
{
  namespace Solver
  {
    namespace Intern
    {
      template<typename Matrix_>
      const Dist::Comm* get_mat_comm(const Matrix_& mat, typename Matrix_::GateRowType*)
      {
        return mat.get_comm();
      }
 
      template<typename Matrix_>
      const Dist::Comm* get_mat_comm(const Matrix_&, ...)
      {
        return nullptr;
      }
 
      template<typename Vector_>
      const Dist::Comm* get_vec_comm(const Vector_& vec, typename Vector_::GateType*)
      {
        return vec.get_comm();
      }
 
      template<typename Vector_>
      const Dist::Comm* get_vec_comm(const Vector_&, ...)
      {
        return nullptr;
      }
    } // namespace Intern
 
    enum class PlotMode
    {
      none = 0,
      iter,
      summary,
      all
    };
 
    inline std::ostream& operator<<(std::ostream& os, PlotMode mode)
    {
      switch(mode)
      {
        case PlotMode::none:
          return os << "none";
        case PlotMode::iter:
          return os << "iter";
        case PlotMode::summary:
          return os << "summary";
        case PlotMode::all:
          return os << "all";
        default:
          return os << "-unknown-";
      }
    }
 
    inline std::istream& operator>>(std::istream& is, PlotMode& mode)
    {
      String s;
      if((is >> s).fail())
        return is;
 
      if(s.compare_no_case("none") == 0)
        mode = PlotMode::none;
      else if(s.compare_no_case("iter") == 0)
        mode = PlotMode::iter;
      else if(s.compare_no_case("summary") == 0)
        mode = PlotMode::summary;
      else if(s.compare_no_case("all") == 0)
        mode = PlotMode::all;
      else
        is.setstate(std::ios_base::failbit);
 
      return is;
    }
 
    template<typename Vector_>
    class IterativeSolver :
      public SolverBase<Vector_>
    {
    public:
      typedef Vector_ VectorType;
      typedef typename VectorType::DataType DataType;
      typedef SolverBase<VectorType> BaseClass;
 
    protected:
      const Dist::Comm* _comm;
      String _plot_name;
      Status _status;
      DataType _tol_rel;
      DataType _tol_abs;
      DataType _tol_abs_low;
      DataType _div_rel;
      DataType _div_abs;
      DataType _stag_rate;
      Index _min_iter;
      Index _max_iter;
      Index _num_iter;
      Index _min_stag_iter;
      Index _num_stag_iter;
      DataType _def_init;
      DataType _def_cur;
      DataType _def_prev;
      Index _iter_digits;
      PlotMode _plot_mode;
      Index _plot_interval;
      bool _force_def_norm_calc;
 
      explicit IterativeSolver(const String& plot_name) :
        BaseClass(),
        _comm(nullptr),
        _plot_name(plot_name),
        _status(Status::undefined),
        _tol_rel(Math::sqrt(Math::eps<DataType>())),
        _tol_abs(DataType(1) / Math::sqr(Math::eps<DataType>())),
        _tol_abs_low(DataType(0)),
        _div_rel(DataType(1) / Math::eps<DataType>()),
        _div_abs(DataType(1) / Math::sqr(Math::eps<DataType>())),
        _stag_rate(DataType(0.95)),
        _min_iter(0),
        _max_iter(100),
        _num_iter(0),
        _min_stag_iter(0),
        _num_stag_iter(0),
        _def_init(0),
        _def_cur(0),
        _def_prev(0),
        _iter_digits(Math::ilog10(_max_iter)),
        _plot_mode(PlotMode::none),
        _plot_interval(1),
        _force_def_norm_calc(false)
      {
      }
 
      explicit IterativeSolver(const String& plot_name, const String& section_name, const PropertyMap* section) :
        IterativeSolver(plot_name)
      {
        auto plot_mode_p = section->get_entry("plot_mode");
        if (plot_mode_p.second && !plot_mode_p.first.parse(this->_plot_mode))
          throw ParseError(section_name + ".plot_mode", plot_mode_p.first, "one of: none, iter, summary, all");
 
        auto plot_name_p = section->get_entry("plot_name");
        if (plot_name_p.second)
          this->_plot_name = plot_name_p.first;
 
        auto plot_interval_p = section->get_entry("plot_interval");
        if (plot_interval_p.second && !plot_interval_p.first.parse(this->_plot_interval))
          throw ParseError(section_name + ".plot_interval", plot_interval_p.first, "an integer >= 0");
 
        auto tol_abs_p = section->get_entry("tol_abs");
        if (tol_abs_p.second && !tol_abs_p.first.parse(this->_tol_abs))
          throw ParseError(section_name + ".tol_abs", tol_abs_p.first, "a floating point number");
 
        auto tol_rel_p = section->get_entry("tol_rel");
        if (tol_rel_p.second && !tol_rel_p.first.parse(this->_tol_rel))
          throw ParseError(section_name + ".tol_rel", tol_rel_p.first, "a floating point number");
 
        auto tol_abs_low_p = section->get_entry("tol_abs_low");
        if (tol_abs_low_p.second && !tol_abs_low_p.first.parse(this->_tol_abs_low))
          throw ParseError(section_name + ".tol_abs_low", tol_abs_low_p.first, "a floating point number");
 
        auto div_abs_p = section->get_entry("div_abs");
        if (div_abs_p.second && !div_abs_p.first.parse(this->_div_abs))
          throw ParseError(section_name + ".div_abs", div_abs_p.first, "a floating point number");
 
        auto div_rel_p = section->get_entry("div_rel");
        if (div_rel_p.second && !div_rel_p.first.parse(this->_div_rel))
          throw ParseError(section_name + ".div_rel", div_rel_p.first, "a floating point number");
 
        auto stag_rate_p = section->get_entry("stag_rate");
        if (stag_rate_p.second && !stag_rate_p.first.parse(this->_stag_rate))
          throw ParseError(section_name + ".stag_rate", stag_rate_p.first, "a floating point number in range [0,1)");
 
        auto max_iter_p = section->get_entry("max_iter");
        if (max_iter_p.second && !max_iter_p.first.parse(this->_max_iter))
          throw ParseError(section_name + ".max_iter", max_iter_p.first, "an integer >= 0");
 
        auto min_iter_p = section->get_entry("min_iter");
        if (min_iter_p.second && !min_iter_p.first.parse(this->_min_iter))
          throw ParseError(section_name + ".min_iter", min_iter_p.first, "an integer >= 0");
 
        auto min_stag_iter_p = section->get_entry("min_stag_iter");
        if (min_stag_iter_p.second && !min_stag_iter_p.first.parse(this->_min_stag_iter))
          throw ParseError(section_name + ".min_stag_iter", min_stag_iter_p.first, "an integer >= 0");
      }
 
    public:
      void set_tol_rel(DataType tol_rel)
      {
        _tol_rel = tol_rel;
      }
 
      void set_tol_abs(DataType tol_abs)
      {
        _tol_abs = tol_abs;
      }
 
      void set_tol_abs_low(DataType tol_abs_low)
      {
        _tol_abs_low = tol_abs_low;
      }
 
      DataType get_tol_rel() const
      {
        return _tol_rel;
      }
 
      DataType get_tol_abs() const
      {
        return _tol_abs;
      }
 
      DataType get_tol_abs_low() const
      {
        return _tol_abs_low;
      }
 
      void set_div_rel(DataType div_rel)
      {
        _div_rel = div_rel;
      }
 
      void set_div_abs(DataType div_abs)
      {
        _div_abs = div_abs;
      }
 
      DataType get_div_rel() const
      {
        return _div_rel;
      }
 
      DataType get_div_abs() const
      {
        return _div_abs;
      }
 
      void set_stag_rate(DataType rate)
      {
        _stag_rate = rate;
      }
 
      DataType get_stag_rate() const
      {
        return _stag_rate;
      }
 
      void set_min_stag_iter(Index min_iter)
      {
        _min_stag_iter = min_iter;
      }
 
      Index get_min_stag_iter() const
      {
        return _min_stag_iter;
      }
 
      void set_min_iter(Index min_iter)
      {
        _min_iter = min_iter;
      }
 
      void set_max_iter(Index max_iter)
      {
        _max_iter = max_iter;
        _iter_digits = Math::ilog10(_max_iter);
      }
 
      Index get_num_iter() const
      {
        return _num_iter;
      }
 
      Index get_min_iter() const
      {
        return _min_iter;
      }
 
      Index get_max_iter() const
      {
        return _max_iter;
      }
 
      virtual void force_defect_norm_calc(bool force)
      {
        _force_def_norm_calc = force;
      }
 
      void set_plot_mode(const PlotMode plot_mode)
      {
        _plot_mode = plot_mode;
      }
 
      void set_plot_interval(const Index plot_interval)
      {
        _plot_interval = plot_interval;
      }
 
      void set_plot_name(const String& plot_name)
      {
        _plot_name = plot_name;
      }
 
      String get_plot_name() const
      {
        return _plot_name;
      }
 
      bool is_converged() const
      {
        return is_converged(_def_cur);
      }
 
      bool is_converged(const DataType def_cur) const
      {
        return (def_cur <= _tol_abs) && ((def_cur <= (_tol_rel * _def_init)) || (def_cur <= _tol_abs_low));
      }
 
      bool is_diverged() const
      {
        return is_diverged(_def_cur);
      }
 
      bool is_diverged(const DataType def_cur) const
      {
        return (def_cur > _div_abs) || (def_cur > (_div_rel * _def_init));
      }
 
      DataType get_def_initial() const
      {
        return _def_init;
      }
 
      DataType get_def_final() const
      {
        return _def_cur;
      }
 
      Status get_status() const
      {
        return _status;
      }
 
      DataType calc_convergence_rate() const
      {
        // no iterations performed?
        if(_num_iter <= Index(0))
          return DataType(0);
 
        // initial defect zero?
        if(_def_init < Math::eps<DataType>())
          return DataType(0);
 
        // compute convergence rate: (def_final / def_initial) ^ (1 / #iter)
        return Math::pow(_def_cur / _def_init, DataType(1) / DataType(_num_iter));
      }
 
      DataType calc_defect_reduction() const
      {
        // avoid division by zero
        if(this->_def_init <= Math::abs(this->_def_cur * Math::eps<DataType>()))
          return DataType(0);
 
        // compute defect reduction (def_final / def_inital)
        return this->_def_cur / this->_def_init;
      }
 
      virtual String get_summary() const
      {
        String msg;
        msg += String(  "Name............: ") + this->get_plot_name();
        msg += String("\nStatus..........: ") + stringify(this->get_status());
        msg += String("\nIterations......: ") + stringify(this->get_num_iter());
        msg += String("\nInitial Defect..: ") + stringify_fp_sci(this->get_def_initial());
        msg += String("\nFinal Defect....: ") + stringify_fp_sci(this->get_def_final());
        msg += String("\nDefect Reduction: ") + stringify_fp_sci(this->calc_defect_reduction());
        msg += String("\nConvergence Rate: ") + stringify_fp_fix(this->calc_convergence_rate());
 
        return msg;
      }
 
      virtual void plot_summary() const
      {
        // Print solver summary
        if(!_plot_summary())
          return;
 
        this->_print_line(this->get_summary());
      }
 
      virtual Status correct(VectorType& vec_sol, const VectorType& vec_rhs) = 0;
 
    protected:
      void _set_comm(const Dist::Comm* comm)
      {
        this->_comm = comm;
      }
 
      template<typename Matrix_>
      void _set_comm_by_matrix(const Matrix_& matrix)
      {
        this->_comm = Intern::get_mat_comm(matrix, nullptr);
      }
 
      void _set_comm_by_vector(const Vector_& vector)
      {
        this->_comm = Intern::get_vec_comm(vector, nullptr);
      }
 
      virtual DataType _calc_def_norm(const VectorType& vec_def, const VectorType& DOXY(vec_sol))
      {
        return vec_def.norm2();
      }
 
      bool _plot_iter(Status st = Status::progress) const
      {
        return ((_num_iter % _plot_interval == 0) || (st != Status::progress))
          && ((_plot_mode == PlotMode::iter) || (_plot_mode == PlotMode::all));
      }
 
      bool _plot_summary() const
      {
        return _plot_mode == PlotMode::summary || _plot_mode == PlotMode::all;
      }
 
      bool _progress() const
      {
        return this->_status == Status::progress;
      }
 
      void _print_line(const String& line) const
      {
        // print message line via comm (if available)
        if(_comm != nullptr)
          _comm->print(line);
        else
          std::cout << line << "\n";
      }
 
      virtual void _plot_iter_line(Index num_iter, DataType def_cur, DataType def_prev)
      {
        // compose message line
        String msg = this->_plot_name
          +  ": " + stringify(num_iter).pad_front(this->_iter_digits)
          + " : " + stringify_fp_sci(def_cur);
 
        // not first iteration?
        if(num_iter > Index(0))
        {
          msg += " / " + stringify_fp_sci(def_cur / this->_def_init);
          msg += " / " + stringify_fp_fix(def_cur / def_prev);
        }
 
        this->_print_line(msg);
      }
 
      virtual Status _set_initial_defect(const VectorType& vec_def, const VectorType& vec_sol)
      {
        // reset iteration variables
        this->_def_init = this->_def_prev = this->_def_cur = DataType(0);
        this->_num_iter = Index(0);
        this->_num_stag_iter = Index(0);
 
        // first, let's see if we have to compute the defect at all
        bool calc_def = _force_def_norm_calc
          || (this->_min_iter < this->_max_iter)
          || (this->_plot_iter() || this->_plot_summary())
          || (this->_min_stag_iter > Index(0));
 
        // if we don't have to compute the initial defect, we can exit here
        if(!calc_def)
          return Status::progress;
 
        // compute new initial defect
        this->_def_init = this->_def_cur = this->_def_prev = this->_calc_def_norm(vec_def, vec_sol);
        Statistics::add_solver_expression(std::make_shared<ExpressionDefect>(this->name(), this->_def_init, this->get_num_iter()));
 
        // plot iteration line?
        if(this->_plot_iter())
          this->_plot_iter_line(Index(0), this->_def_init, this->_def_init);
 
        // ensure that the defect is neither NaN nor infinity
        if(!Math::isfinite(this->_def_init))
          return Status::aborted;
 
        // check against low absolute tolerance; this may be fulfilled if we have only a handful of dofs
        if(this->_def_init < this->_tol_abs_low)
          return Status::success;
 
        // check if the initial defect is zero; we test against eps^2 here
        if(this->_def_init <= Math::sqr(Math::eps<DataType>()))
          return Status::success;
 
        // continue iterating
        return Status::progress;
      }
 
      virtual Status _analyze_defect(Index num_iter, DataType def_cur, DataType def_prev, bool check_stag)
      {
        // ensure that the defect is neither NaN nor infinity
        if(!Math::isfinite(def_cur))
          return Status::aborted;
 
        // is diverged?
        if(this->is_diverged(def_cur))
          return Status::diverged;
 
        // maximum number of iterations performed?
        if(num_iter >= this->_max_iter)
          return Status::max_iter;
 
        // minimum number of iterations performed?
        if(num_iter < this->_min_iter)
          return Status::progress;
 
        // is converged?
        if(this->is_converged(def_cur))
          return Status::success;
 
        // check for stagnation?
        if(check_stag && (this->_min_stag_iter > Index(0)))
        {
          // did this iteration stagnate?
          if(def_cur >= this->_stag_rate * def_prev)
          {
            // increment stagnation count
            if(++this->_num_stag_iter >= this->_min_stag_iter)
              return Status::stagnated;
          }
          else
          {
            // this iteration did not stagnate
            this->_num_stag_iter = Index(0);
          }
        }
 
        // continue iterating
        return Status::progress;
      }
 
      virtual Status _set_new_defect(const VectorType& vec_def, const VectorType& vec_sol)
      {
        // increase iteration count
        ++this->_num_iter;
 
        // store previous defect
        this->_def_prev = this->_def_cur;
 
        // first, let's see if we have to compute the defect at all
        bool calc_def = _force_def_norm_calc
          || (this->_min_iter < this->_max_iter)
          || (this->_plot_iter() || this->_plot_summary())
          || (this->_min_stag_iter > Index(0));
 
        // compute new defect
        if(calc_def)
        {
          this->_def_cur = this->_calc_def_norm(vec_def, vec_sol);
          Statistics::add_solver_expression(std::make_shared<ExpressionDefect>(this->name(), this->_def_cur, this->get_num_iter()));
        }
 
        // analyze defect
        Status status = this->_analyze_defect(this->_num_iter, this->_def_cur, this->_def_prev, true);
 
        // plot defect?
        if(this->_plot_iter(status))
          this->_plot_iter_line(this->_num_iter, this->_def_cur, this->_def_prev);
 
        // return status
        return status;
      }
 
      virtual Status _update_defect(const DataType def_cur_norm)
      {
        // increase iteration count
        ++this->_num_iter;
 
        // store previous defect
        this->_def_prev = this->_def_cur;
 
        // update current defect
        this->_def_cur = def_cur_norm;
        Statistics::add_solver_expression(std::make_shared<ExpressionDefect>(this->name(), this->_def_cur, this->get_num_iter()));
 
        // analyze defect
        Status status = this->_analyze_defect(this->_num_iter, this->_def_cur, this->_def_prev, true);
 
        // plot defect?
        if(this->_plot_iter(status))
          this->_plot_iter_line(this->_num_iter, this->_def_cur, this->_def_prev);
 
        // return status
        return status;
      }
    }; // class IterativeSolver
 
    template<
      typename Vector_,
      typename Matrix_,
      typename Filter_>
    inline Status solve(
      IterativeSolver<Vector_>& solver,
      Vector_& vec_sol,
      const Vector_& vec_rhs,
      const Matrix_& DOXY(matrix),
      const Filter_& DOXY(filter))
    {
      // simply call the 'correct' method
      return solver.correct(vec_sol, vec_rhs);
    }
 
    template<typename Vector_>
    class PreconditionedIterativeSolver :
      public IterativeSolver<Vector_>
    {
    public:
      typedef Vector_ VectorType;
      typedef typename VectorType::DataType DataType;
      typedef IterativeSolver<VectorType> BaseClass;
      typedef SolverBase<VectorType> PrecondType;
 
    protected:
      std::shared_ptr<PrecondType> _precond;
 
      explicit PreconditionedIterativeSolver(const String& plot_name,std::shared_ptr<PrecondType> precond = nullptr) :
        BaseClass(plot_name),
        _precond(precond)
      {
      }
 
      explicit PreconditionedIterativeSolver(const String& plot_name, const String& section_name,
        const PropertyMap* section, std::shared_ptr<PrecondType> precond = nullptr) :
        BaseClass(plot_name, section_name, section),
        _precond(precond)
      {
      }
 
    public:
      virtual ~PreconditionedIterativeSolver()
      {
      }
 
      virtual void init_symbolic() override
      {
        BaseClass::init_symbolic();
        if(_precond)
          _precond->init_symbolic();
      }
 
      virtual void init_numeric() override
      {
        BaseClass::init_numeric();
        if(_precond)
          _precond->init_numeric();
      }
 
      virtual void done_numeric() override
      {
        if(_precond)
          _precond->done_numeric();
        BaseClass::done_numeric();
      }
 
      virtual void done_symbolic() override
      {
        if(_precond)
          _precond->done_symbolic();
        BaseClass::done_symbolic();
      }
 
    protected:
      template<typename Filter_>
      bool _apply_precond(VectorType& vec_cor, const VectorType& vec_def, const Filter_& filter)
      {
        XASSERTM(&vec_cor != &vec_def, "vec_cor and vec_def must not refer to the same object");
        if(this->_precond)
        {
          Statistics::add_solver_expression(std::make_shared<ExpressionCallPrecond>(this->name(), this->_precond->name()));
          return status_success(this->_precond->apply(vec_cor, vec_def));
        }
        else
        {
          vec_cor.copy(vec_def);
          filter.filter_cor(vec_cor);
          return true;
        }
      }
    }; // class PreconditionedIterativeSolver<...>
  } // namespace Solver
} // namespace FEAT