UmfPackSolver.h 3.68 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
    template< typename Vector> 
    int solve(const Matrix& A, Vector& x, Vector& b) 
    { 
67
68
      FUNCNAME("Umfpack_runner::solve()");
      if (!solver) {
69
	solver = new mtl::matrix::umfpack::solver<matrix_type>(A, symmetric_strategy, alloc_init);
70
      } else {
71
	if (!oem.getMultipleRhs()) {
72
73
74
 	  if (store_symbolic)
 	    solver->update_numeric();
 	  else
75
76
	    solver->update();
	}
77
78
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
79
      int code = (*solver)(x, b);
80
      if (oem.getInfo() > 0) {
81
        mtl::dense_vector<typename matrix_type::value_type> r(b); 
82
        r -= A * x; 
83
84
85
        double residual = two_norm(r);
        oem.setResidual(residual);
        oem.setErrorCode(code);
86
        std::cout << "UmfPackSolver: ||b-Ax|| = " << residual << "\n";
87
	TEST_EXIT(residual <= oem.getTolerance())("Tolerance tol = %e could not be reached!\n Set tolerance by '->solver->tolerance:' \n", oem.getTolerance());	
88
      }
89
      return code;
90

91
    }
92

93
94
    ~Umfpack_runner() 
    {
95
      if (!solver)
96
97
98
	delete solver;
    }

99
  private:
100
    mtl::matrix::umfpack::solver<matrix_type> *solver;
Thomas Witkowski's avatar
Thomas Witkowski committed
101

102
    int store_symbolic;
103

104
    int symmetric_strategy;
105
106

    double alloc_init;
107
  };
108
109
110
111

  using namespace MTLTypes;
  class UmfPackSolver : public MTL4Solver<MTLMatrix, MTLVector, Umfpack_runner< MTLMatrix > >
  {
112
  protected:
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136

  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)
    {
    }
    
    
  private:
  };
137
138
139
140
141
142
  
}

#endif // HAVE_UMFPACK

#endif // AMDIS_UMFPACKSOLVER_H