MFEM  v3.4
Finite element discretization library
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.org.
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 "sparsemat.hpp"
24 #include <umfpack.h>
25 #include <klu.h>
26 #endif
27 
28 namespace mfem
29 {
30 
31 /// Abstract base class for iterative solver
32 class IterativeSolver : public Solver
33 {
34 #ifdef MFEM_USE_MPI
35 private:
36  int dot_prod_type; // 0 - local, 1 - global over 'comm'
37  MPI_Comm comm;
38 #endif
39 
40 protected:
41  const Operator *oper;
43 
45  double rel_tol, abs_tol;
46 
47  // stats
48  mutable int final_iter, converged;
49  mutable double final_norm;
50 
51  double Dot(const Vector &x, const Vector &y) const;
52  double Norm(const Vector &x) const { return sqrt(Dot(x, x)); }
53 
54 public:
56 
57 #ifdef MFEM_USE_MPI
58  IterativeSolver(MPI_Comm _comm);
59 #endif
60 
61  void SetRelTol(double rtol) { rel_tol = rtol; }
62  void SetAbsTol(double atol) { abs_tol = atol; }
63  void SetMaxIter(int max_it) { max_iter = max_it; }
64  void SetPrintLevel(int print_lvl);
65 
66  int GetNumIterations() const { return final_iter; }
67  int GetConverged() const { return converged; }
68  double GetFinalNorm() const { return final_norm; }
69 
70  /// This should be called before SetOperator
71  virtual void SetPreconditioner(Solver &pr);
72 
73  /// Also calls SetOperator for the preconditioner if there is one
74  virtual void SetOperator(const Operator &op);
75 };
76 
77 
78 /// Stationary linear iteration: x <- x + B (b - A x)
79 class SLISolver : public IterativeSolver
80 {
81 protected:
82  mutable Vector r, z;
83 
84  void UpdateVectors();
85 
86 public:
87  SLISolver() { }
88 
89 #ifdef MFEM_USE_MPI
90  SLISolver(MPI_Comm _comm) : IterativeSolver(_comm) { }
91 #endif
92 
93  virtual void SetOperator(const Operator &op)
95 
96  virtual void Mult(const Vector &b, Vector &x) const;
97 };
98 
99 /// Stationary linear iteration. (tolerances are squared)
100 void SLI(const Operator &A, const Vector &b, Vector &x,
101  int print_iter = 0, int max_num_iter = 1000,
102  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
103 
104 /// Preconditioned stationary linear iteration. (tolerances are squared)
105 void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
106  int print_iter = 0, int max_num_iter = 1000,
107  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
108 
109 
110 /// Conjugate gradient method
111 class CGSolver : public IterativeSolver
112 {
113 protected:
114  mutable Vector r, d, z;
115 
116  void UpdateVectors();
117 
118 public:
119  CGSolver() { }
120 
121 #ifdef MFEM_USE_MPI
122  CGSolver(MPI_Comm _comm) : IterativeSolver(_comm) { }
123 #endif
124 
125  virtual void SetOperator(const Operator &op)
127 
128  virtual void Mult(const Vector &b, Vector &x) const;
129 };
130 
131 /// Conjugate gradient method. (tolerances are squared)
132 void CG(const Operator &A, const Vector &b, Vector &x,
133  int print_iter = 0, int max_num_iter = 1000,
134  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
135 
136 /// Preconditioned conjugate gradient method. (tolerances are squared)
137 void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x,
138  int print_iter = 0, int max_num_iter = 1000,
139  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
140 
141 
142 /// GMRES method
144 {
145 protected:
146  int m; // see SetKDim()
147 
148 public:
149  GMRESSolver() { m = 50; }
150 
151 #ifdef MFEM_USE_MPI
152  GMRESSolver(MPI_Comm _comm) : IterativeSolver(_comm) { m = 50; }
153 #endif
154 
155  /// Set the number of iteration to perform between restarts, default is 50.
156  void SetKDim(int dim) { m = dim; }
157 
158  virtual void Mult(const Vector &b, Vector &x) const;
159 };
160 
161 /// FGMRES method
163 {
164 protected:
165  int m;
166 
167 public:
168  FGMRESSolver() { m = 50; }
169 
170 #ifdef MFEM_USE_MPI
171  FGMRESSolver(MPI_Comm _comm) : IterativeSolver(_comm) { m = 50; }
172 #endif
173 
174  void SetKDim(int dim) { m = dim; }
175 
176  virtual void Mult(const Vector &b, Vector &x) const;
177 };
178 
179 /// GMRES method. (tolerances are squared)
180 int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M,
181  int &max_iter, int m, double &tol, double atol, int printit);
182 
183 /// GMRES method. (tolerances are squared)
184 void GMRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
185  int print_iter = 0, int max_num_iter = 1000, int m = 50,
186  double rtol = 1e-12, double atol = 1e-24);
187 
188 
189 /// BiCGSTAB method
191 {
192 protected:
193  mutable Vector p, phat, s, shat, t, v, r, rtilde;
194 
195  void UpdateVectors();
196 
197 public:
199 
200 #ifdef MFEM_USE_MPI
201  BiCGSTABSolver(MPI_Comm _comm) : IterativeSolver(_comm) { }
202 #endif
203 
204  virtual void SetOperator(const Operator &op)
206 
207  virtual void Mult(const Vector &b, Vector &x) const;
208 };
209 
210 /// BiCGSTAB method. (tolerances are squared)
211 int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M,
212  int &max_iter, double &tol, double atol, int printit);
213 
214 /// BiCGSTAB method. (tolerances are squared)
215 void BiCGSTAB(const Operator &A, Solver &B, const Vector &b, Vector &x,
216  int print_iter = 0, int max_num_iter = 1000,
217  double rtol = 1e-12, double atol = 1e-24);
218 
219 
220 /// MINRES method
222 {
223 protected:
224  mutable Vector v0, v1, w0, w1, q;
225  mutable Vector u1; // used in the preconditioned version
226 
227 public:
229 
230 #ifdef MFEM_USE_MPI
231  MINRESSolver(MPI_Comm _comm) : IterativeSolver(_comm) { }
232 #endif
233 
234  virtual void SetPreconditioner(Solver &pr)
235  {
237  if (oper) { u1.SetSize(width); }
238  }
239 
240  virtual void SetOperator(const Operator &op);
241 
242  virtual void Mult(const Vector &b, Vector &x) const;
243 };
244 
245 /// MINRES method without preconditioner. (tolerances are squared)
246 void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it = 0,
247  int max_it = 1000, double rtol = 1e-12, double atol = 1e-24);
248 
249 /// MINRES method with preconditioner. (tolerances are squared)
250 void MINRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
251  int print_it = 0, int max_it = 1000,
252  double rtol = 1e-12, double atol = 1e-24);
253 
254 
255 /// Newton's method for solving F(x)=b for a given operator F.
256 /** The method GetGradient() must be implemented for the operator F.
257  The preconditioner is used (in non-iterative mode) to evaluate
258  the action of the inverse gradient of the operator. */
260 {
261 protected:
262  mutable Vector r, c;
263 
264 public:
266 
267 #ifdef MFEM_USE_MPI
268  NewtonSolver(MPI_Comm _comm) : IterativeSolver(_comm) { }
269 #endif
270  virtual void SetOperator(const Operator &op);
271 
272  /// Set the linear solver for inverting the Jacobian.
273  /** This method is equivalent to calling SetPreconditioner(). */
274  virtual void SetSolver(Solver &solver) { prec = &solver; }
275 
276  /// Solve the nonlinear system with right-hand side @a b.
277  /** If `b.Size() != Height()`, then @a b is assumed to be zero. */
278  virtual void Mult(const Vector &b, Vector &x) const;
279 
280  /** @brief This method can be overloaded in derived classes to implement line
281  search algorithms. */
282  /** The base class implementation (NewtonSolver) simply returns 1. A return
283  value of 0 indicates a failure, interrupting the Newton iteration. */
284  virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const
285  { return 1.0; }
286 };
287 
288 /** Adaptive restarted GMRES.
289  m_max and m_min(=1) are the maximal and minimal restart parameters.
290  m_step(=1) is the step to use for going from m_max and m_min.
291  cf(=0.4) is a desired convergence factor. */
292 int aGMRES(const Operator &A, Vector &x, const Vector &b,
293  const Operator &M, int &max_iter,
294  int m_max, int m_min, int m_step, double cf,
295  double &tol, double &atol, int printit);
296 
297 
298 /** SLBQP: (S)ingle (L)inearly Constrained with (B)ounds (Q)uadratic (P)rogram
299 
300  minimize 1/2 ||x - x_t||^2, subject to:
301  lo_i <= x_i <= hi_i
302  sum_i w_i x_i = a
303 */
305 {
306 protected:
308  double a;
309 
310  /// Solve QP at fixed lambda
311  inline double solve(double l, const Vector &xt, Vector &x, int &nclip) const
312  {
313  add(xt, l, w, x);
314  x.median(lo,hi);
315  nclip++;
316  return Dot(w,x)-a;
317  }
318 
319  inline void print_iteration(int it, double r, double l) const;
320 
321 public:
323 
324 #ifdef MFEM_USE_MPI
325  SLBQPOptimizer(MPI_Comm _comm) : IterativeSolver(_comm) {}
326 #endif
327 
328  void SetBounds(const Vector &_lo, const Vector &_hi);
329  void SetLinearConstraint(const Vector &_w, double _a);
330 
331  // For this problem type, we let the target values play the role of the
332  // initial vector xt, from which the operator generates the optimal vector x.
333  virtual void Mult(const Vector &xt, Vector &x) const;
334 
335  /// These are not currently meaningful for this solver and will error out.
336  virtual void SetPreconditioner(Solver &pr);
337  virtual void SetOperator(const Operator &op);
338 };
339 
340 
341 #ifdef MFEM_USE_SUITESPARSE
342 
343 /// Direct sparse solver using UMFPACK
344 class UMFPackSolver : public Solver
345 {
346 protected:
349  void *Numeric;
350  SuiteSparse_long *AI, *AJ;
351 
352  void Init();
353 
354 public:
355  double Control[UMFPACK_CONTROL];
356  mutable double Info[UMFPACK_INFO];
357 
358  /** @brief For larger matrices, if the solver fails, set the parameter @a
359  _use_long_ints = true. */
360  UMFPackSolver(bool _use_long_ints = false)
361  : use_long_ints(_use_long_ints) { Init(); }
362  /** @brief Factorize the given SparseMatrix using the defaults. For larger
363  matrices, if the solver fails, set the parameter @a _use_long_ints =
364  true. */
365  UMFPackSolver(SparseMatrix &A, bool _use_long_ints = false)
366  : use_long_ints(_use_long_ints) { Init(); SetOperator(A); }
367 
368  /** @brief Factorize the given Operator @a op which must be a SparseMatrix.
369 
370  The factorization uses the parameters set in the #Control data member.
371  @note This method calls SparseMatrix::SortColumnIndices() with @a op,
372  modifying the matrix if the column indices are not already sorted. */
373  virtual void SetOperator(const Operator &op);
374 
375  /// Set the print level field in the #Control data member.
376  void SetPrintLevel(int print_lvl) { Control[UMFPACK_PRL] = print_lvl; }
377 
378  virtual void Mult(const Vector &b, Vector &x) const;
379  virtual void MultTranspose(const Vector &b, Vector &x) const;
380 
381  virtual ~UMFPackSolver();
382 };
383 
384 /// Direct sparse solver using KLU
385 class KLUSolver : public Solver
386 {
387 protected:
389  klu_symbolic *Symbolic;
390  klu_numeric *Numeric;
391 
392  void Init();
393 
394 public:
396  : mat(0),Symbolic(0),Numeric(0)
397  { Init(); }
399  : mat(0),Symbolic(0),Numeric(0)
400  { Init(); SetOperator(A); }
401 
402  // Works on sparse matrices only; calls SparseMatrix::SortColumnIndices().
403  virtual void SetOperator(const Operator &op);
404 
405  virtual void Mult(const Vector &b, Vector &x) const;
406  virtual void MultTranspose(const Vector &b, Vector &x) const;
407 
408  virtual ~KLUSolver();
409 
410  mutable klu_common Common;
411 };
412 
413 #endif // MFEM_USE_SUITESPARSE
414 
415 }
416 
417 #endif // MFEM_SOLVERS
Conjugate gradient method.
Definition: solvers.hpp:111
double Info[UMFPACK_INFO]
Definition: solvers.hpp:356
SparseMatrix * mat
Definition: solvers.hpp:348
const Operator * oper
Definition: solvers.hpp:41
SuiteSparse_long * AI
Definition: solvers.hpp:350
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:328
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:702
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1046
NewtonSolver(MPI_Comm _comm)
Definition: solvers.hpp:268
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:295
FGMRES method.
Definition: solvers.hpp:162
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:522
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:117
Direct sparse solver using KLU.
Definition: solvers.hpp:385
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: solvers.hpp:274
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1904
UMFPackSolver(SparseMatrix &A, bool _use_long_ints=false)
Factorize the given SparseMatrix using the defaults. For larger matrices, if the solver fails...
Definition: solvers.hpp:365
MINRES method.
Definition: solvers.hpp:221
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:260
double GetFinalNorm() const
Definition: solvers.hpp:68
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:1009
void SetKDim(int dim)
Definition: solvers.hpp:174
double Norm(const Vector &x) const
Definition: solvers.hpp:52
FGMRESSolver(MPI_Comm _comm)
Definition: solvers.hpp:171
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:344
UMFPackSolver(bool _use_long_ints=false)
For larger matrices, if the solver fails, set the parameter _use_long_ints = true.
Definition: solvers.hpp:360
void SetLinearConstraint(const Vector &_w, double _a)
Definition: solvers.cpp:1476
int dim
Definition: ex3.cpp:47
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:234
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:72
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:99
virtual void MultTranspose(const Vector &b, Vector &x) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: solvers.cpp:1918
Data type sparse matrix.
Definition: sparsemat.hpp:38
double solve(double l, const Vector &xt, Vector &x, int &nclip) const
Solve QP at fixed lambda.
Definition: solvers.hpp:311
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:1204
virtual ~UMFPackSolver()
Definition: solvers.cpp:1851
int GetNumIterations() const
Definition: solvers.hpp:66
void SetMaxIter(int max_it)
Definition: solvers.hpp:63
void median(const Vector &lo, const Vector &hi)
v = median(v,lo,hi) entrywise. Implementation assumes lo <= hi.
Definition: vector.cpp:441
SuiteSparse_long * AJ
Definition: solvers.hpp:350
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition: solvers.hpp:156
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:259
void print_iteration(int it, double r, double l) const
Definition: solvers.cpp:1494
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:444
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:457
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:1873
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1032
BiCGSTABSolver(MPI_Comm _comm)
Definition: solvers.hpp:201
klu_symbolic * Symbolic
Definition: solvers.hpp:389
SLBQPOptimizer(MPI_Comm _comm)
Definition: solvers.hpp:325
CGSolver(MPI_Comm _comm)
Definition: solvers.hpp:122
void SetPrintLevel(int print_lvl)
Set the print level field in the Control data member.
Definition: solvers.hpp:376
Stationary linear iteration: x <- x + B (b - A x)
Definition: solvers.hpp:79
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1241
void SetAbsTol(double atol)
Definition: solvers.hpp:62
void SetRelTol(double rtol)
Definition: solvers.hpp:61
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:355
Abstract base class for iterative solver.
Definition: solvers.hpp:32
virtual void MultTranspose(const Vector &b, Vector &x) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: solvers.cpp:1817
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1488
SLISolver(MPI_Comm _comm)
Definition: solvers.hpp:90
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1785
GMRES method.
Definition: solvers.hpp:143
GMRESSolver(MPI_Comm _comm)
Definition: solvers.hpp:152
void SetBounds(const Vector &_lo, const Vector &_hi)
Definition: solvers.cpp:1470
void UpdateVectors()
Definition: solvers.cpp:111
klu_common Common
Definition: solvers.hpp:410
virtual void Mult(const Vector &xt, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1501
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:204
double Dot(const Vector &x, const Vector &y) const
Definition: solvers.cpp:51
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:125
Vector data type.
Definition: vector.hpp:48
MINRESSolver(MPI_Comm _comm)
Definition: solvers.hpp:231
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:93
KLUSolver(SparseMatrix &A)
Definition: solvers.hpp:398
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:1318
virtual ~KLUSolver()
Definition: solvers.cpp:1932
void UpdateVectors()
Definition: solvers.cpp:288
BiCGSTAB method.
Definition: solvers.hpp:190
SparseMatrix * mat
Definition: solvers.hpp:388
Base class for solvers.
Definition: operator.hpp:268
virtual void SetPreconditioner(Solver &pr)
These are not currently meaningful for this solver and will error out.
Definition: solvers.cpp:1482
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:880
virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const
This method can be overloaded in derived classes to implement line search algorithms.
Definition: solvers.hpp:284
Abstract operator.
Definition: operator.hpp:21
klu_numeric * Numeric
Definition: solvers.hpp:390
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:93
int GetConverged() const
Definition: solvers.hpp:67
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:260
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:844
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1230
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:1688