MFEM  v3.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
solvers.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.googlecode.com.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #ifndef MFEM_SOLVERS
13 #define MFEM_SOLVERS
14 
15 #include "../config/config.hpp"
16 #include "operator.hpp"
17 
18 #ifdef MFEM_USE_MPI
19 #include <mpi.h>
20 #endif
21 
22 #ifdef MFEM_USE_SUITESPARSE
23 #include <umfpack.h>
24 #endif
25 
26 namespace mfem
27 {
28 
30 class IterativeSolver : public Solver
31 {
32 #ifdef MFEM_USE_MPI
33 private:
34  int dot_prod_type; // 0 - local, 1 - global over 'comm'
35  MPI_Comm comm;
36 #endif
37 
38 protected:
39  const Operator *oper;
41 
43  double rel_tol, abs_tol;
44 
45  // stats
46  mutable int final_iter, converged;
47  mutable double final_norm;
48 
49  double Dot(const Vector &x, const Vector &y) const;
50  double Norm(const Vector &x) const { return sqrt(Dot(x, x)); }
51 
52 public:
54 
55 #ifdef MFEM_USE_MPI
56  IterativeSolver(MPI_Comm _comm);
57 #endif
58 
59  void SetRelTol(double rtol) { rel_tol = rtol; }
60  void SetAbsTol(double atol) { abs_tol = atol; }
61  void SetMaxIter(int max_it) { max_iter = max_it; }
62  void SetPrintLevel(int print_lvl);
63 
64  int GetNumIterations() { return final_iter; }
65  int GetConverged() { return converged; }
66  double GetFinalNorm() { return final_norm; }
67 
69  virtual void SetPreconditioner(Solver &pr);
70 
72  virtual void SetOperator(const Operator &op);
73 };
74 
75 
77 class SLISolver : public IterativeSolver
78 {
79 protected:
80  mutable Vector r, z;
81 
82  void UpdateVectors();
83 
84 public:
85  SLISolver() { }
86 
87 #ifdef MFEM_USE_MPI
88  SLISolver(MPI_Comm _comm) : IterativeSolver(_comm) { }
89 #endif
90 
91  virtual void SetOperator(const Operator &op)
93 
94  virtual void Mult(const Vector &x, Vector &y) const;
95 };
96 
98 void SLI(const Operator &A, const Vector &b, Vector &x,
99  int print_iter = 0, int max_num_iter = 1000,
100  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
101 
103 void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
104  int print_iter = 0, int max_num_iter = 1000,
105  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
106 
107 
109 class CGSolver : public IterativeSolver
110 {
111 protected:
112  mutable Vector r, d, z;
113 
114  void UpdateVectors();
115 
116 public:
117  CGSolver() { }
118 
119 #ifdef MFEM_USE_MPI
120  CGSolver(MPI_Comm _comm) : IterativeSolver(_comm) { }
121 #endif
122 
123  virtual void SetOperator(const Operator &op)
125 
126  virtual void Mult(const Vector &x, Vector &y) const;
127 };
128 
130 void CG(const Operator &A, const Vector &b, Vector &x,
131  int print_iter = 0, int max_num_iter = 1000,
132  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
133 
135 void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x,
136  int print_iter = 0, int max_num_iter = 1000,
137  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
138 
139 
142 {
143 protected:
144  int m;
145 
146 public:
147  GMRESSolver() { m = 50; }
148 
149 #ifdef MFEM_USE_MPI
150  GMRESSolver(MPI_Comm _comm) : IterativeSolver(_comm) { m = 50; }
151 #endif
152 
153  void SetKDim(int dim) { m = dim; }
154 
155  virtual void Mult(const Vector &x, Vector &y) const;
156 };
157 
160 {
161 protected:
162  int m;
163 
164 public:
165  FGMRESSolver() { m = 50; }
166 
167 #ifdef MFEM_USE_MPI
168  FGMRESSolver(MPI_Comm _comm) : IterativeSolver(_comm) { m = 50; }
169 #endif
170 
171  void SetKDim(int dim) { m = dim; }
172 
173  virtual void Mult(const Vector &x, Vector &y) const;
174 };
175 
177 int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M,
178  int &max_iter, int m, double &tol, double atol, int printit);
179 
181 void GMRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
182  int print_iter = 0, int max_num_iter = 1000, int m = 50,
183  double rtol = 1e-12, double atol = 1e-24);
184 
185 
188 {
189 protected:
190  mutable Vector p, phat, s, shat, t, v, r, rtilde;
191 
192  void UpdateVectors();
193 
194 public:
196 
197 #ifdef MFEM_USE_MPI
198  BiCGSTABSolver(MPI_Comm _comm) : IterativeSolver(_comm) { }
199 #endif
200 
201  virtual void SetOperator(const Operator &op)
203 
204  virtual void Mult(const Vector &x, Vector &y) const;
205 };
206 
208 int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M,
209  int &max_iter, double &tol, double atol, int printit);
210 
212 void BiCGSTAB(const Operator &A, Solver &B, const Vector &b, Vector &x,
213  int print_iter = 0, int max_num_iter = 1000,
214  double rtol = 1e-12, double atol = 1e-24);
215 
216 
219 {
220 protected:
221  mutable Vector v0, v1, w0, w1, q;
222  mutable Vector u1; // used in the preconditioned version
223 
224 public:
226 
227 #ifdef MFEM_USE_MPI
228  MINRESSolver(MPI_Comm _comm) : IterativeSolver(_comm) { }
229 #endif
230 
231  virtual void SetPreconditioner(Solver &pr)
233 
234  virtual void SetOperator(const Operator &op);
235 
236  virtual void Mult(const Vector &b, Vector &x) const;
237 };
238 
240 void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it = 0,
241  int max_it = 1000, double rtol = 1e-12, double atol = 1e-24);
242 
244 void MINRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
245  int print_it = 0, int max_it = 1000,
246  double rtol = 1e-12, double atol = 1e-24);
247 
248 
255 {
256 protected:
257  mutable Vector r, c;
258 
259 public:
261 
262 #ifdef MFEM_USE_MPI
263  NewtonSolver(MPI_Comm _comm) : IterativeSolver(_comm) { }
264 #endif
265  virtual void SetOperator(const Operator &op);
266 
267  void SetSolver(Solver &solver) { prec = &solver; }
268 
269  virtual void Mult(const Vector &b, Vector &x) const;
270 };
271 
272 
277 int aGMRES(const Operator &A, Vector &x, const Vector &b,
278  const Operator &M, int &max_iter,
279  int m_max, int m_min, int m_step, double cf,
280  double &tol, double &atol, int printit);
281 
282 
290 {
291 protected:
293  double a;
294 
296  inline double solve(double l, const Vector &xt, Vector &x, int &nclip) const
297  {
298  add(xt, l, w, x);
299  x.median(lo,hi);
300  nclip++;
301  return Dot(w,x)-a;
302  }
303 
304  inline void print_iteration(int it, double r, double l) const;
305 
306 public:
308 
309 #ifdef MFEM_USE_MPI
310  SLBQPOptimizer(MPI_Comm _comm) : IterativeSolver(_comm) {}
311 #endif
312 
313  void SetBounds(const Vector &_lo, const Vector &_hi);
314  void SetLinearConstraint(const Vector &_w, double _a);
315 
316  // For this problem type, we let the target values play the role of the
317  // initial vector xt, from which the operator generates the optimal vector x.
318  virtual void Mult(const Vector &xt, Vector &x) const;
319 
321  virtual void SetPreconditioner(Solver &pr);
322  virtual void SetOperator(const Operator &op);
323 };
324 
325 
326 #ifdef MFEM_USE_SUITESPARSE
327 
329 class UMFPackSolver : public Solver
330 {
331 protected:
334  void *Numeric;
335  SuiteSparse_long *AI, *AJ;
336 
337  void Init();
338 
339 public:
340  double Control[UMFPACK_CONTROL];
341  mutable double Info[UMFPACK_INFO];
342 
343  UMFPackSolver(bool _use_long_ints = false)
344  : use_long_ints(_use_long_ints) { Init(); }
345  UMFPackSolver(SparseMatrix &A, bool _use_long_ints = false)
346  : use_long_ints(_use_long_ints) { Init(); SetOperator(A); }
347 
348  // Works on sparse matrices only; calls SparseMatrix::SortColumnIndices().
349  virtual void SetOperator(const Operator &op);
350 
351  void SetPrintLevel(int print_lvl) { Control[UMFPACK_PRL] = print_lvl; }
352 
353  virtual void Mult(const Vector &b, Vector &x) const;
354  virtual void MultTranspose(const Vector &b, Vector &x) const;
355 
356  virtual ~UMFPackSolver();
357 };
358 
359 #endif // MFEM_USE_SUITESPARSE
360 
361 }
362 
363 #endif // MFEM_SOLVERS
void SetSolver(Solver &solver)
Definition: solvers.hpp:267
Conjugate gradient method.
Definition: solvers.hpp:109
virtual void MultTranspose(const Vector &b, Vector &x) const
Action of the transpose operator.
Definition: solvers.cpp:1629
double Info[UMFPACK_INFO]
Definition: solvers.hpp:341
SparseMatrix * mat
Definition: solvers.hpp:333
const Operator * oper
Definition: solvers.hpp:39
SuiteSparse_long * AI
Definition: solvers.hpp:335
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:248
NewtonSolver(MPI_Comm _comm)
Definition: solvers.hpp:263
FGMRES method.
Definition: solvers.hpp:159
UMFPackSolver(SparseMatrix &A, bool _use_long_ints=false)
Definition: solvers.hpp:345
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: solvers.cpp:285
MINRES method.
Definition: solvers.hpp:218
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: solvers.cpp:112
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:227
double Norm(const Vector &x) const
Definition: solvers.hpp:50
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: solvers.cpp:797
int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, double &tol, double atol, int printit)
BiCGSTAB method. (tolerances are squared)
Definition: solvers.cpp:911
void SetKDim(int dim)
Definition: solvers.hpp:171
FGMRESSolver(MPI_Comm _comm)
Definition: solvers.hpp:168
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:329
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: solvers.cpp:485
UMFPackSolver(bool _use_long_ints=false)
Definition: solvers.hpp:343
virtual void Mult(const Vector &b, Vector &x) const
Operator application.
Definition: solvers.cpp:946
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: solvers.cpp:628
void SetLinearConstraint(const Vector &_w, double _a)
Definition: solvers.cpp:1303
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:231
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:71
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:96
double solve(double l, const Vector &xt, Vector &x, int &nclip) const
Solve QP at fixed lambda.
Definition: solvers.hpp:296
Data type sparse matrix.
Definition: sparsemat.hpp:38
void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it, int max_it, double rtol, double atol)
MINRES method without preconditioner. (tolerances are squared)
Definition: solvers.cpp:1074
virtual ~UMFPackSolver()
Definition: solvers.cpp:1663
void SetMaxIter(int max_it)
Definition: solvers.hpp:61
void median(const Vector &lo, const Vector &hi)
v = median(v,lo,hi) entrywise. Implementation assumes lo &lt;= hi.
Definition: vector.cpp:390
SuiteSparse_long * AJ
Definition: solvers.hpp:335
void SetKDim(int dim)
Definition: solvers.hpp:153
void CG(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Conjugate gradient method. (tolerances are squared)
Definition: solvers.cpp:411
void print_iteration(int it, double r, double l) const
Definition: solvers.cpp:1321
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition: solvers.cpp:424
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:934
double Dot(const Vector &x, const Vector &y) const
Definition: solvers.cpp:50
BiCGSTABSolver(MPI_Comm _comm)
Definition: solvers.hpp:198
SLBQPOptimizer(MPI_Comm _comm)
Definition: solvers.hpp:310
virtual void Mult(const Vector &b, Vector &x) const
Operator application.
Definition: solvers.cpp:1111
CGSolver(MPI_Comm _comm)
Definition: solvers.hpp:120
void SetPrintLevel(int print_lvl)
Definition: solvers.hpp:351
Stationary linear iteration: x &lt;- x + B (b - A x)
Definition: solvers.hpp:77
void SetAbsTol(double atol)
Definition: solvers.hpp:60
void SetRelTol(double rtol)
Definition: solvers.hpp:59
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:340
Abstract base class for iterative solver.
Definition: solvers.hpp:30
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1315
SLISolver(MPI_Comm _comm)
Definition: solvers.hpp:88
GMRES method.
Definition: solvers.hpp:141
GMRESSolver(MPI_Comm _comm)
Definition: solvers.hpp:150
void SetBounds(const Vector &_lo, const Vector &_hi)
Definition: solvers.cpp:1297
void UpdateVectors()
Definition: solvers.cpp:106
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:201
virtual void Mult(const Vector &xt, Vector &x) const
Operator application.
Definition: solvers.cpp:1328
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:123
Vector data type.
Definition: vector.hpp:29
MINRESSolver(MPI_Comm _comm)
Definition: solvers.hpp:228
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:90
int aGMRES(const Operator &A, Vector &x, const Vector &b, const Operator &M, int &max_iter, int m_max, int m_min, int m_step, double cf, double &tol, double &atol, int printit)
Definition: solvers.cpp:1168
void UpdateVectors()
Definition: solvers.cpp:278
BiCGSTAB method.
Definition: solvers.hpp:187
Base class for solvers.
Definition: operator.hpp:102
double GetFinalNorm()
Definition: solvers.hpp:66
virtual void SetPreconditioner(Solver &pr)
These are not currently meaningful for this solver and will error out.
Definition: solvers.cpp:1309
Abstract operator.
Definition: operator.hpp:21
virtual void Mult(const Vector &b, Vector &x) const
Operator application.
Definition: solvers.cpp:1597
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:91
void SLI(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Stationary linear iteration. (tolerances are squared)
Definition: solvers.cpp:250
int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M, int &max_iter, int m, double &tol, double atol, int printit)
GMRES method. (tolerances are squared)
Definition: solvers.cpp:761
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1100
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:1505