Solvers.hpp

#pragma once
 
#include <string>
 
// MTL4 includes
#include <boost/numeric/itl/itl.hpp>
 
// more solvers defined in AMDiS
#include <amdis/linearalgebra/mtl/itl/minres.hpp>
#include <amdis/linearalgebra/mtl/itl/gcr.hpp>
#include <amdis/linearalgebra/mtl/itl/fgmres.hpp>
#include <amdis/linearalgebra/mtl/itl/fgmres_householder.hpp>
#include <amdis/linearalgebra/mtl/itl/gmres2.hpp>
#include <amdis/linearalgebra/mtl/itl/gmres_householder.hpp>
#include <amdis/linearalgebra/mtl/itl/preonly.hpp>
 
#include <amdis/CreatorMap.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/linearalgebra/mtl/KrylovRunner.hpp>
#include <amdis/linearalgebra/mtl/Traits.hpp>
#include <amdis/linearalgebra/mtl/UmfpackRunner.hpp>
 
namespace AMDiS
{
  class cg_solver_type
  {
  public:
    explicit cg_solver_type(std::string const& /*prefix*/) {}
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      return itl::cg(A, x, b, precon, iter);
    }
  };
 
 
  class cgs_solver_type
  {
  public:
    explicit cgs_solver_type(std::string const& /*prefix*/) {}
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      return itl::cgs(A, x, b, precon, iter);
    }
  };
 
 
  class bicg_solver_type
  {
  public:
    explicit bicg_solver_type(std::string const& /*prefix*/) {}
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      return itl::bicg(A, x, b, precon, iter);
    }
  };
 
 
  class bicgstab_type
  {
  public:
    explicit bicgstab_type(std::string const& /*prefix*/) {}
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      return itl::bicgstab(A, x, b, precon, iter);
    }
  };
 
 
  class bicgstab2_type
  {
  public:
    explicit bicgstab2_type(std::string const& /*prefix*/) {}
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      return itl::bicgstab_2(A, x, b, precon, iter);
    }
  };
 
 
  class qmr_solver_type
  {
  public:
    explicit qmr_solver_type(std::string const& /*prefix*/) {}
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      itl::pc::identity<LinOp> id(A);
      return itl::qmr(A, x, b, precon, id, iter);
    }
  };
 
 
  class tfqmr_solver_type
  {
  public:
    explicit tfqmr_solver_type(std::string const& /*prefix*/) {}
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      itl::pc::identity<LinOp> id(A);
      return itl::tfqmr(A, x, b, precon, id, iter);
    }
  };
 
 
  class bicgstab_ell_type
  {
    int ell;
  public:
    explicit bicgstab_ell_type(std::string const& prefix) : ell(3)
    {
      Parameters::get(prefix + "->ell", ell);
    }
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      itl::pc::identity<LinOp> id(A);
      return itl::bicgstab_ell(A, x, b, precon, id, iter, ell);
    }
  };
 
 
  class gmres_type
  {
 
  public:
    explicit gmres_type(std::string const& prefix) : restart(30), ortho(GRAM_SCHMIDT)
    {
      Parameters::get(prefix + "->restart", restart);
      Parameters::get(prefix + "->orthogonalization", ortho);
    }
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      itl::pc::identity<LinOp> id(A);
      switch (ortho)
      {
      default:
      case GRAM_SCHMIDT:
        return itl::gmres2(A, x, b, id, precon, iter, restart);
        break;
      case HOUSEHOLDER:
        return itl::gmres_householder(A, x, b, precon, iter, restart);
        break;
      }
    }
 
  private:
    static constexpr int GRAM_SCHMIDT = 1;
    static constexpr int HOUSEHOLDER = 2;
 
    int restart;
    int ortho;
  };
 
 
  class idr_s_type
  {
    int s;
  public:
    explicit idr_s_type(std::string const& prefix) : s(30)
    {
      Parameters::get(prefix + "->s", s);
    }
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const&, I& iter)
    {
      itl::pc::identity<LinOp> id(A);
      return itl::idr_s(A, x, b, id, id, iter, s);
    }
  };
 
 
  class minres_solver_type
  {
  public:
    explicit minres_solver_type(std::string const& /*prefix*/) {}
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      return itl::minres(A, x, b, precon, iter);
    }
  };
 
 
  class gcr_type
  {
    int restart;
 
  public:
    explicit gcr_type(std::string const& prefix) : restart(30)
    {
      Parameters::get(prefix + "->restart", restart);
    }
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      itl::pc::identity<LinOp> id(A);
      return itl::gcr(A, x, b, id, precon, iter, restart);
    }
  };
 
 
  class fgmres_type
  {
  public:
    explicit fgmres_type(std::string const& prefix)
      : restart(30), orthogonalization(GRAM_SCHMIDT)
    {
      Parameters::get(prefix + "->restart", restart);
      Parameters::get(prefix + "->orthogonalization", orthogonalization);
    }
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      itl::pc::identity<LinOp> id(A);
      switch (orthogonalization)
      {
      default:
      case GRAM_SCHMIDT:
        return itl::fgmres(A, x, b, id, precon, iter, restart);
        break;
      case HOUSEHOLDER:
        return itl::fgmres_householder(A, x, b, precon, iter, restart);
        break;
      }
    }
 
  private:
    static constexpr int GRAM_SCHMIDT = 1;
    static constexpr int HOUSEHOLDER = 2;
 
    int restart;
    int orthogonalization;
  };
 
 
  class preonly_type
  {
  public:
    explicit preonly_type(std::string const& /*prefix*/) {}
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
    {
      return itl::preonly(A, x, b, precon, iter);
    }
  };
 
 
  // ===========================================================================
 
 
  template <class M, class X, class Y>
  class DefaultCreators< LinearSolverInterface<M,X,Y> >
  {
    using SolverBase = LinearSolverInterface<M,X,Y>;
 
    template <class ITLSolver>
    using SolverCreator = typename KrylovRunner<M,X,Y,ITLSolver>::Creator;
 
#if HAVE_SUITESPARSE_UMFPACK
    using UmfpackSolverCreator = typename UmfpackRunner<M,X,Y>::Creator;
#endif
 
    using Map = CreatorMap<SolverBase>;
 
  public:
    static void init()
    {
      auto cg = new SolverCreator<cg_solver_type>;
      Map::addCreator("cg", cg);
 
      auto cgs = new SolverCreator<cgs_solver_type>;
      Map::addCreator("cgs", cgs);
 
      auto bicg = new SolverCreator<bicg_solver_type>;
      Map::addCreator("bicg", bicg);
 
      auto bcgs = new SolverCreator<bicgstab_type>;
      Map::addCreator("bicgstab", bcgs);
      Map::addCreator("bcgs", bcgs);
 
      auto bcgs2 = new SolverCreator<bicgstab2_type>;
      Map::addCreator("bicgstab2", bcgs2);
      Map::addCreator("bcgs2", bcgs2);
 
      auto qmr = new SolverCreator<qmr_solver_type>;
      Map::addCreator("qmr", qmr);
 
      auto tfqmr = new SolverCreator<tfqmr_solver_type>;
      Map::addCreator("tfqmr", tfqmr);
 
      auto bcgsl = new SolverCreator<bicgstab_ell_type>;
      Map::addCreator("bicgstab_ell", bcgsl);
      Map::addCreator("bcgsl", bcgsl);
 
      auto gmres = new SolverCreator<gmres_type>;
      Map::addCreator("gmres", gmres);
 
      auto fgmres = new SolverCreator<fgmres_type>;
      Map::addCreator("fgmres", fgmres);
 
      auto minres = new SolverCreator<minres_solver_type>;
      Map::addCreator("minres", minres);
 
      auto idrs = new SolverCreator<idr_s_type>;
      Map::addCreator("idrs", idrs);
 
      auto gcr = new SolverCreator<gcr_type>;
      Map::addCreator("gcr", gcr);
 
      auto preonly = new SolverCreator<preonly_type>;
      Map::addCreator("preonly", preonly);
 
      // default iterative solver
      Map::addCreator("default", gmres);
 
      init_direct(std::is_same<typename Dune::FieldTraits<typename M::value_type>::real_type, double>{});
    }
 
    static void init_direct(std::false_type) {}
    static void init_direct(std::true_type)
    {
#if HAVE_SUITESPARSE_UMFPACK
      auto umfpack = new UmfpackSolverCreator;
      Map::addCreator("umfpack", umfpack);
 
      // default direct solvers
      Map::addCreator("direct", umfpack);
#endif
    }
  };
 
} // end namespace AMDiS