UmfpackRunner.hpp 4.14 KB
Newer Older
1
2
#pragma once

3
4
5
6
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

7
8
#ifdef HAVE_UMFPACK

9
#include <algorithm>
10
#include <iostream>
11
#include <string>
12
13
14
15

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

16
17
#include <amdis/linear_algebra/RunnerInterface.hpp>
#include <amdis/linear_algebra/SolverInfo.hpp>
18

19
#include <amdis/linear_algebra/mtl/Copy.hpp>
20

21
22
namespace AMDiS
{
23
24
  /**
   * \ingroup Solver
25
26
   * \class AMDiS::UmfpackRunner
   * \brief \implements RunnerInterface
27
   * Wrapper for the external
28
   *   [UMFPACK solver](http://faculty.cse.tamu.edu/davis/suitesparse.html)
29
30
31
   *
   * This is a direct solver for large sparse matrices.
   */
32
33
  template <class Matrix, class Vector>
  class UmfpackRunner;
34

35
  template <class Matrix, class Vector>
36
37
  class UmfpackRunnerBase
      : public RunnerInterface<Matrix, Vector, Vector>
38
  {
39
  protected:
40
41
    using Super = RunnerInterface<Matrix, Vector>;
    using PreconBase = typename Super::PreconBase;
42

43
    using FullMatrix = BaseMatrix_t<Matrix>;
44
    using SolverType = mtl::mat::umfpack::solver<FullMatrix>;
45

46
47
  public:
    /// Constructor. Reads UMFPACK parameters from initfile
48
    UmfpackRunnerBase(std::string prefix)
49
50
51
52
53
54
55
    {
      Parameters::get(prefix + "->store symbolic", store_symbolic); // ?
      Parameters::get(prefix + "->symmetric strategy", symmetric_strategy);
      Parameters::get(prefix + "->alloc init", alloc_init);
    }

    /// Destructor
56
    ~UmfpackRunnerBase() = default;
57
58

    /// Implementation of \ref RunnerBase::solve()
59
    virtual int solve(Matrix const& A, Vector& x, Vector const& b,
60
                      SolverInfo& solverInfo) override
61
    {
62
      AMDIS_FUNCNAME("Umfpack_Runner::solve()");
63
64
      test_exit(bool(solver), "The umfpack solver was not initialized\n");

65
66
      int code = 0;
      try {
67
	      code = (*solver)(x, b);
68
      } catch (mtl::mat::umfpack::error& e) {
69
	      error_exit("UMFPACK_ERROR(solve, ", e.code, ") = ", e.what());
70
      }
71

72
73
74
      auto r = Vector(b);
      if (two_norm(x) != 0)
        r -= A * x;
75

76
77
      solverInfo.setAbsResidual(two_norm(r));
      solverInfo.setError(code);
78

79
80
81
82
83
84
      return code;
    }

    /// Implementation of \ref RunnerInterface::adjoint_solve()
    virtual void exit() override {}

85
  protected:
86
    std::shared_ptr<SolverType> solver;
87

88
89
90
    int store_symbolic = 0;
    int symmetric_strategy = 0;
    double alloc_init = 0.7;
91
  };
92

93
  // specialization for block-matrices
94
  template <class SubMatrix, std::size_t _N, std::size_t _M, class Vector>
95
  class UmfpackRunner<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>
96
      : public UmfpackRunnerBase<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>
97
98
99
  {
    using Matrix = BlockMTLMatrix<SubMatrix, _N, _M>;
    using Super = UmfpackRunnerBase<Matrix, Vector>;
100

101
    using SolverType = typename Super::SolverType;
102

103
104
105
106
107
108
109
110
111
  public:
    UmfpackRunner(std::string prefix)
      : Super(prefix)
    {}

    /// Implementation of \ref RunnerInterface::init()
    virtual void init(Matrix const& A) override
    {
      copy(A, fullMatrix);
112

113
      try {
114
115
        Super::solver.reset(new SolverType(fullMatrix, Super::symmetric_strategy, Super::alloc_init));
      } catch (mtl::mat::umfpack::error const& e) {
116
        error_exit("UMFPACK_ERROR(factorize, ", e.code, ") = ", e.what());
117
118
      }
    }
119

120
121
  private:
    typename Super::FullMatrix fullMatrix;
122
  };
123

124
125
  // specialization for full-matrices
  template <class T, class Param, class Vector>
126
  class UmfpackRunner<mtl::compressed2D<T, Param>, Vector>
127
      : public UmfpackRunnerBase<mtl::compressed2D<T, Param>, Vector>
128
129
130
  {
    using FullMatrix = mtl::compressed2D<T, Param>;
    using Super = UmfpackRunnerBase<FullMatrix, Vector>;
131

132
    using SolverType = typename Super::SolverType;
133

134
135
136
137
138
139
140
141
142
  public:
    UmfpackRunner(std::string prefix)
      : Super(prefix)
    {}

    /// Implementation of \ref RunnerInterface::init()
    virtual void init(FullMatrix const& fullMatrix) override
    {
      try {
143
144
        Super::solver.reset(new SolverType(fullMatrix, Super::symmetric_strategy, Super::alloc_init));
      } catch (mtl::mat::umfpack::error const& e) {
145
        error_exit("UMFPACK_ERROR(factorize, ", e.code, ") = ", e.what());
146
147
148
      }
    }
  };
149
150
151
}

#endif // HAVE_UMFPACK