UmfpackRunner.hpp

#pragma once

#ifdef HAVE_UMFPACK

#include <algorithm>
#include <string>

#include <dune/istl/solver.hh>

#include <boost/numeric/mtl/operation/two_norm.hpp>
#include <boost/numeric/mtl/interface/umfpack_solve.hpp>

#include <amdis/linearalgebra/mtl/Traits.hpp>
#include <amdis/linearalgebra/LinearSolverInterface.hpp>
#include <amdis/Output.hpp>

namespace AMDiS
{
  template <class M, class X, class Y>
  class UmfpackRunner
      : public LinearSolverInterface<M,X,Y>
  {
    using Self = UmfpackRunner;
    using SolverType = mtl::mat::umfpack::solver<M>;

  public:
    struct Creator final : CreatorInterfaceName<LinearSolverInterface<M,X,Y>>
    {
      std::unique_ptr<LinearSolverInterface<M,X,Y>>
      createWithString(std::string prefix) final
      {
        return std::make_unique<Self>(prefix);
      }
    };

  public:
    UmfpackRunner(std::string const& prefix)
    {
      Parameters::get(prefix + "->store symbolic", storeSymbolic_);
      Parameters::get(prefix + "->symmetric strategy", symmetricStrategy_);
      Parameters::get(prefix + "->alloc init", allocInit_);
    }

    void init(M const& A) override
    {
      try {
        if (bool(solver_) && storeSymbolic_)
          solver_->update_numeric();
        else
          solver_ = std::make_shared<SolverType>(A, symmetricStrategy_, allocInit_);
      } catch (mtl::mat::umfpack::error const& e) {
        umfpack_error_exit("factorize", e.code);
      }
    }

    void finish() override {}

    void apply(X& x, Y const& b, Dune::InverseOperatorResult& stat) override
    {
      test_exit(bool(solver_), "UmfpackRunner must be initialized.");
      Dune::Timer t;
      int code = 0;
      try {
        code = (*solver_)(x, b);
      } catch (mtl::mat::umfpack::error& e) {
        umfpack_error_exit("solve", e.code);
      }

      stat.iterations = 1;
      stat.converged = (code == UMFPACK_OK);
      stat.elapsed = t.elapsed();
    }

  protected:
    static void umfpack_error_exit(std::string solutionPhase, int code)
    {
      static std::map<int, std::string> error_message = {
        {UMFPACK_OK, "UMFPACK was successful"},
        {UMFPACK_WARNING_singular_matrix, "Matrix is singular. There are exact zeros on the diagonal."},
        {UMFPACK_WARNING_determinant_underflow, "The determinant is nonzero, but smaller in magnitude than the smallest positive floating-point number."},
        {UMFPACK_WARNING_determinant_overflow, "The determinant is larger in magnitude than the largest positive floating-point number (IEEE Inf)."},
        {UMFPACK_ERROR_out_of_memory, "Not enough memory. The ANSI C malloc or realloc routine failed."},
        {UMFPACK_ERROR_invalid_Numeric_object, "Invalid Numeric object."},
        {UMFPACK_ERROR_invalid_Symbolic_object, "Invalid Symbolic object."},
        {UMFPACK_ERROR_argument_missing, "Some required arguments are missing."},
        {UMFPACK_ERROR_n_nonpositive, "The number of rows or columns of the matrix must be greater than zero."},
        {UMFPACK_ERROR_invalid_matrix, "The matrix is invalid."},
        {UMFPACK_ERROR_different_pattern, "The pattern of the matrix has changed between the symbolic and numeric factorization."},
        {UMFPACK_ERROR_invalid_system, "The sys argument provided to one of the solve routines is invalid."},
        {UMFPACK_ERROR_invalid_permutation, "The permutation vector provided as input is invalid."},
        {UMFPACK_ERROR_file_IO, "Error in file IO"},
        {UMFPACK_ERROR_ordering_failed, "The ordering method failed."},
        {UMFPACK_ERROR_internal_error, "An internal error has occurred, of unknown cause."}
      };

      error_exit("UMFPACK_ERROR({}, {}) = {}", solutionPhase, code, error_message[code]);
    }

  protected:
    std::shared_ptr<SolverType> solver_;

    bool storeSymbolic_ = false;
    int symmetricStrategy_ = 0;
    double allocInit_ = 0.7;
  };

} // end namespace AMDiS

#endif // HAVE_UMFPACK