UmfPackSolver.h 3.63 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
57
      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
58
    }
59

60
61
62
    template< typename Vector> 
    int solve(const Matrix& A, Vector& x, Vector& b) 
    { 
63
64
65
     FUNCNAME("(UmfPackSolver.h)::solve()");
     if (!solver) {
        if (!symmetric_strategy)
66
67
68
69
	  solver = new mtl::matrix::umfpack::solver<matrix_type>(A);
	else
	  solver = new mtl::matrix::umfpack::solver<matrix_type>(A, UMFPACK_STRATEGY_SYMMETRIC);
      } else {
70
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
81
      if( oem.getInfo() > 0) {
        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
	if (residual > oem.getTolerance()) {
88
	    ERROR_EXIT("Tolerance tol=%e could not be reached!\n Set tolerance by '->solver->tolerance:' \n", oem.getTolerance());	  
89
90
	}
      }
91
      return code;
92

93
    }
94

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

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

104
105
    int store_symbolic;
    int symmetric_strategy;
106
  };
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135

  using namespace MTLTypes;
  class UmfPackSolver : public MTL4Solver<MTLMatrix, MTLVector, Umfpack_runner< MTLMatrix > >
  {
    protected:

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

#endif // HAVE_UMFPACK

#endif // AMDIS_UMFPACKSOLVER_H