UmfpackRunner.hpp 3.99 KB
Newer Older
1
2
3
4
#pragma once

#ifdef HAVE_UMFPACK

5
6
#include <algorithm>
#include <string>
7
8
9
10

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

11
#include <amdis/linearalgebra/mtl/Traits.hpp>
12
13
#include <amdis/linearalgebra/RunnerInterface.hpp>
#include <amdis/linearalgebra/SolverInfo.hpp>
14
#include <amdis/Output.hpp>
15

16
17
namespace AMDiS
{
18
  /**
19
20
21
22
23
24
25
26
  * \ingroup Solver
  * \class AMDiS::UmfpackRunner
  * \brief \implements RunnerInterface
  * Wrapper for the external
  *   [UMFPACK solver](http://faculty.cse.tamu.edu/davis/suitesparse.html)
  *
  * This is a direct solver for large sparse matrices.
  */
Praetorius, Simon's avatar
Praetorius, Simon committed
27
  template <class Mat, class Vec>
28
  class UmfpackRunner
Praetorius, Simon's avatar
Praetorius, Simon committed
29
      : public RunnerInterface<Mat,Vec>
30
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
31
32
33
    using M = typename Mat::BaseMatrix;
    using X = typename Vec::BaseVector;
    using Y = typename Vec::BaseVector;
34

Praetorius, Simon's avatar
Praetorius, Simon committed
35
    using SolverType = mtl::mat::umfpack::solver<typename Mat::BaseMatrix>;
36

37
38
  public:
    /// Constructor. Reads UMFPACK parameters from initfile
39
    UmfpackRunner(std::string const& prefix)
40
    {
41
      Parameters::get(prefix + "->store symbolic", storeSymbolic_);
42
43
      Parameters::get(prefix + "->symmetric strategy", symmetricStrategy_);
      Parameters::get(prefix + "->alloc init", allocInit_);
44
45
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
46
    void init(M const& A) override
47
48
    {
      try {
49
50
51
        if (bool(solver_) && storeSymbolic_)
          solver_->update_numeric();
        else
Praetorius, Simon's avatar
Praetorius, Simon committed
52
          solver_ = std::make_shared<SolverType>(A, symmetricStrategy_, allocInit_);
53
      } catch (mtl::mat::umfpack::error const& e) {
54
        umfpack_error_exit("factorize", e.code);
55
56
57
58
59
      }
    }

    void exit() override {}

Praetorius, Simon's avatar
Praetorius, Simon committed
60
61
    /// Implements \ref RunnerInterface::solveImpl()
    int solve(M const& A, X& x, Y const& b, SolverInfo& solverInfo) override
62
    {
63
      test_exit(bool(solver_), "UmfpackRunner must be initialized.");
64
65
      int code = 0;
      try {
66
        code = (*solver_)(x, b);
67
      } catch (mtl::mat::umfpack::error& e) {
68
        umfpack_error_exit("solve", e.code);
69
      }
70

71
      solverInfo.setAbsResidual(two_norm(b - A*x));
72
73
74
      return code;
    }

75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
    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]);
    }

99
  protected:
100
    std::shared_ptr<SolverType> solver_;
101

102
    bool storeSymbolic_ = false;
103
104
    int symmetricStrategy_ = 0;
    double allocInit_ = 0.7;
105
  };
106
107

} // end namespace AMDiS
108
109

#endif // HAVE_UMFPACK