Solvers.hpp

#pragma once

#include <memory>

#include <Eigen/IterativeLinearSolvers>
#include <Eigen/MetisSupport>
#include <Eigen/SparseLU>
#include <Eigen/SuperLUSupport>
#include <Eigen/UmfPackSupport>
#include <unsupported/Eigen/IterativeSolvers>

#include <amdis/CreatorMap.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/Output.hpp>

#include <amdis/linearalgebra/LinearSolverInterface.hpp>
#include <amdis/linearalgebra/eigen/DirectRunner.hpp>
#include <amdis/linearalgebra/eigen/IterativeRunner.hpp>
#include <amdis/linearalgebra/eigen/Traits.hpp>

namespace AMDiS
{
  template <class M, class X, class Y, template <class> class IterativeSolver>
  struct IterativeSolverCreator
      : public CreatorInterfaceName<LinearSolverInterface<M,X,Y>>
  {
    using SolverBase = LinearSolverInterface<M,X,Y>;
    using T = typename M::Scalar;

    template <class Precon>
    using Solver = IterativeRunner<M,X,Y,IterativeSolver<Precon>>;

    template <class Ordering>
    using IncompleteCholesky =
      Solver<Eigen::IncompleteCholesky<T, Eigen::Lower|Eigen::Upper, Ordering>>;

    std::unique_ptr<SolverBase> createWithString(std::string prefix) override
    {
      // get creator string for preconditioner
      std::string precon = "no";
      Parameters::get(prefix + "->precon->name", precon);

      if (precon == "diag" || precon == "jacobi") {
        return std::make_unique<Solver<Eigen::DiagonalPreconditioner<T>>>(prefix);
      }
      else if (precon == "ilu") {
        return std::make_unique<Solver<Eigen::IncompleteLUT<T>>>(prefix);
      }
      else if (precon == "ic") {
        std::string ordering = "amd";
        Parameters::get(prefix + "->precon->ordering", ordering);

        if (ordering == "amd") {
          using AMD = Eigen::AMDOrdering<typename M::StorageIndex>;
          return std::make_unique<IncompleteCholesky<AMD>>(prefix);
        }
        else {
          using NATURAL = Eigen::NaturalOrdering<typename M::StorageIndex>;
          return std::make_unique<IncompleteCholesky<NATURAL>>(prefix);
        }
      }
      else {
        return std::make_unique<Solver<Eigen::IdentityPreconditioner>>(prefix);
      }
    }
  };


  template <class M, class X, class Y>
  class DefaultCreators< LinearSolverInterface<M,X,Y> >
  {
    template <template <class> class IterativeSolver>
    using SolverCreator
      = IterativeSolverCreator<M,X,Y,IterativeSolver>;

    template <template <class> class DirectSolver>
    using DirectSolverCreator
      = typename DirectRunner<M,X,Y,DirectSolver>::Creator;

    template <class Precon>
    using CG = Eigen::ConjugateGradient<M, Eigen::Lower|Eigen::Upper, Precon>;

    template <class Precon>
    using BCGS = Eigen::BiCGSTAB<M, Precon>;

    template <class Precon>
    using MINRES = Eigen::MINRES<M, Eigen::Lower|Eigen::Upper, Precon>;

    template <class Precon>
    using GMRES = Eigen::GMRES<M, Precon>;

    template <class Precon>
    using DGMRES = Eigen::DGMRES<M, Precon>;
    // @}

    using SolverBase = LinearSolverInterface<M,X,Y>;
    using Map = CreatorMap<SolverBase>;

  public:
    static void init()
    {
      // conjugate gradient solver
      auto cg = new SolverCreator<CG>;
      Map::addCreator("cg", cg);

      // bi-conjugate gradient stabilized solver
      auto bicgstab = new SolverCreator<BCGS>;
      Map::addCreator("bicgstab", bicgstab);
      Map::addCreator("bcgs", bicgstab);

      // Minimal Residual Algorithm (for symmetric matrices)
      auto minres = new SolverCreator<MINRES>;
      Map::addCreator("minres", minres);

      // Generalized Minimal Residual Algorithm
      auto gmres = new SolverCreator<GMRES>;
      Map::addCreator("gmres", gmres);

      // Restarted GMRES with deflation.
      auto dgmres = new SolverCreator<DGMRES>;
      Map::addCreator("dgmres", dgmres);

      // default iterative solver
      Map::addCreator("default", gmres);

      init_direct(std::is_same<typename Dune::FieldTraits<typename M::Scalar>::real_type, double>{});
    }

    static void init_direct(std::false_type) {}
    static void init_direct(std::true_type)
    {
#if HAVE_SUITESPARSE_UMFPACK
      // sparse LU factorization and solver based on UmfPack
      auto umfpack = new DirectSolverCreator<Eigen::UmfPackLU>;
      Map::addCreator("umfpack", umfpack);
#endif

#if HAVE_SUPERLU
      // sparse direct LU factorization and solver based on the SuperLU library
      auto superlu = new DirectSolverCreator<Eigen::SuperLU>;
      Map::addCreator("superlu", superlu);
#endif

      // default direct solvers
#if HAVE_SUITESPARSE_UMFPACK
      Map::addCreator("direct", umfpack);
#elif HAVE_SUPERLU
      Map::addCreator("direct", superlu);
#endif
    }
  };

} // end namespace AMDiS