UmfPackSolver.h 3.66 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
53
    typedef Matrix matrix_type ;
    OEMSolver& oem;
    Umfpack_runner(OEMSolver* oem_):
      oem(*oem_),
      solver(NULL),
      store_symbolic(0),
      symmetric_strategy(0)
54
    {
55
56
      FUNCNAME("Umfpack_runner::Umfpack_runner()");

57
58
59
      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);
Thomas Witkowski's avatar
Thomas Witkowski committed
60
    }
61

62
63
64
    template< typename Vector> 
    int solve(const Matrix& A, Vector& x, Vector& b) 
    { 
65
66
      FUNCNAME("Umfpack_runner::solve()");
      if (!solver) {
67
        if (!symmetric_strategy)
68
69
70
71
	  solver = new mtl::matrix::umfpack::solver<matrix_type>(A);
	else
	  solver = new mtl::matrix::umfpack::solver<matrix_type>(A, UMFPACK_STRATEGY_SYMMETRIC);
      } else {
72
	if (!oem.getMultipleRhs()) {
73
74
75
 	  if (store_symbolic)
 	    solver->update_numeric();
 	  else
76
77
	    solver->update();
	}
78
79
      }

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

92
    }
93

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

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

103
104
    int store_symbolic;
    int symmetric_strategy;
105
  };
106
107
108
109

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

  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:
  };
135
136
137
138
139
140
  
}

#endif // HAVE_UMFPACK

#endif // AMDIS_UMFPACKSOLVER_H