OEMSolver.h 8.4 KB
Newer Older
1
2
3
4
5
6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9
10
11
// ==  Institut fr Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12
13
14
15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
// ==                                                                        ==
// ============================================================================

/** \file OEMSolver.h */

/**
 * \defgroup Solver Solver module
 * @{ <img src="solver.png"> @}
 * 
 * \brief
 * Contains all classes needed for solving linear and non linear equation
 * systems.
 */

#ifndef AMDIS_OEM_SOLVER_H
#define AMDIS_OEM_SOLVER_H

#include "Global.h"
35
36
37
38
39
40
41
#include "AMDiS_fwd.h"
#include "Parameters.h"
#include "SolverMatrix.h"
#include "DOFVector.h"
#include "SystemVector.h"
#include "DOFMatrix.h"
#include "ITL_Preconditioner.h"
42
43
44
45
46
47
48
49
50
51
52
53

namespace AMDiS {

  /**
   * \ingroup Solver
   * 
   *\brief
   * Solver for linear equation systems.
   */
  class OEMSolver
  {
  public:
54
55
56
57

    typedef DOFMatrix::base_matrix_type   matrix_type;
    typedef DOFMatrix::value_type         value_type;

58
    /// The constructor reads needed parameters and sets solvers \ref name.
59
60
61
62
63
64
65
66
    OEMSolver(std::string str) 
      : name(str),
	tolerance(0),
	relative(0),
	max_iter(1000),
	info(0),
	residual(0),
	print_cycle(100),
67
68
	iterations(-1),
	error(-1),
69
	multipleRhs(false),
70
71
72
	leftPrecon(NULL),
	rightPrecon(NULL)
    {}
73

74
    ///
75
    virtual ~OEMSolver() 
Thomas Witkowski's avatar
Thomas Witkowski committed
76
    {}
77

78
79
80
    void initParameters()
    {
      FUNCNAME("OEMSolver::initParameters()");
81
      
82
83
84
85
86
87
      GET_PARAMETER(0, name + "->tolerance", "%f", &tolerance);
      GET_PARAMETER(0, name + "->relative tolerance", "%f", &relative);
      GET_PARAMETER(0, name + "->max iteration", "%d", &max_iter);
      GET_PARAMETER(0, name + "->print cycle", "%d", &print_cycle);
      GET_PARAMETER(0, name + "->info", "%d", &info);
    }
88

89
90
91
92
93
94
95
    /** Set left Preconditioner
     *
     *  It is assumed that the solver owns the preconditioner from this call.
     *  That means that the preconditioner is deleted by the solver when not needed.
     */
    void setLeftPrecon(ITL_BasePreconditioner* p) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
96
      leftPrecon = p;
97
    }
98

99
100
101
102
103
104
105
    /** Set right Preconditioner
     *
     *  It is assumed that the solver owns the preconditioner from this call.
     *  That means that the preconditioner is deleted by the solver when not needed.
     */
    void setRightPrecon(ITL_BasePreconditioner* p) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
106
      rightPrecon = p;
107
    }
108

109
    /// Linear System to be solved in the derived class
Thomas Witkowski's avatar
Thomas Witkowski committed
110
111
112
    virtual int solveSystem(const DOFMatrix::base_matrix_type& A,
			    mtl::dense_vector<value_type>& x,
			    const mtl::dense_vector<value_type>& b) = 0;
113
    
114

115
116
    /// Solve a linear system for a scalar problem.
    int solveSystem(const SolverMatrix<DOFMatrix>& A,
Thomas Witkowski's avatar
Thomas Witkowski committed
117
118
		    DOFVector<double>& x, 
		    DOFVector<double>& b) 
119
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
120
      mtl::dense_vector<value_type> xx(x.getUsedSize()), bb(b.getUsedSize());
Thomas Witkowski's avatar
Thomas Witkowski committed
121
122
123
124

      // Copy rhs vector
      int counter = 0;
      DOFVector<double>::Iterator it_b(&b, USED_DOFS);
125
      DOFVector<double>::Iterator it_x(&x, USED_DOFS);
126

127
128
129
130
131
      for (it_b.reset(), it_x.reset(); !it_b.end(); ++it_b, ++it_x) {	
	bb[counter] = *it_b;
	xx[counter] = *it_x;
	counter++;
      }
132

Thomas Witkowski's avatar
Thomas Witkowski committed
133
134
135
136
137
138
139
140
      int r = solveSystem(A.getMatrix(), xx, bb);

      // Copy solution vector to DOFVector
      counter = 0;
      for (it_x.reset(); !it_x.end(); ++it_x)
	*it_x = xx[counter++];

      return r;
141
    }
142

143
    /// Solve a linear system for a vectorial problem.
144
145
146
    virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
			    SystemVector& x, 
			    SystemVector& b) 
147
    {
148
149
      int ns = x.getSize(),  // Number of systems.
	size = x.getUsedSize(); // Size of all DOFVectors
150
151
152
153

      // Copy vectors
      mtl::dense_vector<value_type>   xx(size), bb(size);
      for (int i = 0, counter = 0; i < ns; i++) {
154
155
156
157
158
159
160
161
	DOFVector<double>::Iterator it_b(b.getDOFVector(i), USED_DOFS);	  
	DOFVector<double>::Iterator it_x(x.getDOFVector(i), USED_DOFS);
	
	for (it_b.reset(), it_x.reset(); !it_b.end(); ++it_b, ++it_x) {	
	  bb[counter] = *it_b;
	  xx[counter] = *it_x;
	  counter++;
	}
162
      }
163
      
164
      // Solver on DOFVector for single system
165
      int r = solveSystem(A.getMatrix(), xx, bb);
166
      
167
      for (int i = 0, counter = 0; i < ns; i++) {
168
169
170
	DOFVector<double>::Iterator it(x.getDOFVector(i), USED_DOFS);	  
	for (it.reset(); !it.end(); ++it)	
	  *it = xx[counter++];	
171
      }
172
      
173
      return r;
174
    }
175
176
177
178
179

    /** \name getting methods
     * \{ 
     */

180
    /// Returns solvers \ref name.
181
    inline std::string getName() const
182
    { 
183
      return name; 
184
    }
185

186
    /// Returns \ref tolerance
187
    inline double getTolerance() const
188
    {
189
      return tolerance;
190
    }  
191

192
    /// Returns \ref max_iter
193
    inline int getMaxIterations() const
194
    {
195
      return max_iter;
196
    }
197

198
    /// Returns \ref residual
199
    inline double getResidual() const
200
    {
201
      return residual;
202
    }
203
204
205
206
207
208
209

    /** \} */ 

    /** \name setting methods
     * \{ 
     */

210
    /// Sets \ref tolerance
211
212
    inline void setTolerance(double tol) 
    {
213
      tolerance = tol;
214
    }
215

216
    /// Sets \ref relative
Thomas Witkowski's avatar
Thomas Witkowski committed
217
    inline void setRelative(double rel) 
218
    {
219
      relative = rel;
220
    }
221

222
    /// Sets \ref max_iter
223
224
    inline void setMaxIterations(int i) 
    {
225
      max_iter = i;
226
    }
227

228
229
230
    /// Returns number of iterations in last run of an iterative solver
    inline int getIterations() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
231
      return iterations;
232
233
234
235
236
    }

    /// Returns error code in last run of an iterative solver
    inline int getErrorCode() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
237
      return error;
238
239
    }

240
    /// Sets \ref info
241
242
    inline void setInfo(int i) 
    {
243
      info = i;
244
    }
245

246
247
248
249
250
    void setMultipleRhs(bool b)
    {
      multipleRhs = b;
    }

251
252
253
    /** \} */

  protected:
254
    /// solvers name.
255
    std::string name;
256

257
    /// Solver tolerance |r|. Set in OEMSolver's constructor. 
258
259
    double tolerance;

260
    /// Relative solver tolerance |r|/|r0|. Set in OEMSolver's constructor. 
261
    double relative;
262

263
    /// maximal number of iterations. Set in OEMSolver's constructor.
264
265
    int max_iter;

266
    /// info level during solving the system. Set in OEMSolver's constructor.
267
268
    int info;

269
    /// current residual.
270
271
    double residual;

272
273
    /// Print cycle, after how many iterations the residuum norm is logged.
    int print_cycle;
274

275
276
277
278
279
280
    /// How many iterations were performed in last solver (not set by UmfPack)
    int iterations;

    /// Error code  in last solver (not set by UmfPack)
    int error;

281
282
283
284
285
286
287
    /** \brief
     * If true, the solver has to solve multiple right hand sides with the same
     * matrix. Some solvers, e.g. direct once, may take advanteges from this knowledge,
     * as they can do the factorization of the matrix only once.
     */
    bool multipleRhs;

288
    /// Left preconditioner (ignored by some iterative solvers and by UmfPack)
289
    ITL_BasePreconditioner *leftPrecon;
290

291
    /// Right preconditioner (ignored by some iterative solvers and by UmfPack)
292
    ITL_BasePreconditioner *rightPrecon;
293
294
295
296
297
298
299
300
  };

  /** 
   * \ingroup Solver
   *
   * \brief
   * Interface for creators of concrete OEMSolvers. 
   */
301
  class OEMSolverCreator : public CreatorInterface<OEMSolver>
302
303
  { 
  public:
304
    virtual ~OEMSolverCreator() {}
305

306
    /// Sets \ref problem
307
308
    void setName(std::string str) 
    { 
309
      name = str; 
310
    }
311
312
313
314
315
316

  protected:
    /** \brief
     * Pointer to the problem the solver will belong to. Needed as parameter
     * when constructing an OEMSolver
     */
317
    std::string name;
318
319
320
  };
}

321
#include "ITL_Solver.h"
322
323

#endif // AMDIS_OEM_SOLVER_H