chebyshev.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/base_header.hpp>
#include <kernel/solver/iterative.hpp>
 
namespace FEAT
{
  namespace Solver
  {
    template<typename Matrix_, typename Filter_>
    class Chebyshev :
      public IterativeSolver<typename Matrix_::VectorTypeR>
    {
    public:
      typedef Matrix_ MatrixType;
      typedef Filter_ FilterType;
      typedef typename MatrixType::VectorTypeL VectorType;
      typedef typename MatrixType::DataType DataType;
      typedef IterativeSolver<VectorType> BaseClass;
 
    protected:
      const MatrixType& _system_matrix;
      const FilterType& _system_filter;
      VectorType _vec_def;
      VectorType _vec_cor;
      DataType _min_ev;
      DataType _max_ev;
      DataType _fraction_min_ev;
      DataType _fraction_max_ev;
 
    public:
      explicit Chebyshev(const MatrixType& matrix, const FilterType& filter,
        const DataType fraction_min_ev = DataType(0.5), const DataType fraction_max_ev = DataType(0.8)) :
        BaseClass("Chebyshev"),
        _system_matrix(matrix),
        _system_filter(filter),
        _min_ev(0),
        _max_ev(0),
        _fraction_min_ev(fraction_min_ev),
        _fraction_max_ev(fraction_max_ev)
      {
        // set communicator by system matrix
        this->_set_comm_by_matrix(matrix);
      }
 
      explicit Chebyshev(const String& section_name, const PropertyMap* section,
        const MatrixType& matrix, const FilterType& filter) :
        BaseClass("Chebyshev", section_name, section),
        _system_matrix(matrix),
        _system_filter(filter),
        _min_ev(DataType(0)),
        _max_ev(DataType(0)),
        _fraction_min_ev(DataType(0.5)),
        _fraction_max_ev(DataType(0.8))
      {
        // set communicator by system matrix
        this->_set_comm_by_matrix(matrix);
        auto fmin_p = section->query("fraction_min_ev");
        if(fmin_p.second && !fmin_p.first.parse(this->_fraction_min_ev))
          throw ParseError(section_name + ".fraction_min_ev", fmin_p.first, "a float");
 
        auto fmax_p = section->query("fraction_max_ev");
        if(fmax_p.second && !fmax_p.first.parse(this->_fraction_max_ev))
          throw ParseError(section_name + ".fraction_max_ev", fmax_p.first, "a float");
      }
 
      virtual ~Chebyshev()
      {
      }
 
      virtual String name() const override
      {
        return "Chebyshev";
      }
 
      void set_fraction_min_ev(DataType fraction_min_ev)
      {
        _fraction_min_ev = fraction_min_ev;
      }
 
      void set_fraction_max_ev(DataType fraction_max_ev)
      {
        _fraction_max_ev = fraction_max_ev;
      }
 
      virtual void init_symbolic() override
      {
        BaseClass::init_symbolic();
        _vec_def = this->_system_matrix.create_vector_r();
        _vec_cor = this->_system_matrix.create_vector_r();
      }
 
      virtual void done_symbolic() override
      {
        this->_vec_cor.clear();
        this->_vec_def.clear();
        BaseClass::done_symbolic();
      }
 
      virtual void init_numeric() override
      {
        BaseClass::init_numeric();
 
        //PowerMethod for maximum eigenvalue
        //http://blogs.sas.com/content/iml/2012/05/09/the-power-method.html
        //https://en.wikipedia.org/wiki/Power_iteration
        VectorType& v(this->_vec_cor);
        VectorType& z(this->_vec_def);
        const MatrixType& matrix(this->_system_matrix);
        v.format(DataType(3));
        DataType tolerance(DataType(1e-4));
        Index max_iters(40);
        v.scale(v, DataType(1) / v.norm2());
        DataType lambda_old(0), lambda(0);
        for (Index i(0) ; i < max_iters ; ++i)
        {
          matrix.apply(z, v);
          v.scale(z, DataType(1) / z.norm2());
          lambda = v.dot(z);
          if (Math::abs((lambda - lambda_old) / lambda) < tolerance)
            break;
          lambda_old = lambda;
        }
 
        _min_ev = lambda * _fraction_min_ev;
        _max_ev = lambda * _fraction_max_ev;
      }
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
      {
        // save defect
        this->_vec_def.copy(vec_def);
 
        // clear solution vector
        vec_cor.format();
 
        // apply
        this->_status = _apply_intern(vec_cor, vec_def);
        this->plot_summary();
        return this->_status;
      }
 
      virtual Status correct(VectorType& vec_sol, const VectorType& vec_rhs) override
      {
        // compute defect
        this->_system_matrix.apply(this->_vec_def, vec_sol, vec_rhs, -DataType(1));
        this->_system_filter.filter_def(this->_vec_def);
 
        // 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)
      {
        Statistics::add_solver_expression(std::make_shared<ExpressionStartSolve>(this->name()));
 
        VectorType& vec_def(this->_vec_def);
        VectorType& vec_cor(this->_vec_cor);
        const MatrixType& matrix(this->_system_matrix);
        const FilterType& filter(this->_system_filter);
 
        // compute initial defect
        Status status = this->_set_initial_defect(vec_def, vec_sol);
 
        DataType d = (_max_ev + _min_ev) / DataType(2);
        DataType c = (_max_ev - _min_ev) / DataType(2);
        DataType alpha(0);
        DataType beta(0);
        vec_cor.scale(vec_def, (DataType(1) / d));
 
        // start iterating
        while(status == Status::progress)
        {
          switch (this->get_num_iter())
          {
          case 0:
            alpha = DataType(1) / d;
            break;
          case 1:
            alpha = DataType(2) * d * (DataType(1) / ((DataType(2) * d * d) - (c * c)));
            break;
          default:
            alpha = DataType(1) / (d - ((alpha * c * c) / DataType(4)));
          }
 
          beta = (alpha * d) - DataType(1);
          vec_cor.scale(vec_cor, beta);
          vec_cor.axpy(vec_def, alpha);
 
          // update solution vector
          vec_sol.axpy(vec_cor);
 
          // compute new defect vector
          matrix.apply(vec_def, vec_sol, vec_rhs, -DataType(1));
          filter.filter_def(vec_def);
 
          // compute new defect norm
          status = this->_set_new_defect(vec_def, vec_sol);
        }
 
        // return our status
        Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), status, this->get_num_iter()));
        return status;
      }
    }; // class Chebyshev<...>
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<Chebyshev<Matrix_, Filter_>> new_chebyshev(
      const Matrix_& matrix, const Filter_& filter,
      const typename Matrix_::DataType fraction_min_ev = typename Matrix_::DataType(0.03),
      const typename Matrix_::DataType fraction_max_ev = typename Matrix_::DataType(1.1))
    {
      return std::make_shared<Chebyshev<Matrix_, Filter_>>(matrix, filter, fraction_min_ev, fraction_max_ev);
    }
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<Chebyshev<Matrix_, Filter_>> new_chebyshev(
      const String& section_name, const PropertyMap* section,
      const Matrix_& matrix, const Filter_& filter)
    {
      return std::make_shared<Chebyshev<Matrix_, Filter_>>(section_name, section, matrix, filter);
    }
  } // namespace Solver
} // namespace FEAT