UmfPackSolver.h 4.15 KB
Newer Older
1
2
3
4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6
7
// ==                                                                        ==
// ============================================================================
8
9
10
11
12
13
14
15
16
17
18
19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


20
21
22
23
24
25

/** \file UmfPackSolver.h */

#ifndef AMDIS_UMFPACKSOLVER_H
#define AMDIS_UMFPACKSOLVER_H

26
#if defined HAVE_UMFPACK && defined MTL_HAS_UMFPACK
27

28
29
#include <iostream>
#include <boost/numeric/mtl/operation/two_norm.hpp>
30
#include <boost/numeric/mtl/interface/umfpack_solve.hpp>
31
#include "MTL4Solver.h"
32

33
34
35
36
37
38
39
40
41
42
43
namespace AMDiS {

  /**
   * \ingroup Solver
   * 
   * \brief
   * Wrapper for the external UMFPACK solver:
   *   http://www.cise.ufl.edu/research/sparse/umfpack/
   *
   * This is a direct solver for large sparse matrices.
   */
44
45
  template< typename Matrix >
  struct Umfpack_runner 
46
  {
47
48
49
50
51
52
    typedef Matrix matrix_type ;
    OEMSolver& oem;
    Umfpack_runner(OEMSolver* oem_):
      oem(*oem_),
      solver(NULL),
      store_symbolic(0),
53
54
      symmetric_strategy(0),
      alloc_init(0.7)
55
    {
56
57
      FUNCNAME("Umfpack_runner::Umfpack_runner()");

58
59
60
      TEST_EXIT_DBG(oem_ != NULL)("Need real OEMSolver\n"); 
      Parameters::get(oem.getName() + "->store symbolic", store_symbolic);
      Parameters::get(oem.getName() + "->symmetric strategy", symmetric_strategy);
61
      Parameters::get(oem.getName() + "->alloc init", alloc_init);
Thomas Witkowski's avatar
Thomas Witkowski committed
62
    }
63

64
65
66
    void init(const Matrix& A)
    {
      FUNCNAME("Umfpack_runner::init()")
67
68
69
70
71
72

      if (solver != NULL) {
	delete solver;
	solver = NULL;
      }

73
74
75
      solver = new mtl::matrix::umfpack::solver<matrix_type>(A, symmetric_strategy, alloc_init);
    }

76
77
78
    template< typename Vector> 
    int solve(const Matrix& A, Vector& x, Vector& b) 
    { 
79
      FUNCNAME("Umfpack_runner::solve()");
80
81
      TEST_EXIT(solver != NULL)("The umfpack solver was not initialized\n");
#if 0 
82
83
84
85
      if (store_symbolic)
	solver->update_numeric();
      else
	solver->update();
86
#endif
87

Thomas Witkowski's avatar
Thomas Witkowski committed
88
      int code = (*solver)(x, b);
89
      if (oem.getInfo() > 0) {
90
        mtl::dense_vector<typename matrix_type::value_type> r(b); 
91
        r -= A * x; 
92
93
94
        double residual = two_norm(r);
        oem.setResidual(residual);
        oem.setErrorCode(code);
95
        std::cout << "UmfPackSolver: ||b-Ax|| = " << residual << "\n";
96
	TEST_EXIT(residual <= oem.getTolerance())("Tolerance tol = %e could not be reached!\n Set tolerance by '->solver->tolerance:' \n", oem.getTolerance());	
97
      }
98
      return code;
99
    }
100

101
102
103
    void exit()
    {}

104
105
    ~Umfpack_runner() 
    {
106
      if (solver != NULL) {	
107
	delete solver;
108
109
	solver == NULL;
      }
110
111
    }

112
  private:
113
    mtl::matrix::umfpack::solver<matrix_type> *solver;
Thomas Witkowski's avatar
Thomas Witkowski committed
114

115
    int store_symbolic;
116

117
    int symmetric_strategy;
118
119

    double alloc_init;
120
  };
121
122
123
124

  using namespace MTLTypes;
  class UmfPackSolver : public MTL4Solver<MTLMatrix, MTLVector, Umfpack_runner< MTLMatrix > >
  {
125
  protected:
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

  public:
    /// Creator class used in the OEMSolverMap.
    class Creator : public OEMSolverCreator
    {
    public:
      virtual ~Creator() {}
      
      /// Returns a new UmfPackSolver object.
      OEMSolver* create() 
      { 
	return new UmfPackSolver(this->name); 
      }
    };

    /// Constructor
    UmfPackSolver(std::string name) 
      : MTL4Solver< MTLMatrix, MTLVector, Umfpack_runner< MTLMatrix > >(name)
144
145
    {}
    
Naumann, Andreas's avatar
Naumann, Andreas committed
146
147
148
    int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
		    SystemVector& x, 
		    SystemVector& b,
149
150
151
		    VectorialMapper& m,
		    bool createMatrixData,
		    bool storeMatrixData) 
Naumann, Andreas's avatar
Naumann, Andreas committed
152
    {
153
154
      return solve(A, x, b, m, createMatrixData, storeMatrixData);
    }       
155
  };
156
157
158
159
160
161
  
}

#endif // HAVE_UMFPACK

#endif // AMDIS_UMFPACKSOLVER_H