bicgstabl.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>
 
namespace FEAT
{
  namespace Solver
  {
 
    enum class BiCGStabLPreconVariant
    {
      left,
      right
    };
 
 
    inline std::ostream& operator<<(std::ostream& os, BiCGStabLPreconVariant precon_variant)
    {
      switch(precon_variant)
      {
        case BiCGStabLPreconVariant::left:
          return os << "left";
        case BiCGStabLPreconVariant::right:
          return os << "right";
        default:
          return os << "-unknown-";
      }
    }
 
    inline std::istream& operator>>(std::istream& is, BiCGStabLPreconVariant& precon_variant)
    {
      String s;
      if((is >> s).fail())
        return is;
 
      if(s.compare_no_case("left") == 0)
        precon_variant = BiCGStabLPreconVariant::left;
      else if(s.compare_no_case("right") == 0)
        precon_variant = BiCGStabLPreconVariant::right;
      else
        is.setstate(std::ios_base::failbit);
 
      return is;
    }
 
 
    template<typename Matrix_, typename Filter_>
    class BiCGStabL :
      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;
 
        VectorType _vec_r0_tilde;
        VectorType _vec_u0;
        VectorType _vec_pc;
 
        std::vector<VectorType> _vec_uj_hat;
 
        std::vector<VectorType> _vec_rj_hat;
 
        BiCGStabLPreconVariant _precon_variant;
 
        int _l;
 
      public:
        explicit BiCGStabL(const MatrixType& matrix, const FilterType& filter,
          int l = 2, //default parameter 2 (1 would lead to standard BiCGStab algorithm)
          std::shared_ptr<PrecondType> precond = nullptr,
          BiCGStabLPreconVariant precon_variant = BiCGStabLPreconVariant::left)
           :
          BaseClass("BiCGStab(" + stringify(l) +")", precond),
          _system_matrix(matrix),
          _system_filter(filter),
          _precon_variant(precon_variant),
          _l(l)
        {
          // we always need to compute defect norms
          this->_force_def_norm_calc = true;
 
          // set communicator by system matrix
          this->_set_comm_by_matrix(matrix);
          XASSERT(precon_variant == BiCGStabLPreconVariant::left || precon_variant == BiCGStabLPreconVariant::right);
          XASSERTM(_l > 0, "bicgstabl polynomial degree must be greater than zero!");
        }
 
 
        explicit BiCGStabL(const String& section_name, const PropertyMap* section,
          const MatrixType& matrix, const FilterType& filter,
          std::shared_ptr<PrecondType> precond = nullptr)
           :
          BaseClass("BiCGStabL", section_name, section, precond),
          _system_matrix(matrix),
          _system_filter(filter),
          _precon_variant(BiCGStabLPreconVariant::left)
        {
          // we always need to compute defect norms
          this->_force_def_norm_calc = true;
 
          // set communicator by system matrix
          this->_set_comm_by_matrix(matrix);
          // Check if we have set _p_variant
          auto p_variant_p = section->query("precon_variant");
          if(p_variant_p.second && !p_variant_p.first.parse(this->_precon_variant))
            throw ParseError(section_name + ".precon_variant", p_variant_p.first, "one of: left, right");
 
          //Check if we have set the solver parameter _l
          std::pair<String , bool> l_pair = section->query("polynomial_degree");
 
          int l_pm = 2; //use default l, if the parameter is not given
          if (l_pair.second && (!l_pair.first.parse(l_pm) || (l_pm <= 0)))
            throw ParseError(section_name + ".polynomial_degree", l_pair.first, "a positive integer");
 
          this->_l = l_pm;
          this->set_plot_name("BiCGStab(" + stringify(_l) + ")");
        }
 
        virtual String name() const override
        {
          return "BiCGStabL";
        }
 
        virtual void force_defect_norm_calc(bool DOXY(force)) override
        {
          // force is already applied in constructor
        }
 
        virtual void init_symbolic() override
        {
          BaseClass::init_symbolic();
 
          // create all temporary vectors
          _vec_r0_tilde = this->_system_matrix.create_vector_r();
          _vec_u0       = this->_system_matrix.create_vector_r();
          _vec_pc       = this->_system_matrix.create_vector_r();
 
          _vec_uj_hat.resize(Index(_l+1));
          _vec_rj_hat.resize(Index(_l+1));
          for (int i=0; i <= _l; i++)
          {
            _vec_rj_hat.at(Index(i)) = this->_system_matrix.create_vector_r();
            _vec_uj_hat.at(Index(i)) = this->_system_matrix.create_vector_r();
          }
        }
 
        virtual void done_symbolic() override
        {
          _vec_r0_tilde.clear();
          _vec_u0.clear();
          _vec_pc.clear();
          for (int i=0; i <= _l; i++)
          {
            _vec_rj_hat.at(Index(i)).clear();
            _vec_uj_hat.at(Index(i)).clear();
          }
          BaseClass::done_symbolic();
        }
 
        virtual Status correct(VectorType& vec_sol, const VectorType& vec_rhs) override
        {
          // compute initial defect
          this->_system_matrix.apply(this->_vec_rj_hat.at(0), vec_sol, vec_rhs, -DataType(1));
          this->_system_filter.filter_def(this->_vec_rj_hat.at(0));
 
          // apply solver
          this->_status = _apply_intern(vec_sol, vec_rhs);
 
          // plot summary
          this->plot_summary();
 
          // return status
          return this->_status;
        }
 
        virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
        {
          // save rhs vector as initial defect
          this->_vec_rj_hat.at(0).copy(vec_def);
 
          // format solution vector
          vec_cor.format();
 
          // apply solver
          this->_status = _apply_intern(vec_cor, vec_def);
 
          // plot summary
          this->plot_summary();
 
          // return status
          return this->_status;
        }
 
 
        virtual void set_precon_variant(BiCGStabLPreconVariant precon_variant)
        {
          XASSERT(precon_variant == BiCGStabLPreconVariant::left || precon_variant == BiCGStabLPreconVariant::right);
          _precon_variant = precon_variant;
        }
 
      protected:
        Status _apply_intern(VectorType& vec_sol, const VectorType& vec_rhs)
        {
          IterationStats pre_iter(*this);
          Statistics::add_solver_expression(std::make_shared<ExpressionStartSolve>(this->name()));
 
 
          int l (_l);
 
          //choose  vec_r0_tilde as vec_r0_hat (vec_rj_hat[0])
          _vec_r0_tilde.copy(_vec_rj_hat.at(0));
          _vec_uj_hat.at(0).format();
 
          std::vector<DataType>  tau               ( Index( l*(l-1) ) , DataType(0));             //vector for tau_ij (j from 0 to l-1 , i from 0 to j-1, so LESS THAN HALF of the vector gets used)
                                                                                                  //tau_ij = tau[j*(l-1) + i] (column major ordering)
 
          std::vector<DataType>  sigma             ( Index(l) , DataType(0) );                    //0-offset in implementation; 1-offset in pseudo-code
          std::vector<DataType>  gamma_j           ( Index(l) , DataType(0) );                    //0-offset in implementation; 1-offset in pseudo-code
          std::vector<DataType>  gamma_prime       ( Index(l) , DataType(0) );                    //0-offset in implementation; 1-offset in pseudo-code
          std::vector<DataType>  gamma_prime_prime ( Index(l) , DataType(0) );                    //0-offset in implementation; 1-offset in pseudo-code
 
          DataType rho_0(1), alpha(0), omega(1), rho_1(0), beta(0), gamma(0);
 
 
          const MatrixType& mat_sys(this->_system_matrix);
          const FilterType& fil_sys(this->_system_filter);
          Status status(Status::progress);
 
          // Compute initial defect
          status = this->_set_initial_defect(_vec_rj_hat.at(0), vec_sol);
 
          if (_precon_variant ==  BiCGStabLPreconVariant::left)
          {
            _vec_pc.copy(_vec_rj_hat.at(0));
            if(!this->_apply_precond(_vec_rj_hat.at(0), _vec_pc, fil_sys))
            {
              Statistics::add_solver_expression(
                std::make_shared<ExpressionEndSolve>(this->name(), Status::aborted, this->get_num_iter()));
              return Status::aborted;
            }
          }
 
 
          pre_iter.destroy();
          while(status == Status::progress)
          {
            IterationStats stat(*this);
            rho_0 *= (-omega);
 
 
            for( int j(0); j < l ; j++)
            {
              rho_1 = _vec_rj_hat.at( Index(j) ).dot(_vec_r0_tilde);
              beta  = alpha * rho_1 / rho_0;
              rho_0 = rho_1;
              for ( int i(0) ; i <= j ; i++)
              {
                _vec_uj_hat.at(Index(i)).scale(_vec_uj_hat.at(Index(i)), -beta);
                _vec_uj_hat.at(Index(i)).axpy(_vec_rj_hat.at( Index(i)));
              }
 
 
 
              if (_precon_variant ==  BiCGStabLPreconVariant::left)
              {
                mat_sys.apply(_vec_pc, _vec_uj_hat.at(Index(j)) );
                fil_sys.filter_def(_vec_pc);
                if(!this->_apply_precond(_vec_uj_hat.at( Index(j+1) ), _vec_pc, fil_sys))
                {
                  Statistics::add_solver_expression(
                    std::make_shared<ExpressionEndSolve>(this->name(), Status::aborted, this->get_num_iter()));
                  return Status::aborted;
                }
              }
              else
              {
                if(!this->_apply_precond(_vec_pc, _vec_uj_hat.at( Index(j) ), fil_sys))
                {
                  Statistics::add_solver_expression(
                    std::make_shared<ExpressionEndSolve>(this->name(), Status::aborted, this->get_num_iter()));
                  return Status::aborted;
                }
                mat_sys.apply(_vec_uj_hat.at( Index(j+1) ), _vec_pc);
                fil_sys.filter_def(_vec_uj_hat.at( Index(j+1) ));
              }
 
              gamma = _vec_uj_hat.at( Index(j+1) ).dot(_vec_r0_tilde);
              alpha = rho_0 / gamma;
 
              for ( int i(0) ; i <= j ; i++)
              {
                _vec_rj_hat.at(Index(i)).axpy(_vec_uj_hat.at( Index(i+1)), -alpha);
              }
 
              if (_precon_variant ==  BiCGStabLPreconVariant::left)
              {
 
                mat_sys.apply( _vec_pc, _vec_rj_hat.at(Index(j)) );
                fil_sys.filter_def(_vec_pc);
                if(!this->_apply_precond(_vec_rj_hat.at( Index(j+1) ) , _vec_pc, fil_sys))
                {
                  Statistics::add_solver_expression(
                    std::make_shared<ExpressionEndSolve>(this->name(), Status::aborted, this->get_num_iter()));
                  return Status::aborted;
                }
              }
              else
              {
                if(!this->_apply_precond(_vec_pc, _vec_rj_hat.at(Index(j)), fil_sys))
                {
                  Statistics::add_solver_expression(
                    std::make_shared<ExpressionEndSolve>(this->name(), Status::aborted, this->get_num_iter()));
                  return Status::aborted;
                }
                mat_sys.apply(_vec_rj_hat.at(Index(j+1)), _vec_pc);
                fil_sys.filter_def( _vec_rj_hat.at(Index(j+1)) );
              }
 
              vec_sol.axpy(_vec_uj_hat.at(0), alpha);
            }
            for ( int j(0) ; j < l ; j++)       //changed counter in comparison to pseudo code
            {
 
              for ( int i(0); i<j ; i++)        //changed counter
              {
                tau.at( Index(j*(l-1) + i) ) = ( _vec_rj_hat.at( Index(i+1) ).dot( _vec_rj_hat.at( Index(j+1) ) ) ) / sigma.at(Index(i));
                _vec_rj_hat.at(Index(j+1)).axpy(_vec_rj_hat.at(Index(i+1)), -tau.at( Index(j*(l-1) + i)));
              }
              sigma.at(Index(j)) = _vec_rj_hat.at( Index(j+1) ).dot( _vec_rj_hat.at(Index(j+1)) );
              gamma_prime.at(Index(j)) = ( _vec_rj_hat.at(0).dot( _vec_rj_hat.at(Index(j+1)) ) ) / sigma.at(Index(j));
            }
            gamma_j.at(Index(l-1)) = gamma_prime.at(Index(l-1));
            omega = gamma_j.at(Index(l-1));
 
            for ( int j(l-2) ; j >= 0 ; j--)  //changed counter
            {
              gamma_j.at(Index(j)) = gamma_prime.at(Index(j));
              for ( int i(j+1) ; i < l ; i++)
              {
                gamma_j.at(Index(j)) -= tau.at( Index(i*(l-1) + j) ) * gamma_j.at(Index(i));
              }
            }
 
            for ( int j(0); j < l-1 ; j++)                                      //changed counter
            {
              gamma_prime_prime.at( Index(j) ) = gamma_j.at( Index(j+1) );
              for ( int i(j+1) ; i < l-1 ; i++)                                  //changed counter
              {
                gamma_prime_prime.at(Index(j)) += tau.at(Index(i*(l-1) + j)) * gamma_j.at(Index(i+1));
              }
            }
                      vec_sol.axpy(_vec_rj_hat.at(0),         gamma_j.at(0));
            _vec_rj_hat.at(0).axpy(_vec_rj_hat.at(Index(l)), -gamma_prime.at(Index(l-1)));
            _vec_uj_hat.at(0).axpy(_vec_uj_hat.at(Index(l)), -gamma_j.at(Index(l-1)));
 
            for ( int j(1) ; j<l; j++)
            {
              _vec_uj_hat.at(0).axpy(_vec_uj_hat.at(Index(j)), -gamma_j.at(Index(j-1)));
                        vec_sol.axpy(_vec_rj_hat.at(Index(j)), gamma_prime_prime.at(Index(j-1)));
              _vec_rj_hat.at(0).axpy(_vec_rj_hat.at(Index(j)), -gamma_prime.at(Index(j-1)));
            }
            // Compute defect norm
            if (_precon_variant ==  BiCGStabLPreconVariant::left)
            {
              mat_sys.apply(_vec_pc, vec_sol);
              fil_sys.filter_def(_vec_pc);
              _vec_pc.axpy( vec_rhs, -1);
              fil_sys.filter_def(_vec_pc);
              status = this->_set_new_defect(_vec_pc, vec_sol);
            }
            else
            {
              status = this->_set_new_defect(_vec_rj_hat.at(0), vec_sol);
            }
 
            if(status != Status::progress)
            {
              stat.destroy();
              Statistics::add_solver_expression(
                std::make_shared<ExpressionEndSolve>(this->name(), status, this->get_num_iter()));
              if (_precon_variant == BiCGStabLPreconVariant::right)
              {
                if(!this->_apply_precond(_vec_pc, vec_sol, fil_sys))
                {
                  Statistics::add_solver_expression(
                    std::make_shared<ExpressionEndSolve>(this->name(), Status::aborted, this->get_num_iter()));
                  return Status::aborted;
                }
                vec_sol.copy(_vec_pc);
              }
              return status;
            }
 
          }
          // we should never reach this point...
          Statistics::add_solver_expression(
            std::make_shared<ExpressionEndSolve>(this->name(), Status::undefined, this->get_num_iter()));
          return Status::undefined;
        }
    };
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<BiCGStabL<Matrix_, Filter_>> new_bicgstabl(
      const Matrix_& matrix, const Filter_& filter, int l = 2,
      std::shared_ptr<SolverBase<typename Matrix_::VectorTypeL>> precond = nullptr,
      BiCGStabLPreconVariant precon_variant = BiCGStabLPreconVariant::left)
    {
      return std::make_shared<BiCGStabL<Matrix_, Filter_>>(matrix, filter, l, precond, precon_variant);
    }
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<BiCGStabL<Matrix_, Filter_>> new_bicgstabl(
      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<BiCGStabL<Matrix_, Filter_>>(section_name, section, matrix, filter, precond);
    }
 
  } // namespace Solver
} // namespace FEAT