UmfpackRunner.hpp 3.29 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#pragma once

#ifdef HAVE_UMFPACK

#include <iostream>

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

#include <dune/amdis/linear_algebra/RunnerInterface.hpp>
#include <dune/amdis/linear_algebra/SolverInfo.hpp>

#include <dune/amdis/linear_algebra/mtl/Copy.hpp>

namespace AMDiS 
{
  using namespace mtl::mat::umfpack;
  
  /**
   * \ingroup Solver
   * \class AMDiS::UmfPackSolver
   * \brief \implements LinearSolver
   * Wrapper for the external UMFPACK solver:
   *   http://www.cise.ufl.edu/research/sparse/umfpack/
   *
   * This is a direct solver for large sparse matrices.
   */
  template <class Matrix, class Vector>
  class UmfpackRunner : public RunnerInterface<Matrix, Vector, Vector>
  {
    using Super = RunnerInterface<Matrix, Vector>;
    using PreconBase = typename Super::PreconBase;
    
    // extract the matrix type used by UMFPACK
    using FullMatrix = typename mtl::mat::umfpack::matrix_copy<
        typename Matrix::BaseMatrix, 
        typename Matrix::value_type, 
        typename Matrix::BaseMatrix::parameters::orientation
      >::matrix_type;
    using SolverType = mtl::mat::umfpack::solver<FullMatrix>;
    
  public:
    /// Constructor. Reads UMFPACK parameters from initfile
    UmfpackRunner(std::string prefix)
      : solver(NULL)
      , store_symbolic(0)
      , symmetric_strategy(0)
      , alloc_init(0.7)
    {
      Parameters::get(prefix + "->store symbolic", store_symbolic); // ?
      Parameters::get(prefix + "->symmetric strategy", symmetric_strategy);
      Parameters::get(prefix + "->alloc init", alloc_init);
    }

    /// Destructor
    ~UmfpackRunner() = default;

    /// Implementation of \ref RunnerInterface::init()
    virtual void init(Matrix const& A) override
    {
      AMDIS_MSG("copy matrix...");
      copy(A, fullMatrix);
      AMDIS_MSG("size(fullMatrix) = " << num_rows(fullMatrix) << " x " << num_cols(fullMatrix));
      
      {
        mtl::io::matrix_market_ostream out("matrix.mtx");
        out << fullMatrix;
      }
      
      AMDIS_MSG("create factorization...");
      try {
	solver.reset(new SolverType(fullMatrix, symmetric_strategy, alloc_init));
      } catch (mtl::mat::umfpack::error& e) {
	AMDIS_ERROR_EXIT("UMFPACK_ERROR(factorize, " << e.code << ") = " << e.what());
      }
    }

    /// Implementation of \ref RunnerBase::solve()
    virtual int solve(Matrix const& A, Vector& x, Vector const& b, 
                      SolverInfo& solverInfo) override
    { 
      AMDIS_FUNCNAME("Umfpack_Runner::solve()");
      AMDIS_TEST_EXIT(solver, "The umfpack solver was not initialized\n");
      
      AMDIS_MSG("solve system...");
      int code = 0;
      try {
	code = (*solver)(x, b);
      } catch (mtl::mat::umfpack::error& e) {
	AMDIS_ERROR_EXIT("UMFPACK_ERROR(solve, " << e.code << ") = " << e.what());
      }
      
      auto r = Vector(b);
      if (two_norm(x) != 0)
        r -= A * x;
      
      solverInfo.setAbsResidual(two_norm(r));
      solverInfo.setError(code);
      
      return code;
    }

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

  private:
    FullMatrix fullMatrix;
    std::shared_ptr<SolverType> solver;

    int store_symbolic;
    int symmetric_strategy;
    double alloc_init;
  };
  
}

#endif // HAVE_UMFPACK