MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
solvers.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_SOLVERS
13 #define MFEM_SOLVERS
14 
15 #include "../config/config.hpp"
16 #include "densemat.hpp"
17 #include "handle.hpp"
18 #include <memory>
19 
20 #ifdef MFEM_USE_MPI
21 #include <mpi.h>
22 #endif
23 
24 #ifdef MFEM_USE_SUITESPARSE
25 #include "sparsemat.hpp"
26 #include <umfpack.h>
27 #include <klu.h>
28 #endif
29 
30 namespace mfem
31 {
32 
33 class BilinearForm;
34 
35 /// Abstract base class for an iterative solver monitor
37 {
38 protected:
39  /// The last IterativeSolver to which this monitor was attached.
41 
42 public:
44 
46 
47  /// Monitor the residual vector r
48  virtual void MonitorResidual(int it, double norm, const Vector &r,
49  bool final)
50  {
51  }
52 
53  /// Monitor the solution vector x
54  virtual void MonitorSolution(int it, double norm, const Vector &x,
55  bool final)
56  {
57  }
58 
59  /** @brief This method is invoked by IterativeSolver::SetMonitor, informing
60  the monitor which IterativeSolver is using it. */
61  void SetIterativeSolver(const IterativeSolver &solver)
62  { iter_solver = &solver; }
63 };
64 
65 /// Abstract base class for iterative solver
66 class IterativeSolver : public Solver
67 {
68 public:
69  /** @brief Settings for the output behavior of the IterativeSolver.
70 
71  By default, all output is suppressed. The construction of the desired
72  print level can be achieved through a builder pattern, for example
73 
74  PrintLevel().Errors().Warnings()
75 
76  constructs the print level with only errors and warnings enabled.
77  */
78  struct PrintLevel
79  {
80  /** @brief If a fatal problem has been detected the failure will be
81  reported to @ref mfem::err. */
82  bool errors = false;
83  /** @brief If a non-fatal problem has been detected some context-specific
84  information will be reported to @ref mfem::out */
85  bool warnings = false;
86  /** @brief Detailed information about each iteration will be reported to
87  @ref mfem::out */
88  bool iterations = false;
89  /** @brief A summary of the solver process will be reported after the last
90  iteration to @ref mfem::out */
91  bool summary = false;
92  /** @brief Information about the first and last iteration will be printed
93  to @ref mfem::out */
94  bool first_and_last = false;
95 
96  /// Initializes the print level to suppress
97  PrintLevel() = default;
98 
99  /** @name Builder
100  These methods are utilized to construct PrintLevel objects through a
101  builder approach by chaining the function calls in this group. */
102  ///@{
103  PrintLevel &None() { *this = PrintLevel(); return *this; }
104  PrintLevel &Warnings() { warnings=true; return *this; }
105  PrintLevel &Errors() { errors=true; return *this; }
106  PrintLevel &Iterations() { iterations=true; return *this; }
107  PrintLevel &FirstAndLast() { first_and_last=true; return *this; }
108  PrintLevel &Summary() { summary=true; return *this; }
110  ///@}
111  };
112 
113 #ifdef MFEM_USE_MPI
114 private:
115  int dot_prod_type; // 0 - local, 1 - global over 'comm'
116  MPI_Comm comm = MPI_COMM_NULL;
117 #endif
118 
119 protected:
120  const Operator *oper;
123 
124  /// @name Reporting (protected attributes and member functions)
125  ///@{
126 
127  /** @brief (DEPRECATED) Legacy print level definition, which is left for
128  compatibility with custom iterative solvers.
129  @deprecated #print_options should be used instead. */
130  int print_level = -1;
131 
132  /** @brief Output behavior for the iterative solver.
133 
134  This primarily controls the output behavior of the iterative solvers
135  provided by this library. This member must be synchronized with
136  #print_level to ensure compatibility with custom iterative solvers. */
138 
139  /// Convert a legacy print level integer to a PrintLevel object
141 
142  /// @brief Use some heuristics to guess a legacy print level corresponding to
143  /// the given PrintLevel.
144  static int GuessLegacyPrintLevel(PrintLevel);
145  ///@}
146 
147  /// @name Convergence (protected attributes)
148  ///@{
149 
150  /// Limit for the number of iterations the solver is allowed to do
151  int max_iter;
152 
153  /// Relative tolerance.
154  double rel_tol;
155 
156  /// Absolute tolerance.
157  double abs_tol;
158 
159  ///@}
160 
161  /// @name Solver statistics (protected attributes)
162  ///@{
163 
164  mutable int final_iter;
165  mutable bool converged;
166  mutable double final_norm;
167 
168  ///@}
169 
170  double Dot(const Vector &x, const Vector &y) const;
171  double Norm(const Vector &x) const { return sqrt(Dot(x, x)); }
172  void Monitor(int it, double norm, const Vector& r, const Vector& x,
173  bool final=false) const;
174 
175 public:
176  IterativeSolver();
177 
178 #ifdef MFEM_USE_MPI
179  IterativeSolver(MPI_Comm comm_);
180 #endif
181 
182  /** @name Convergence
183  @brief Termination criteria for the iterative solvers.
184 
185  @details While the convergence criterion is solver specific, most of the
186  provided iterative solvers use one of the following criteria
187 
188  \f$ ||r||_X \leq tol_{rel}||r_0||_X \f$,
189 
190  \f$ ||r||_X \leq tol_{abs} \f$,
191 
192  \f$ ||r||_X \leq \max\{ tol_{abs}, tol_{rel} ||r_0||_X \} \f$,
193 
194  where X denotes the space in which the norm is measured. The choice of
195  X depends on the specific iterative solver.
196  */
197  ///@{
198  void SetRelTol(double rtol) { rel_tol = rtol; }
199  void SetAbsTol(double atol) { abs_tol = atol; }
200  void SetMaxIter(int max_it) { max_iter = max_it; }
201  ///@}
202 
203  /** @name Reporting
204  These options control the internal reporting behavior into ::mfem::out
205  and ::mfem::err of the iterative solvers.
206  */
207  ///@{
208 
209  /// @brief Legacy method to set the level of verbosity of the solver output.
210  /** This is the old way to control what information will be printed to
211  ::mfem::out and ::mfem::err. The behavior for the print level for all
212  iterative solvers is:
213 
214  - -1: Suppress all outputs.
215  - 0: Print information about all detected issues (e.g. no convergence).
216  - 1: Same as level 0, but with detailed information about each
217  iteration.
218  - 2: Print detected issues and a summary when the solver terminates.
219  - 3: Same as 2, but print also the first and last iterations.
220  - >3: Custom print options which are dependent on the specific solver.
221 
222  In parallel, only rank 0 produces output.
223 
224  @note It is recommended to use @ref SetPrintLevel(PrintLevel) instead.
225 
226  @note Some derived classes, like KINSolver, redefine this method and use
227  their own set of print level constants. */
228  virtual void SetPrintLevel(int print_lvl);
229 
230  /// @brief Set the level of verbosity of the solver output.
231  /** In parallel, only rank 0 produces outputs. Errors are output to
232  ::mfem::err and all other information to ::mfem::out.
233 
234  @note Not all subclasses of IterativeSolver support all possible options.
235 
236  @note Some derived classes, like KINSolver, disable this method in favor
237  of SetPrintLevel(int).
238 
239  @sa PrintLevel for possible options.
240  */
241  virtual void SetPrintLevel(PrintLevel);
242  ///@}
243 
244  /// @name Solver statistics
245  ///@{
246  int GetNumIterations() const { return final_iter; }
247  bool GetConverged() const { return converged; }
248  double GetFinalNorm() const { return final_norm; }
249  ///@}
250 
251  /// This should be called before SetOperator
252  virtual void SetPreconditioner(Solver &pr);
253 
254  /// Also calls SetOperator for the preconditioner if there is one
255  virtual void SetOperator(const Operator &op) override;
256 
257  /// Set the iterative solver monitor
259  { monitor = &m; m.SetIterativeSolver(*this); }
260 
261 #ifdef MFEM_USE_MPI
262  /** @brief Return the associated MPI communicator, or MPI_COMM_NULL if no
263  communicator is set. */
264  MPI_Comm GetComm() const
265  { return dot_prod_type == 0 ? MPI_COMM_NULL : comm; }
266 #endif
267 };
268 
269 
270 /// Jacobi smoothing for a given bilinear form (no matrix necessary).
271 /** Useful with tensorized, partially assembled operators. Can also be defined
272  by given diagonal vector. This is basic Jacobi iteration; for tolerances,
273  iteration control, etc. wrap with SLISolver. */
275 {
276 public:
277  /** @brief Default constructor: the diagonal will be computed by subsequent
278  calls to SetOperator() using the Operator method AssembleDiagonal. */
279  /** In this case the array of essential tdofs will be empty. */
280  OperatorJacobiSmoother(const double damping=1.0);
281 
282  /** Setup a Jacobi smoother with the diagonal of @a a obtained by calling
283  a.AssembleDiagonal(). It is assumed that the underlying operator acts as
284  the identity on entries in ess_tdof_list, corresponding to (assembled)
285  DIAG_ONE policy or ConstrainedOperator in the matrix-free setting.
286 
287  @note For objects created with this constructor, calling SetOperator()
288  will only set the internal Operator pointer to the given new Operator
289  without any other changes to the object. This is done to preserve the
290  original behavior of this class. */
292  const Array<int> &ess_tdof_list,
293  const double damping=1.0);
294 
295  /** Application is by the *inverse* of the given vector. It is assumed that
296  the underlying operator acts as the identity on entries in ess_tdof_list,
297  corresponding to (assembled) DIAG_ONE policy or ConstrainedOperator in
298  the matrix-free setting.
299 
300  @note For objects created with this constructor, calling SetOperator()
301  will only set the internal Operator pointer to the given new Operator
302  without any other changes to the object. This is done to preserve the
303  original behavior of this class. */
305  const Array<int> &ess_tdof_list,
306  const double damping=1.0);
307 
309 
310  /// Replace diagonal entries with their absolute values.
311  void SetPositiveDiagonal(bool pos_diag = true) { use_abs_diag = pos_diag; }
312 
313  void Mult(const Vector &x, Vector &y) const;
314  void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); }
315 
316  /** @brief Recompute the diagonal using the method AssembleDiagonal of the
317  given new Operator, @a op. */
318  /** Note that (Par)BilinearForm operators are treated similar to the way they
319  are treated in the constructor that takes a BilinearForm parameter.
320  Specifically, this means that the OperatorJacobiSmoother will work with
321  true-dof vectors even though the size of the BilinearForm may be
322  different.
323 
324  When the new Operator, @a op, is not a (Par)BilinearForm, any previously
325  set array of essential true-dofs will be thrown away because in this case
326  any essential b.c. will be handled by the AssembleDiagonal method. */
327  void SetOperator(const Operator &op);
328 
329 private:
330  Vector dinv;
331  const double damping;
332  const Array<int> *ess_tdof_list; // not owned; may be NULL
333  mutable Vector residual;
334  /// Uses absolute values of the diagonal entries.
335  bool use_abs_diag = false;
336 
337  const Operator *oper; // not owned
338 
339  // To preserve the original behavior, some constructors set this flag to
340  // false to disallow updating the OperatorJacobiSmoother with SetOperator.
341  const bool allow_updates;
342 
343 public:
344  void Setup(const Vector &diag);
345 };
346 
347 /// Chebyshev accelerated smoothing with given vector, no matrix necessary
348 /** Potentially useful with tensorized operators, for example. This is just a
349  very basic Chebyshev iteration, if you want tolerances, iteration control,
350  etc. wrap this with SLISolver. */
352 {
353 public:
354  /** Application is by *inverse* of the given vector. It is assumed the
355  underlying operator acts as the identity on entries in ess_tdof_list,
356  corresponding to (assembled) DIAG_ONE policy or ConstrainedOperator in
357  the matrix-free setting. The estimated largest eigenvalue of the
358  diagonally preconditoned operator must be provided via
359  max_eig_estimate. */
360  OperatorChebyshevSmoother(const Operator &oper_, const Vector &d,
361  const Array<int>& ess_tdof_list,
362  int order, double max_eig_estimate);
363 
364  /// Deprecated: see pass-by-reference version above
365  MFEM_DEPRECATED
366  OperatorChebyshevSmoother(const Operator* oper_, const Vector &d,
367  const Array<int>& ess_tdof_list,
368  int order, double max_eig_estimate);
369 
370  /** Application is by *inverse* of the given vector. It is assumed the
371  underlying operator acts as the identity on entries in ess_tdof_list,
372  corresponding to (assembled) DIAG_ONE policy or ConstrainedOperator in
373  the matrix-free setting. The largest eigenvalue of the diagonally
374  preconditoned operator is estimated internally via a power method. The
375  accuracy of the estimated eigenvalue may be controlled via
376  power_iterations and power_tolerance. */
377 #ifdef MFEM_USE_MPI
378  OperatorChebyshevSmoother(const Operator &oper_, const Vector &d,
379  const Array<int>& ess_tdof_list,
380  int order, MPI_Comm comm = MPI_COMM_NULL,
381  int power_iterations = 10,
382  double power_tolerance = 1e-8);
383 
384  /// Deprecated: see pass-by-reference version above
385  MFEM_DEPRECATED
386  OperatorChebyshevSmoother(const Operator* oper_, const Vector &d,
387  const Array<int>& ess_tdof_list,
388  int order, MPI_Comm comm = MPI_COMM_NULL,
389  int power_iterations = 10,
390  double power_tolerance = 1e-8);
391 #else
392  OperatorChebyshevSmoother(const Operator &oper_, const Vector &d,
393  const Array<int>& ess_tdof_list,
394  int order, int power_iterations = 10,
395  double power_tolerance = 1e-8);
396 
397  /// Deprecated: see pass-by-reference version above
398  MFEM_DEPRECATED
399  OperatorChebyshevSmoother(const Operator* oper_, const Vector &d,
400  const Array<int>& ess_tdof_list,
401  int order, int power_iterations = 10,
402  double power_tolerance = 1e-8);
403 #endif
404 
406 
407  void Mult(const Vector&x, Vector &y) const;
408 
409  void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); }
410 
411  void SetOperator(const Operator &op_)
412  {
413  oper = &op_;
414  }
415 
416  void Setup();
417 
418 private:
419  const int order;
420  double max_eig_estimate;
421  const int N;
422  Vector dinv;
423  const Vector &diag;
424  Array<double> coeffs;
425  const Array<int>& ess_tdof_list;
426  mutable Vector residual;
427  mutable Vector helperVector;
428  const Operator* oper;
429 };
430 
431 
432 /// Stationary linear iteration: x <- x + B (b - A x)
434 {
435 protected:
436  mutable Vector r, z;
437 
438  void UpdateVectors();
439 
440 public:
441  SLISolver() { }
442 
443 #ifdef MFEM_USE_MPI
444  SLISolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
445 #endif
446 
447  virtual void SetOperator(const Operator &op)
449 
450  virtual void Mult(const Vector &b, Vector &x) const;
451 };
452 
453 /// Stationary linear iteration. (tolerances are squared)
454 void SLI(const Operator &A, const Vector &b, Vector &x,
455  int print_iter = 0, int max_num_iter = 1000,
456  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
457 
458 /// Preconditioned stationary linear iteration. (tolerances are squared)
459 void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
460  int print_iter = 0, int max_num_iter = 1000,
461  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
462 
463 
464 /// Conjugate gradient method
465 class CGSolver : public IterativeSolver
466 {
467 protected:
468  mutable Vector r, d, z;
469 
470  void UpdateVectors();
471 
472 public:
473  CGSolver() { }
474 
475 #ifdef MFEM_USE_MPI
476  CGSolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
477 #endif
478 
479  virtual void SetOperator(const Operator &op)
481 
482  virtual void Mult(const Vector &b, Vector &x) const;
483 };
484 
485 /// Conjugate gradient method. (tolerances are squared)
486 void CG(const Operator &A, const Vector &b, Vector &x,
487  int print_iter = 0, int max_num_iter = 1000,
488  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
489 
490 /// Preconditioned conjugate gradient method. (tolerances are squared)
491 void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x,
492  int print_iter = 0, int max_num_iter = 1000,
493  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
494 
495 
496 /// GMRES method
498 {
499 protected:
500  int m; // see SetKDim()
501 
502 public:
503  GMRESSolver() { m = 50; }
504 
505 #ifdef MFEM_USE_MPI
506  GMRESSolver(MPI_Comm comm_) : IterativeSolver(comm_) { m = 50; }
507 #endif
508 
509  /// Set the number of iteration to perform between restarts, default is 50.
510  void SetKDim(int dim) { m = dim; }
511 
512  virtual void Mult(const Vector &b, Vector &x) const;
513 };
514 
515 /// FGMRES method
517 {
518 protected:
519  int m;
520 
521 public:
522  FGMRESSolver() { m = 50; }
523 
524 #ifdef MFEM_USE_MPI
525  FGMRESSolver(MPI_Comm comm_) : IterativeSolver(comm_) { m = 50; }
526 #endif
527 
528  void SetKDim(int dim) { m = dim; }
529 
530  virtual void Mult(const Vector &b, Vector &x) const;
531 };
532 
533 /// GMRES method. (tolerances are squared)
534 int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M,
535  int &max_iter, int m, double &tol, double atol, int printit);
536 
537 /// GMRES method. (tolerances are squared)
538 void GMRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
539  int print_iter = 0, int max_num_iter = 1000, int m = 50,
540  double rtol = 1e-12, double atol = 1e-24);
541 
542 
543 /// BiCGSTAB method
545 {
546 protected:
547  mutable Vector p, phat, s, shat, t, v, r, rtilde;
548 
549  void UpdateVectors();
550 
551 public:
553 
554 #ifdef MFEM_USE_MPI
555  BiCGSTABSolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
556 #endif
557 
558  virtual void SetOperator(const Operator &op)
560 
561  virtual void Mult(const Vector &b, Vector &x) const;
562 };
563 
564 /// BiCGSTAB method. (tolerances are squared)
565 int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M,
566  int &max_iter, double &tol, double atol, int printit);
567 
568 /// BiCGSTAB method. (tolerances are squared)
569 void BiCGSTAB(const Operator &A, Solver &B, const Vector &b, Vector &x,
570  int print_iter = 0, int max_num_iter = 1000,
571  double rtol = 1e-12, double atol = 1e-24);
572 
573 
574 /// MINRES method
576 {
577 protected:
578  mutable Vector v0, v1, w0, w1, q;
579  mutable Vector u1; // used in the preconditioned version
580 
581 public:
583 
584 #ifdef MFEM_USE_MPI
585  MINRESSolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
586 #endif
587 
588  virtual void SetPreconditioner(Solver &pr)
589  {
591  if (oper) { u1.SetSize(width); }
592  }
593 
594  virtual void SetOperator(const Operator &op);
595 
596  virtual void Mult(const Vector &b, Vector &x) const;
597 };
598 
599 /// MINRES method without preconditioner. (tolerances are squared)
600 void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it = 0,
601  int max_it = 1000, double rtol = 1e-12, double atol = 1e-24);
602 
603 /// MINRES method with preconditioner. (tolerances are squared)
604 void MINRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
605  int print_it = 0, int max_it = 1000,
606  double rtol = 1e-12, double atol = 1e-24);
607 
608 
609 /// Newton's method for solving F(x)=b for a given operator F.
610 /** The method GetGradient() must be implemented for the operator F.
611  The preconditioner is used (in non-iterative mode) to evaluate
612  the action of the inverse gradient of the operator. */
614 {
615 protected:
616  mutable Vector r, c;
617  mutable Operator *grad;
618 
619  // Adaptive linear solver rtol variables
620 
621  // Method to determine rtol, 0 means the adaptive algorithm is deactivated.
622  int lin_rtol_type = 0;
623  // rtol to use in first iteration
624  double lin_rtol0;
625  // Maximum rtol
626  double lin_rtol_max;
627  // Function norm ||F(x)|| of the previous iterate
628  mutable double fnorm_last = 0.0;
629  // Linear residual norm of the previous iterate
630  mutable double lnorm_last = 0.0;
631  // Forcing term (linear residual rtol) from the previous iterate
632  mutable double eta_last = 0.0;
633  // Eisenstat-Walker factor gamma
634  double gamma;
635  // Eisenstat-Walker factor alpha
636  double alpha;
637 
638  /** @brief Method for the adaptive linear solver rtol invoked before the
639  linear solve. */
640  void AdaptiveLinRtolPreSolve(const Vector &x,
641  const int it,
642  const double fnorm) const;
643 
644  /** @brief Method for the adaptive linear solver rtol invoked after the
645  linear solve. */
646  void AdaptiveLinRtolPostSolve(const Vector &x,
647  const Vector &b,
648  const int it,
649  const double fnorm) const;
650 
651 public:
653 
654 #ifdef MFEM_USE_MPI
655  NewtonSolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
656 #endif
657  virtual void SetOperator(const Operator &op);
658 
659  /// Set the linear solver for inverting the Jacobian.
660  /** This method is equivalent to calling SetPreconditioner(). */
661  virtual void SetSolver(Solver &solver) { prec = &solver; }
662 
663  /// Solve the nonlinear system with right-hand side @a b.
664  /** If `b.Size() != Height()`, then @a b is assumed to be zero. */
665  virtual void Mult(const Vector &b, Vector &x) const;
666 
667  /** @brief This method can be overloaded in derived classes to implement line
668  search algorithms. */
669  /** The base class implementation (NewtonSolver) simply returns 1. A return
670  value of 0 indicates a failure, interrupting the Newton iteration. */
671  virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const
672  { return 1.0; }
673 
674  /** @brief This method can be overloaded in derived classes to perform
675  computations that need knowledge of the newest Newton state. */
676  virtual void ProcessNewState(const Vector &x) const { }
677 
678  /// Enable adaptive linear solver relative tolerance algorithm.
679  /** Compute a relative tolerance for the Krylov method after each nonlinear
680  iteration, based on the algorithm presented in [1].
681 
682  The maximum linear solver relative tolerance @a rtol_max should be < 1. For
683  @a type 1 the parameters @a alpha and @a gamma are ignored. For @a type 2
684  @a alpha has to be between 0 and 1 and @a gamma between 1 and 2.
685 
686  [1] Eisenstat, Stanley C., and Homer F. Walker. "Choosing the forcing terms
687  in an inexact Newton method."
688  */
689  void SetAdaptiveLinRtol(const int type = 2,
690  const double rtol0 = 0.5,
691  const double rtol_max = 0.9,
692  const double alpha = 0.5 * (1.0 + sqrt(5.0)),
693  const double gamma = 1.0);
694 };
695 
696 /** L-BFGS method for solving F(x)=b for a given operator F, by minimizing
697  the norm of F(x) - b. Requires only the action of the operator F. */
698 class LBFGSSolver : public NewtonSolver
699 {
700 protected:
701  int m = 10;
703 
705  {
706  for (int i = 0; i < skArray.Size(); i++)
707  {
708  delete skArray[i];
709  delete ykArray[i];
710  }
711  }
712 
714  {
716  skArray.SetSize(m);
717  ykArray.SetSize(m);
718  for (int i = 0; i < m; i++)
719  {
720  skArray[i] = new Vector(width);
721  ykArray[i] = new Vector(width);
722  skArray[i]->UseDevice(true);
723  ykArray[i]->UseDevice(true);
724  }
725  }
726 
727 public:
729 
730 #ifdef MFEM_USE_MPI
731  LBFGSSolver(MPI_Comm comm_) : NewtonSolver(comm_) { }
732 #endif
733 
734  virtual void SetOperator(const Operator &op)
735  {
738  }
739 
740  void SetHistorySize(int dim)
741  {
742  m = dim;
744  }
745 
746  /// Solve the nonlinear system with right-hand side @a b.
747  /** If `b.Size() != Height()`, then @a b is assumed to be zero. */
748  virtual void Mult(const Vector &b, Vector &x) const;
749 
750  virtual void SetPreconditioner(Solver &pr)
751  { MFEM_WARNING("L-BFGS won't use the given preconditioner."); }
752  virtual void SetSolver(Solver &solver)
753  { MFEM_WARNING("L-BFGS won't use the given solver."); }
754 
756 };
757 
758 
759 /** Adaptive restarted GMRES.
760  m_max and m_min(=1) are the maximal and minimal restart parameters.
761  m_step(=1) is the step to use for going from m_max and m_min.
762  cf(=0.4) is a desired convergence factor. */
763 int aGMRES(const Operator &A, Vector &x, const Vector &b,
764  const Operator &M, int &max_iter,
765  int m_max, int m_min, int m_step, double cf,
766  double &tol, double &atol, int printit);
767 
768 
769 /** Defines operators and constraints for the following optimization problem:
770  *
771  * Find x that minimizes the objective function F(x), subject to
772  * C(x) = c_e,
773  * d_lo <= D(x) <= d_hi,
774  * x_lo <= x <= x_hi.
775  *
776  * The operators F, C, D must take input of the same size (same width).
777  * Gradients of F, C, D might be needed, depending on the OptimizationSolver.
778  * When used with Hiop, gradients of C and D must be DenseMatrices.
779  * F always returns a scalar value, see CalcObjective(), CalcObjectiveGrad().
780  * C and D can have arbitrary heights.
781  * C and D can be NULL, meaning that their constraints are not used.
782  *
783  * When used in parallel, all Vectors are assumed to be true dof vectors, and
784  * the operators are expected to be defined for tdof vectors. */
786 {
787 protected:
788  /// Not owned, some can remain unused (NULL).
789  const Operator *C, *D;
790  const Vector *c_e, *d_lo, *d_hi, *x_lo, *x_hi;
791 
792 public:
793  const int input_size;
794 
795  /// In parallel, insize is the number of the local true dofs.
796  OptimizationProblem(int insize, const Operator *C_, const Operator *D_);
797 
798  /// Objective F(x). In parallel, the result should be reduced over tasks.
799  virtual double CalcObjective(const Vector &x) const = 0;
800  /// The result grad is expected to enter with the correct size.
801  virtual void CalcObjectiveGrad(const Vector &x, Vector &grad) const
802  { MFEM_ABORT("The objective gradient is not implemented."); }
803 
804  void SetEqualityConstraint(const Vector &c);
805  void SetInequalityConstraint(const Vector &dl, const Vector &dh);
806  void SetSolutionBounds(const Vector &xl, const Vector &xh);
807 
808  const Operator *GetC() const { return C; }
809  const Operator *GetD() const { return D; }
810  const Vector *GetEqualityVec() const { return c_e; }
811  const Vector *GetInequalityVec_Lo() const { return d_lo; }
812  const Vector *GetInequalityVec_Hi() const { return d_hi; }
813  const Vector *GetBoundsVec_Lo() const { return x_lo; }
814  const Vector *GetBoundsVec_Hi() const { return x_hi; }
815 
816  int GetNumConstraints() const;
817 };
818 
819 /// Abstract solver for OptimizationProblems.
821 {
822 protected:
824 
825 public:
827 #ifdef MFEM_USE_MPI
828  OptimizationSolver(MPI_Comm comm_): IterativeSolver(comm_), problem(NULL) { }
829 #endif
830  virtual ~OptimizationSolver() { }
831 
832  /** This function is virtual as solvers might need to perform some initial
833  * actions (e.g. validation) with the OptimizationProblem. */
835  { problem = &prob; }
836 
837  virtual void Mult(const Vector &xt, Vector &x) const = 0;
838 
839  virtual void SetPreconditioner(Solver &pr)
840  { MFEM_ABORT("Not meaningful for this solver."); }
841  virtual void SetOperator(const Operator &op)
842  { MFEM_ABORT("Not meaningful for this solver."); }
843 };
844 
845 /** SLBQP optimizer:
846  * (S)ingle (L)inearly Constrained with (B)ounds (Q)uadratic (P)rogram
847  *
848  * Minimize || x-x_t ||, subject to
849  * sum w_i x_i = a,
850  * x_lo <= x <= x_hi.
851  */
853 {
854 protected:
856  double a;
857 
858  /// Solve QP at fixed lambda
859  inline double solve(double l, const Vector &xt, Vector &x, int &nclip) const
860  {
861  add(xt, l, w, x);
862  if (problem == NULL) { x.median(lo,hi); }
863  else
864  {
867  }
868  nclip++;
869  if (problem == NULL) { return Dot(w, x) - a; }
870  else
871  {
872  Vector c(1);
873  // Includes parallel communication.
874  problem->GetC()->Mult(x, c);
875 
876  return c(0) - (*problem->GetEqualityVec())(0);
877  }
878  }
879 
880  inline void print_iteration(int it, double r, double l) const;
881 
882 public:
884 
885 #ifdef MFEM_USE_MPI
886  SLBQPOptimizer(MPI_Comm comm_) : OptimizationSolver(comm_) { }
887 #endif
888 
889  /** Setting an OptimizationProblem will overwrite the Vectors given by
890  * SetBounds and SetLinearConstraint. The objective function remains
891  * unchanged. */
892  virtual void SetOptimizationProblem(const OptimizationProblem &prob);
893 
894  void SetBounds(const Vector &lo_, const Vector &hi_);
895  void SetLinearConstraint(const Vector &w_, double a_);
896 
897  /** We let the target values play the role of the initial vector xt, from
898  * which the operator generates the optimal vector x. */
899  virtual void Mult(const Vector &xt, Vector &x) const;
900 };
901 
902 /** Block ILU solver:
903  * Performs a block ILU(k) approximate factorization with specified block
904  * size. Currently only k=0 is supported. This is useful as a preconditioner
905  * for DG-type discretizations, where the system matrix has a natural
906  * (elemental) block structure.
907  *
908  * In the case of DG discretizations, the block size should usually be set to
909  * either ndofs_per_element or vdim*ndofs_per_element (if the finite element
910  * space has Ordering::byVDIM). The block size must evenly divide the size of
911  * the matrix.
912  *
913  * Renumbering the blocks is also supported by specifying a reordering method.
914  * Currently greedy minimum discarded fill ordering and no reordering are
915  * supported. Renumbering the blocks can lead to a much better approximate
916  * factorization.
917  */
918 class BlockILU : public Solver
919 {
920 public:
921 
922  /// The reordering method used by the BlockILU factorization.
923  enum class Reordering
924  {
926  NONE
927  };
928 
929  /** Create an "empty" BlockILU solver. SetOperator must be called later to
930  * actually form the factorization
931  */
932  BlockILU(int block_size_,
934  int k_fill_ = 0);
935 
936  /** Create a block ILU approximate factorization for the matrix @a op.
937  * @a op should be of type either SparseMatrix or HypreParMatrix. In the
938  * case that @a op is a HypreParMatrix, the ILU factorization is performed
939  * on the diagonal blocks of the parallel decomposition.
940  */
941  BlockILU(const Operator &op, int block_size_ = 1,
943  int k_fill_ = 0);
944 
945  /** Perform the block ILU factorization for the matrix @a op.
946  * As in the constructor, @a op must either be a SparseMatrix or
947  * HypreParMatrix
948  */
949  void SetOperator(const Operator &op);
950 
951  /// Solve the system `LUx = b`, where `L` and `U` are the block ILU factors.
952  void Mult(const Vector &b, Vector &x) const;
953 
954  /** Get the I array for the block CSR representation of the factorization.
955  * Similar to SparseMatrix::GetI(). Mostly used for testing.
956  */
957  int *GetBlockI() { return IB.GetData(); }
958 
959  /** Get the J array for the block CSR representation of the factorization.
960  * Similar to SparseMatrix::GetJ(). Mostly used for testing.
961  */
962  int *GetBlockJ() { return JB.GetData(); }
963 
964  /** Get the data array for the block CSR representation of the factorization.
965  * Similar to SparseMatrix::GetData(). Mostly used for testing.
966  */
967  double *GetBlockData() { return AB.Data(); }
968 
969 private:
970  /// Set up the block CSR structure corresponding to a sparse matrix @a A
971  void CreateBlockPattern(const class SparseMatrix &A);
972 
973  /// Perform the block ILU factorization
974  void Factorize();
975 
976  int block_size;
977 
978  /// Fill level for block ILU(k) factorizations. Only k=0 is supported.
979  int k_fill;
980 
981  Reordering reordering;
982 
983  /// Temporary vector used in the Mult() function.
984  mutable Vector y;
985 
986  /// Permutation and inverse permutation vectors for the block reordering.
987  Array<int> P, Pinv;
988 
989  /** Block CSR storage of the factorization. The block upper triangular part
990  * stores the U factor. The L factor implicitly has identity on the diagonal
991  * blocks, and the rest of L is given by the strictly block lower triangular
992  * part.
993  */
994  Array<int> IB, ID, JB;
995  DenseTensor AB;
996 
997  /// DB(i) stores the LU factorization of the i'th diagonal block
998  mutable DenseTensor DB;
999  /// Pivot arrays for the LU factorizations given by #DB
1000  mutable Array<int> ipiv;
1001 };
1002 
1003 
1004 /// Monitor that checks whether the residual is zero at a given set of dofs.
1005 /** This monitor is useful for checking if the initial guess, rhs, operator, and
1006  preconditioner are properly setup for solving in the subspace with imposed
1007  essential boundary conditions. */
1009 {
1010 protected:
1011  const Array<int> *ess_dofs_list; ///< Not owned
1012 
1013 public:
1014  ResidualBCMonitor(const Array<int> &ess_dofs_list_)
1015  : ess_dofs_list(&ess_dofs_list_) { }
1016 
1017  void MonitorResidual(int it, double norm, const Vector &r,
1018  bool final) override;
1019 };
1020 
1021 
1022 #ifdef MFEM_USE_SUITESPARSE
1023 
1024 /// Direct sparse solver using UMFPACK
1025 class UMFPackSolver : public Solver
1026 {
1027 protected:
1030  void *Numeric;
1031  SuiteSparse_long *AI, *AJ;
1032 
1033  void Init();
1034 
1035 public:
1036  double Control[UMFPACK_CONTROL];
1037  mutable double Info[UMFPACK_INFO];
1038 
1039  /** @brief For larger matrices, if the solver fails, set the parameter @a
1040  use_long_ints_ = true. */
1041  UMFPackSolver(bool use_long_ints_ = false)
1042  : use_long_ints(use_long_ints_) { Init(); }
1043  /** @brief Factorize the given SparseMatrix using the defaults. For larger
1044  matrices, if the solver fails, set the parameter @a use_long_ints_ =
1045  true. */
1046  UMFPackSolver(SparseMatrix &A, bool use_long_ints_ = false)
1047  : use_long_ints(use_long_ints_) { Init(); SetOperator(A); }
1048 
1049  /** @brief Factorize the given Operator @a op which must be a SparseMatrix.
1050 
1051  The factorization uses the parameters set in the #Control data member.
1052  @note This method calls SparseMatrix::SortColumnIndices() with @a op,
1053  modifying the matrix if the column indices are not already sorted. */
1054  virtual void SetOperator(const Operator &op);
1055 
1056  /// Set the print level field in the #Control data member.
1057  void SetPrintLevel(int print_lvl) { Control[UMFPACK_PRL] = print_lvl; }
1058 
1059  virtual void Mult(const Vector &b, Vector &x) const;
1060  virtual void MultTranspose(const Vector &b, Vector &x) const;
1061 
1062  virtual ~UMFPackSolver();
1063 };
1064 
1065 /// Direct sparse solver using KLU
1066 class KLUSolver : public Solver
1067 {
1068 protected:
1070  klu_symbolic *Symbolic;
1071  klu_numeric *Numeric;
1072 
1073  void Init();
1074 
1075 public:
1077  : mat(0),Symbolic(0),Numeric(0)
1078  { Init(); }
1080  : mat(0),Symbolic(0),Numeric(0)
1081  { Init(); SetOperator(A); }
1082 
1083  // Works on sparse matrices only; calls SparseMatrix::SortColumnIndices().
1084  virtual void SetOperator(const Operator &op);
1085 
1086  virtual void Mult(const Vector &b, Vector &x) const;
1087  virtual void MultTranspose(const Vector &b, Vector &x) const;
1088 
1089  virtual ~KLUSolver();
1090 
1091  mutable klu_common Common;
1092 };
1093 
1094 #endif // MFEM_USE_SUITESPARSE
1095 
1096 /// Block diagonal solver for A, each block is inverted by direct solver
1098 {
1099  SparseMatrix& block_dof;
1100  mutable Array<int> local_dofs;
1101  mutable Vector sub_rhs;
1102  mutable Vector sub_sol;
1103  std::unique_ptr<DenseMatrixInverse[]> block_solvers;
1104 public:
1105  /// block_dof is a boolean matrix, block_dof(i, j) = 1 if j-th dof belongs to
1106  /// i-th block, block_dof(i, j) = 0 otherwise.
1107  DirectSubBlockSolver(const SparseMatrix& A, const SparseMatrix& block_dof);
1108  virtual void Mult(const Vector &x, Vector &y) const;
1109  virtual void SetOperator(const Operator &op) { }
1110 };
1111 
1112 /// Solver S such that I - A * S = (I - A * S1) * (I - A * S0).
1113 /// That is, S = S0 + S1 - S1 * A * S0.
1114 class ProductSolver : public Solver
1115 {
1116  OperatorPtr A;
1117  OperatorPtr S0;
1118  OperatorPtr S1;
1119 public:
1121  bool ownA, bool ownS0, bool ownS1)
1122  : Solver(A_->NumRows()), A(A_, ownA), S0(S0_, ownS0), S1(S1_, ownS1) { }
1123  virtual void Mult(const Vector &x, Vector &y) const;
1124  virtual void MultTranspose(const Vector &x, Vector &y) const;
1125  virtual void SetOperator(const Operator &op) { }
1126 };
1127 
1128 /// Solver wrapper which orthogonalizes the input and output vector
1129 /**
1130  * OrthoSolver wraps an existing Solver and orthogonalizes the input vector
1131  * before passing it to the Mult() method of the Solver. This is a convenience
1132  * implementation to handle e.g. a Poisson problem with pure Neumann boundary
1133  * conditions, where this procedure removes the Nullspace.
1134  */
1135 class OrthoSolver : public Solver
1136 {
1137 private:
1138 #ifdef MFEM_USE_MPI
1139  MPI_Comm mycomm;
1140  mutable HYPRE_BigInt global_size;
1141  const bool parallel;
1142 #else
1143  mutable int global_size;
1144 #endif
1145 
1146 public:
1147  OrthoSolver();
1148 #ifdef MFEM_USE_MPI
1149  OrthoSolver(MPI_Comm mycomm_);
1150 #endif
1151 
1152  /// Set the solver used by the OrthoSolver.
1153  /** The action of the OrthoSolver is given by P * s * P where P is the
1154  projection to the subspace of vectors with zero sum. Calling this method
1155  is required before calling SetOperator() or Mult(). */
1156  void SetSolver(Solver &s);
1157 
1158  /// Set the Operator that is the OrthoSolver is to invert (approximately).
1159  /** The Operator @a op is simply forwarded to the solver object given by
1160  SetSolver() which needs to be called before this method. Calling this
1161  method is optional when the solver already has an associated Operator. */
1162  virtual void SetOperator(const Operator &op);
1163 
1164  /** @brief Perform the action of the OrthoSolver: P * solver * P where P is
1165  the projection to the subspace of vectors with zero sum. */
1166  /** @note The projection P can be written as P = I - 1 1^T / (1^T 1) where
1167  I is the identity matrix and 1 is the column-vector with all components
1168  equal to 1. */
1169  void Mult(const Vector &b, Vector &x) const;
1170 
1171 private:
1172  Solver *solver = nullptr;
1173 
1174  mutable Vector b_ortho;
1175 
1176  void Orthogonalize(const Vector &v, Vector &v_ortho) const;
1177 };
1178 
1179 #ifdef MFEM_USE_MPI
1180 /** This smoother does relaxations on an auxiliary space (determined by a map
1181  from the original space to the auxiliary space provided by the user).
1182  The smoother on the auxiliary space is a HypreSmoother. Its options can be
1183  modified through GetSmoother.
1184  For example, the space can be the nullspace of div/curl, in which case the
1185  smoother can be used to construct a Hiptmair smoother. */
1186 class AuxSpaceSmoother : public Solver
1187 {
1188  OperatorPtr aux_map_;
1189  OperatorPtr aux_system_;
1190  OperatorPtr aux_smoother_;
1191  void Mult(const Vector &x, Vector &y, bool transpose) const;
1192 public:
1193  AuxSpaceSmoother(const HypreParMatrix &op, HypreParMatrix *aux_map,
1194  bool op_is_symmetric = true, bool own_aux_map = false);
1195  virtual void Mult(const Vector &x, Vector &y) const { Mult(x, y, false); }
1196  virtual void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y, true); }
1197  virtual void SetOperator(const Operator &op) { }
1198  HypreSmoother& GetSmoother() { return *aux_smoother_.As<HypreSmoother>(); }
1199 };
1200 #endif // MFEM_USE_MPI
1201 
1202 }
1203 
1204 #endif // MFEM_SOLVERS
const Array< int > * ess_dofs_list
Not owned.
Definition: solvers.hpp:1011
MINRESSolver(MPI_Comm comm_)
Definition: solvers.hpp:585
const Vector * GetInequalityVec_Hi() const
Definition: solvers.hpp:812
GMRESSolver(MPI_Comm comm_)
Definition: solvers.hpp:506
Conjugate gradient method.
Definition: solvers.hpp:465
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
CGSolver(MPI_Comm comm_)
Definition: solvers.hpp:476
HypreSmoother & GetSmoother()
Definition: solvers.hpp:1198
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:3222
double Info[UMFPACK_INFO]
Definition: solvers.hpp:1037
void Mult(const Vector &b, Vector &x) const
Solve the system LUx = b, where L and U are the block ILU factors.
Definition: solvers.cpp:2991
Chebyshev accelerated smoothing with given vector, no matrix necessary.
Definition: solvers.hpp:351
SparseMatrix * mat
Definition: solvers.hpp:1029
int GetNumIterations() const
Definition: solvers.hpp:246
const Operator * oper
Definition: solvers.hpp:120
MPI_Comm GetComm() const
Return the associated MPI communicator, or MPI_COMM_NULL if no communicator is set.
Definition: solvers.hpp:264
Array< Vector * > ykArray
Definition: solvers.hpp:702
void SetPositiveDiagonal(bool pos_diag=true)
Replace diagonal entries with their absolute values.
Definition: solvers.hpp:311
SuiteSparse_long * AI
Definition: solvers.hpp:1031
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:734
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:711
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:3363
void SetHistorySize(int dim)
Definition: solvers.hpp:740
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:841
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
void Monitor(int it, double norm, const Vector &r, const Vector &x, bool final=false) const
Definition: solvers.cpp:190
virtual void ProcessNewState(const Vector &x) const
This method can be overloaded in derived classes to perform computations that need knowledge of the n...
Definition: solvers.hpp:676
bool summary
A summary of the solver process will be reported after the last iteration to mfem::out.
Definition: solvers.hpp:91
bool first_and_last
Information about the first and last iteration will be printed to mfem::out.
Definition: solvers.hpp:94
OperatorJacobiSmoother(const double damping=1.0)
Default constructor: the diagonal will be computed by subsequent calls to SetOperator() using the Ope...
Definition: solvers.cpp:200
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
class IterativeSolver * iter_solver
The last IterativeSolver to which this monitor was attached.
Definition: solvers.hpp:40
T * GetData()
Returns the data.
Definition: array.hpp:112
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:300
FGMRES method.
Definition: solvers.hpp:516
int * GetBlockI()
Definition: solvers.hpp:957
void InitializeStorageVectors()
Definition: solvers.hpp:713
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1159
Direct sparse solver using KLU.
Definition: solvers.hpp:1066
bool errors
If a fatal problem has been detected the failure will be reported to mfem::err.
Definition: solvers.hpp:82
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: solvers.hpp:661
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:3324
void SetOperator(const Operator &op)
Definition: solvers.cpp:2760
MINRES method.
Definition: solvers.hpp:575
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:1046
const Operator * GetC() const
Definition: solvers.hpp:808
bool GetConverged() const
Definition: solvers.hpp:247
Reordering
The reordering method used by the BlockILU factorization.
Definition: solvers.hpp:923
virtual ~IterativeSolverMonitor()
Definition: solvers.hpp:45
ResidualBCMonitor(const Array< int > &ess_dofs_list_)
Definition: solvers.hpp:1014
int print_level
(DEPRECATED) Legacy print level definition, which is left for compatibility with custom iterative sol...
Definition: solvers.hpp:130
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.hpp:1125
PrintLevel()=default
Initializes the print level to suppress.
virtual double CalcObjective(const Vector &x) const =0
Objective F(x). In parallel, the result should be reduced over tasks.
SLISolver(MPI_Comm comm_)
Definition: solvers.hpp:444
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:313
virtual ~OptimizationSolver()
Definition: solvers.hpp:830
void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: solvers.hpp:409
double Norm(const Vector &x) const
Definition: solvers.hpp:171
const Vector * d_hi
Definition: solvers.hpp:790
virtual void SetOperator(const Operator &op)
Set the Operator that is the OrthoSolver is to invert (approximately).
Definition: solvers.cpp:3433
double GetFinalNorm() const
Definition: solvers.hpp:248
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:1565
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:525
void SetKDim(int dim)
Definition: solvers.hpp:528
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:179
const Vector * x_lo
Definition: solvers.hpp:790
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:1025
Operator * grad
Definition: solvers.hpp:617
const Vector * x_hi
Definition: solvers.hpp:790
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:971
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1609
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:588
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
const Vector * GetBoundsVec_Lo() const
Definition: solvers.hpp:813
double solve(double l, const Vector &xt, Vector &x, int &nclip) const
Solve QP at fixed lambda.
Definition: solvers.hpp:859
int GetNumConstraints() const
Definition: solvers.cpp:2335
Data type sparse matrix.
Definition: sparsemat.hpp:50
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:1775
virtual ~UMFPackSolver()
Definition: solvers.cpp:3257
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:274
double b
Definition: lissajous.cpp:42
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
OperatorChebyshevSmoother(const Operator &oper_, const Vector &d, const Array< int > &ess_tdof_list, int order, double max_eig_estimate)
Definition: solvers.cpp:328
void SetLinearConstraint(const Vector &w_, double a_)
Definition: solvers.cpp:2362
void median(const Vector &lo, const Vector &hi)
v = median(v,lo,hi) entrywise. Implementation assumes lo &lt;= hi.
Definition: vector.cpp:523
DirectSubBlockSolver(const SparseMatrix &A, const SparseMatrix &block_dof)
Definition: solvers.cpp:3348
SuiteSparse_long * AJ
Definition: solvers.hpp:1031
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition: solvers.hpp:510
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:613
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:889
void print_iteration(int it, double r, double l) const
Definition: solvers.cpp:2368
void Setup(const Vector &diag)
Definition: solvers.cpp:277
Parallel smoothers in hypre.
Definition: hypre.hpp:962
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:904
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:3279
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1588
double Dot(const Vector &x, const Vector &y) const
Definition: solvers.cpp:55
Array< Vector * > skArray
Definition: solvers.hpp:702
prob_type prob
Definition: ex25.cpp:153
klu_symbolic * Symbolic
Definition: solvers.hpp:1070
void SetIterativeSolver(const IterativeSolver &solver)
This method is invoked by IterativeSolver::SetMonitor, informing the monitor which IterativeSolver is...
Definition: solvers.hpp:61
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.hpp:1197
const Vector * GetBoundsVec_Hi() const
Definition: solvers.hpp:814
BlockILU(int block_size_, Reordering reordering_=Reordering::MINIMUM_DISCARDED_FILL, int k_fill_=0)
Definition: solvers.cpp:2742
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1814
void SetPrintLevel(int print_lvl)
Set the print level field in the Control data member.
Definition: solvers.hpp:1057
void AdaptiveLinRtolPostSolve(const Vector &x, const Vector &b, const int it, const double fnorm) const
Method for the adaptive linear solver rtol invoked after the linear solve.
Definition: solvers.cpp:1982
Stationary linear iteration: x &lt;- x + B (b - A x)
Definition: solvers.hpp:433
void SetInequalityConstraint(const Vector &dl, const Vector &dh)
Definition: solvers.cpp:2317
void SetEqualityConstraint(const Vector &c)
Definition: solvers.cpp:2309
Abstract base class for an iterative solver monitor.
Definition: solvers.hpp:36
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:839
void SetOperator(const Operator &op_)
Set/update the solver for the given operator.
Definition: solvers.hpp:411
ProductSolver(Operator *A_, Solver *S0_, Solver *S1_, bool ownA, bool ownS0, bool ownS1)
Definition: solvers.hpp:1120
void SetAbsTol(double atol)
Definition: solvers.hpp:199
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: solvers.cpp:3395
PrintLevel print_options
Output behavior for the iterative solver.
Definition: solvers.hpp:137
void SetRelTol(double rtol)
Definition: solvers.hpp:198
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:1036
Monitor that checks whether the residual is zero at a given set of dofs.
Definition: solvers.hpp:1008
Abstract base class for iterative solver.
Definition: solvers.hpp:66
UMFPackSolver(bool use_long_ints_=false)
For larger matrices, if the solver fails, set the parameter use_long_ints_ = true.
Definition: solvers.hpp:1041
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:2002
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1373
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: operator.hpp:70
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3310
void DeleteStorageVectors()
Definition: solvers.hpp:704
void SetAdaptiveLinRtol(const int type=2, const double rtol0=0.5, const double rtol_max=0.9, const double alpha=0.5 *(1.0+sqrt(5.0)), const double gamma=1.0)
Enable adaptive linear solver relative tolerance algorithm.
Definition: solvers.cpp:1922
int max_iter
Limit for the number of iterations the solver is allowed to do.
Definition: solvers.hpp:151
LBFGSSolver(MPI_Comm comm_)
Definition: solvers.hpp:731
HYPRE_Int HYPRE_BigInt
double * GetBlockData()
Definition: solvers.hpp:967
GMRES method.
Definition: solvers.hpp:497
void Mult(const Vector &b, Vector &x) const
Perform the action of the OrthoSolver: P * solver * P where P is the projection to the subspace of ve...
Definition: solvers.cpp:3443
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: solvers.hpp:752
static int GuessLegacyPrintLevel(PrintLevel)
Use some heuristics to guess a legacy print level corresponding to the given PrintLevel.
Definition: solvers.cpp:149
const Vector * GetEqualityVec() const
Definition: solvers.hpp:810
FGMRESSolver(MPI_Comm comm_)
Definition: solvers.hpp:525
virtual void SetOptimizationProblem(const OptimizationProblem &prob)
Definition: solvers.cpp:2343
void SetMonitor(IterativeSolverMonitor &m)
Set the iterative solver monitor.
Definition: solvers.hpp:258
const Operator * GetD() const
Definition: solvers.hpp:809
double a
Definition: lissajous.cpp:41
void UpdateVectors()
Definition: solvers.cpp:519
Abstract solver for OptimizationProblems.
Definition: solvers.hpp:820
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
Solver wrapper which orthogonalizes the input and output vector.
Definition: solvers.hpp:1135
virtual ~LBFGSSolver()
Definition: solvers.hpp:755
klu_common Common
Definition: solvers.hpp:1091
double rel_tol
Relative tolerance.
Definition: solvers.hpp:154
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:671
OptimizationProblem(int insize, const Operator *C_, const Operator *D_)
In parallel, insize is the number of the local true dofs.
Definition: solvers.cpp:2299
void SetSolutionBounds(const Vector &xl, const Vector &xh)
Definition: solvers.cpp:2327
const Vector * GetInequalityVec_Lo() const
Definition: solvers.hpp:811
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.hpp:1195
AuxSpaceSmoother(const HypreParMatrix &op, HypreParMatrix *aux_map, bool op_is_symmetric=true, bool own_aux_map=false)
Definition: solvers.cpp:3500
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:479
virtual void Mult(const Vector &xt, Vector &x) const =0
Operator application: y=A(x).
BiCGSTABSolver(MPI_Comm comm_)
Definition: solvers.hpp:555
void SetOperator(const Operator &op)
Recompute the diagonal using the method AssembleDiagonal of the given new Operator, op.
Definition: solvers.cpp:241
int dim
Definition: ex24.cpp:53
SLBQPOptimizer(MPI_Comm comm_)
Definition: solvers.hpp:886
void AdaptiveLinRtolPreSolve(const Vector &x, const int it, const double fnorm) const
Method for the adaptive linear solver rtol invoked before the linear solve.
Definition: solvers.cpp:1935
PrintLevel FromLegacyPrintLevel(int)
Convert a legacy print level integer to a PrintLevel object.
Definition: solvers.cpp:114
int * GetBlockJ()
Definition: solvers.hpp:962
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:558
virtual void Mult(const Vector &xt, Vector &x) const
Definition: solvers.cpp:2377
void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: solvers.hpp:314
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:479
virtual void CalcObjectiveGrad(const Vector &x, Vector &grad) const
The result grad is expected to enter with the correct size.
Definition: solvers.hpp:801
Vector data type.
Definition: vector.hpp:60
void SetSolver(Solver &s)
Set the solver used by the OrthoSolver.
Definition: solvers.cpp:3424
void MonitorResidual(int it, double norm, const Vector &r, bool final) override
Monitor the residual vector r.
Definition: solvers.cpp:3042
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: solvers.hpp:1196
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:750
const Vector * d_lo
Definition: solvers.hpp:790
virtual void MonitorSolution(int it, double norm, const Vector &x, bool final)
Monitor the solution vector x.
Definition: solvers.hpp:54
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
KLUSolver(SparseMatrix &A)
Definition: solvers.hpp:1079
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:2144
virtual ~KLUSolver()
Definition: solvers.cpp:3338
RefCoord s[3]
double * Data()
Definition: densemat.hpp:1070
void UpdateVectors()
Definition: solvers.cpp:702
BiCGSTAB method.
Definition: solvers.hpp:544
bool warnings
If a non-fatal problem has been detected some context-specific information will be reported to mfem::...
Definition: solvers.hpp:85
bool iterations
Detailed information about each iteration will be reported to mfem::out.
Definition: solvers.hpp:88
OptimizationSolver(MPI_Comm comm_)
Definition: solvers.hpp:828
SparseMatrix * mat
Definition: solvers.hpp:1069
Base class for solvers.
Definition: operator.hpp:655
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.hpp:1109
const OptimizationProblem * problem
Definition: solvers.hpp:823
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Block diagonal solver for A, each block is inverted by direct solver.
Definition: solvers.hpp:1097
void SetBounds(const Vector &lo_, const Vector &hi_)
Definition: solvers.cpp:2356
const Operator * D
Definition: solvers.hpp:789
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:953
virtual void SetOptimizationProblem(const OptimizationProblem &prob)
Definition: solvers.hpp:834
NewtonSolver(MPI_Comm comm_)
Definition: solvers.hpp:655
const Vector * c_e
Definition: solvers.hpp:790
IterativeSolverMonitor * monitor
Definition: solvers.hpp:122
klu_numeric * Numeric
Definition: solvers.hpp:1071
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:3378
virtual void MonitorResidual(int it, double norm, const Vector &r, bool final)
Monitor the residual vector r.
Definition: solvers.hpp:48
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3189
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:447
Settings for the output behavior of the IterativeSolver.
Definition: solvers.hpp:78
const Operator * C
Not owned, some can remain unused (NULL).
Definition: solvers.hpp:789
double abs_tol
Absolute tolerance.
Definition: solvers.hpp:157
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:670
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:1335
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1803
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:3094