polynomial_precond.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/base.hpp>
 
namespace FEAT
{
  namespace Solver
  {
    template<typename Matrix_, typename Filter_>
    class PolynomialPrecond :
      public SolverBase<typename Matrix_::VectorTypeL>
    {
    public:
      typedef Matrix_ MatrixType;
      typedef Filter_ FilterType;
      typedef typename MatrixType::VectorTypeL VectorType;
      typedef typename MatrixType::DataType DataType;
      typedef SolverBase<VectorType> BaseClass;
 
    protected:
      const MatrixType& _matrix;
      const FilterType& _filter;
      Index _m;
      DataType _omega;
      VectorType _inv_diag;
      VectorType _aux1, _aux2, _aux3;
 
    public:
      explicit PolynomialPrecond(const MatrixType& matrix, const FilterType& filter, Index m, DataType omega = DataType(1)) :
        _matrix(matrix),
        _filter(filter),
        _m(m),
        _omega(omega)
      {
      }
 
      explicit PolynomialPrecond(const String& section_name, const PropertyMap* section,
      const MatrixType& matrix, const FilterType& filter) :
        BaseClass(section_name, section),
        _matrix(matrix),
        _filter(filter),
        _m(1),
        _omega(1)
        {
          auto omega_p = section->query("omega");
          if(omega_p.second)
          {
            set_omega(DataType(std::stod(omega_p.first)));
          }
          else
          {
            XABORTM(name()+" config section is missing the mandatory omega key!");
          }
 
          auto m_p = section->query("m");
          if(m_p.second)
          {
            set_m(std::stoul(m_p.first));
          }
          else
          {
            XABORTM(name()+" config section is missing the mandatory m key!");
          }
        }
 
      virtual ~PolynomialPrecond()
      {
      }
 
      virtual String name() const override
      {
        return "Polynomial";
      }
 
      virtual void init_symbolic() override
      {
        _inv_diag = _matrix.create_vector_r();
        _aux1 = _matrix.create_vector_r();
        _aux2 = _matrix.create_vector_r();
        _aux3 = _matrix.create_vector_r();
      }
 
      virtual void done_symbolic() override
      {
        _inv_diag.clear();
        _aux1.clear();
        _aux2.clear();
        _aux3.clear();
      }
 
      virtual void init_numeric() override
      {
        // extract matrix diagonal
        _matrix.extract_diag(_inv_diag);
 
        // invert diagonal elements
        _inv_diag.component_invert(_inv_diag, _omega);
      }
 
      void set_omega(DataType omega)
      {
        XASSERT(omega > DataType(0));
        _omega = omega;
      }
 
      void set_m(Index m)
      {
        XASSERT(m > Index(0));
        _m = m;
      }
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
      {
        /*
         * preconditioner is given by
         *   \f$ M^-1 = \left(I + (I - \tilde M^{-1}A) + ... + (I - \tilde M^{-1 }A)^m\right) \tilde M^{-1} \f$
         *
         * the preconditioner only works, if
         *   ||I - \tilde M^{-1} A||_2 < 1.
         */
 
        vec_cor.component_product(_inv_diag, vec_def);
        _aux3.copy(vec_cor);
 
        for (Index i = 1; i <= _m; ++i)
        {
          _matrix.apply(_aux1, vec_cor);
          _filter.filter_def(_aux1);
          _aux2.component_product(_inv_diag, _aux1);
          vec_cor.axpy(_aux3);
          vec_cor.axpy(_aux2, DataType(-1.0));
        }
 
        this->_filter.filter_cor(vec_cor);
 
        return Status::success;
      }
    }; // class PolynomialPrecond<...>
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<PolynomialPrecond<Matrix_, Filter_>> new_polynomial_precond(
      const Matrix_& matrix, const Filter_& filter,
      const Index m, const typename Matrix_::DataType omega = typename Matrix_::DataType(1))
    {
      return std::make_shared<PolynomialPrecond<Matrix_, Filter_>>(matrix, filter, m, omega);
    }
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<PolynomialPrecond<Matrix_, Filter_>> new_polynomial_precond(
      const String& section_name, const PropertyMap* section,
      const Matrix_& matrix, const Filter_& filter)
    {
      return std::make_shared<PolynomialPrecond<Matrix_, Filter_>>(section_name, section, matrix, filter);
    }
  } // namespace Solver
} // namespace FEAT