MFEM  v4.3.0
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-2021, 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 ItertiveSolver::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 #ifdef MFEM_USE_MPI
69 private:
70  int dot_prod_type; // 0 - local, 1 - global over 'comm'
71  MPI_Comm comm;
72 #endif
73 
74 protected:
75  const Operator *oper;
78 
80  double rel_tol, abs_tol;
81 
82  // stats
83  mutable int final_iter, converged;
84  mutable double final_norm;
85 
86  double Dot(const Vector &x, const Vector &y) const;
87  double Norm(const Vector &x) const { return sqrt(Dot(x, x)); }
88  void Monitor(int it, double norm, const Vector& r, const Vector& x,
89  bool final=false) const;
90 
91 public:
93 
94 #ifdef MFEM_USE_MPI
95  IterativeSolver(MPI_Comm comm_);
96 #endif
97 
98  void SetRelTol(double rtol) { rel_tol = rtol; }
99  void SetAbsTol(double atol) { abs_tol = atol; }
100  void SetMaxIter(int max_it) { max_iter = max_it; }
101  void SetPrintLevel(int print_lvl);
102 
103  int GetNumIterations() const { return final_iter; }
104  int GetConverged() const { return converged; }
105  double GetFinalNorm() const { return final_norm; }
106 
107  /// This should be called before SetOperator
108  virtual void SetPreconditioner(Solver &pr);
109 
110  /// Also calls SetOperator for the preconditioner if there is one
111  virtual void SetOperator(const Operator &op);
112 
113  /// Set the iterative solver monitor
115  { monitor = &m; m.SetIterativeSolver(*this); }
116 
117 #ifdef MFEM_USE_MPI
118  /** @brief Return the associated MPI communicator, or MPI_COMM_NULL if no
119  communicator is set. */
120  MPI_Comm GetComm() const
121  { return dot_prod_type == 0 ? MPI_COMM_NULL : comm; }
122 #endif
123 };
124 
125 
126 /// Jacobi smoothing for a given bilinear form (no matrix necessary).
127 /** Useful with tensorized, partially assembled operators. Can also be defined
128  by given diagonal vector. This is basic Jacobi iteration; for tolerances,
129  iteration control, etc. wrap with SLISolver. */
131 {
132 public:
133  /** @brief Default constructor: the diagonal will be computed by subsequent
134  calls to SetOperator() using the Operator method AssembleDiagonal. */
135  /** In this case the array of essential tdofs will be empty. */
136  OperatorJacobiSmoother(const double damping=1.0);
137 
138  /** Setup a Jacobi smoother with the diagonal of @a a obtained by calling
139  a.AssembleDiagonal(). It is assumed that the underlying operator acts as
140  the identity on entries in ess_tdof_list, corresponding to (assembled)
141  DIAG_ONE policy or ConstrainedOperator in the matrix-free setting.
142 
143  @note For objects created with this constructor, calling SetOperator()
144  will only set the internal Operator pointer to the given new Operator
145  without any other changes to the object. This is done to preserve the
146  original behavior of this class. */
148  const Array<int> &ess_tdof_list,
149  const double damping=1.0);
150 
151  /** Application is by the *inverse* of the given vector. It is assumed that
152  the underlying operator acts as the identity on entries in ess_tdof_list,
153  corresponding to (assembled) DIAG_ONE policy or ConstrainedOperator in
154  the matrix-free setting.
155 
156  @note For objects created with this constructor, calling SetOperator()
157  will only set the internal Operator pointer to the given new Operator
158  without any other changes to the object. This is done to preserve the
159  original behavior of this class. */
161  const Array<int> &ess_tdof_list,
162  const double damping=1.0);
163 
165 
166  void Mult(const Vector &x, Vector &y) const;
167  void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); }
168 
169  /** @brief Recompute the diagonal using the method AssembleDiagonal of the
170  given new Operator, @a op. */
171  /** Note that (Par)BilinearForm operators are treated similar to the way they
172  are treated in the constructor that takes a BilinearForm parameter.
173  Specifically, this means that the OperatorJacobiSmoother will work with
174  true-dof vectors even though the size of the BilinearForm may be
175  different.
176 
177  When the new Operator, @a op, is not a (Par)BilinearForm, any previously
178  set array of essential true-dofs will be thrown away because in this case
179  any essential b.c. will be handled by the AssembleDiagonal method. */
180  void SetOperator(const Operator &op);
181 
182 private:
183  Vector dinv;
184  const double damping;
185  const Array<int> *ess_tdof_list; // not owned; may be NULL
186  mutable Vector residual;
187 
188  const Operator *oper; // not owned
189 
190  // To preserve the original behavior, some constructors set this flag to
191  // false to disallow updating the OperatorJacobiSmoother with SetOperator.
192  const bool allow_updates;
193 
194 public:
195  void Setup(const Vector &diag);
196 };
197 
198 /// Chebyshev accelerated smoothing with given vector, no matrix necessary
199 /** Potentially useful with tensorized operators, for example. This is just a
200  very basic Chebyshev iteration, if you want tolerances, iteration control,
201  etc. wrap this with SLISolver. */
203 {
204 public:
205  /** Application is by *inverse* of the given vector. It is assumed the
206  underlying operator acts as the identity on entries in ess_tdof_list,
207  corresponding to (assembled) DIAG_ONE policy or ConstrainedOperator in
208  the matrix-free setting. The estimated largest eigenvalue of the
209  diagonally preconditoned operator must be provided via
210  max_eig_estimate. */
211  OperatorChebyshevSmoother(const Operator &oper_, const Vector &d,
212  const Array<int>& ess_tdof_list,
213  int order, double max_eig_estimate);
214 
215  /// Deprecated: see pass-by-reference version above
216  MFEM_DEPRECATED
217  OperatorChebyshevSmoother(const Operator* oper_, const Vector &d,
218  const Array<int>& ess_tdof_list,
219  int order, double max_eig_estimate);
220 
221  /** Application is by *inverse* of the given vector. It is assumed the
222  underlying operator acts as the identity on entries in ess_tdof_list,
223  corresponding to (assembled) DIAG_ONE policy or ConstrainedOperator in
224  the matrix-free setting. The largest eigenvalue of the diagonally
225  preconditoned operator is estimated internally via a power method. The
226  accuracy of the estimated eigenvalue may be controlled via
227  power_iterations and power_tolerance. */
228 #ifdef MFEM_USE_MPI
229  OperatorChebyshevSmoother(const Operator &oper_, const Vector &d,
230  const Array<int>& ess_tdof_list,
231  int order, MPI_Comm comm = MPI_COMM_NULL,
232  int power_iterations = 10,
233  double power_tolerance = 1e-8);
234 
235  /// Deprecated: see pass-by-reference version above
236  MFEM_DEPRECATED
237  OperatorChebyshevSmoother(const Operator* oper_, const Vector &d,
238  const Array<int>& ess_tdof_list,
239  int order, MPI_Comm comm = MPI_COMM_NULL,
240  int power_iterations = 10,
241  double power_tolerance = 1e-8);
242 #else
243  OperatorChebyshevSmoother(const Operator &oper_, const Vector &d,
244  const Array<int>& ess_tdof_list,
245  int order, int power_iterations = 10,
246  double power_tolerance = 1e-8);
247 
248  /// Deprecated: see pass-by-reference version above
249  MFEM_DEPRECATED
250  OperatorChebyshevSmoother(const Operator* oper_, const Vector &d,
251  const Array<int>& ess_tdof_list,
252  int order, int power_iterations = 10,
253  double power_tolerance = 1e-8);
254 #endif
255 
257 
258  void Mult(const Vector&x, Vector &y) const;
259 
260  void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); }
261 
262  void SetOperator(const Operator &op_)
263  {
264  oper = &op_;
265  }
266 
267  void Setup();
268 
269 private:
270  const int order;
271  double max_eig_estimate;
272  const int N;
273  Vector dinv;
274  const Vector &diag;
275  Array<double> coeffs;
276  const Array<int>& ess_tdof_list;
277  mutable Vector residual;
278  mutable Vector helperVector;
279  const Operator* oper;
280 };
281 
282 
283 /// Stationary linear iteration: x <- x + B (b - A x)
285 {
286 protected:
287  mutable Vector r, z;
288 
289  void UpdateVectors();
290 
291 public:
292  SLISolver() { }
293 
294 #ifdef MFEM_USE_MPI
295  SLISolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
296 #endif
297 
298  virtual void SetOperator(const Operator &op)
300 
301  virtual void Mult(const Vector &b, Vector &x) const;
302 };
303 
304 /// Stationary linear iteration. (tolerances are squared)
305 void SLI(const Operator &A, const Vector &b, Vector &x,
306  int print_iter = 0, int max_num_iter = 1000,
307  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
308 
309 /// Preconditioned stationary linear iteration. (tolerances are squared)
310 void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
311  int print_iter = 0, int max_num_iter = 1000,
312  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
313 
314 
315 /// Conjugate gradient method
316 class CGSolver : public IterativeSolver
317 {
318 protected:
319  mutable Vector r, d, z;
320 
321  void UpdateVectors();
322 
323 public:
324  CGSolver() { }
325 
326 #ifdef MFEM_USE_MPI
327  CGSolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
328 #endif
329 
330  virtual void SetOperator(const Operator &op)
332 
333  virtual void Mult(const Vector &b, Vector &x) const;
334 };
335 
336 /// Conjugate gradient method. (tolerances are squared)
337 void CG(const Operator &A, const Vector &b, Vector &x,
338  int print_iter = 0, int max_num_iter = 1000,
339  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
340 
341 /// Preconditioned conjugate gradient method. (tolerances are squared)
342 void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x,
343  int print_iter = 0, int max_num_iter = 1000,
344  double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
345 
346 
347 /// GMRES method
349 {
350 protected:
351  int m; // see SetKDim()
352 
353 public:
354  GMRESSolver() { m = 50; }
355 
356 #ifdef MFEM_USE_MPI
357  GMRESSolver(MPI_Comm comm_) : IterativeSolver(comm_) { m = 50; }
358 #endif
359 
360  /// Set the number of iteration to perform between restarts, default is 50.
361  void SetKDim(int dim) { m = dim; }
362 
363  virtual void Mult(const Vector &b, Vector &x) const;
364 };
365 
366 /// FGMRES method
368 {
369 protected:
370  int m;
371 
372 public:
373  FGMRESSolver() { m = 50; }
374 
375 #ifdef MFEM_USE_MPI
376  FGMRESSolver(MPI_Comm comm_) : IterativeSolver(comm_) { m = 50; }
377 #endif
378 
379  void SetKDim(int dim) { m = dim; }
380 
381  virtual void Mult(const Vector &b, Vector &x) const;
382 };
383 
384 /// GMRES method. (tolerances are squared)
385 int GMRES(const Operator &A, Vector &x, const Vector &b, Solver &M,
386  int &max_iter, int m, double &tol, double atol, int printit);
387 
388 /// GMRES method. (tolerances are squared)
389 void GMRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
390  int print_iter = 0, int max_num_iter = 1000, int m = 50,
391  double rtol = 1e-12, double atol = 1e-24);
392 
393 
394 /// BiCGSTAB method
396 {
397 protected:
398  mutable Vector p, phat, s, shat, t, v, r, rtilde;
399 
400  void UpdateVectors();
401 
402 public:
404 
405 #ifdef MFEM_USE_MPI
406  BiCGSTABSolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
407 #endif
408 
409  virtual void SetOperator(const Operator &op)
411 
412  virtual void Mult(const Vector &b, Vector &x) const;
413 };
414 
415 /// BiCGSTAB method. (tolerances are squared)
416 int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M,
417  int &max_iter, double &tol, double atol, int printit);
418 
419 /// BiCGSTAB method. (tolerances are squared)
420 void BiCGSTAB(const Operator &A, Solver &B, const Vector &b, Vector &x,
421  int print_iter = 0, int max_num_iter = 1000,
422  double rtol = 1e-12, double atol = 1e-24);
423 
424 
425 /// MINRES method
427 {
428 protected:
429  mutable Vector v0, v1, w0, w1, q;
430  mutable Vector u1; // used in the preconditioned version
431 
432 public:
434 
435 #ifdef MFEM_USE_MPI
436  MINRESSolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
437 #endif
438 
439  virtual void SetPreconditioner(Solver &pr)
440  {
442  if (oper) { u1.SetSize(width); }
443  }
444 
445  virtual void SetOperator(const Operator &op);
446 
447  virtual void Mult(const Vector &b, Vector &x) const;
448 };
449 
450 /// MINRES method without preconditioner. (tolerances are squared)
451 void MINRES(const Operator &A, const Vector &b, Vector &x, int print_it = 0,
452  int max_it = 1000, double rtol = 1e-12, double atol = 1e-24);
453 
454 /// MINRES method with preconditioner. (tolerances are squared)
455 void MINRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
456  int print_it = 0, int max_it = 1000,
457  double rtol = 1e-12, double atol = 1e-24);
458 
459 
460 /// Newton's method for solving F(x)=b for a given operator F.
461 /** The method GetGradient() must be implemented for the operator F.
462  The preconditioner is used (in non-iterative mode) to evaluate
463  the action of the inverse gradient of the operator. */
465 {
466 protected:
467  mutable Vector r, c;
468  mutable Operator *grad;
469 
470  // Adaptive linear solver rtol variables
471 
472  // Method to determine rtol, 0 means the adaptive algorithm is deactivated.
473  int lin_rtol_type = 0;
474  // rtol to use in first iteration
475  double lin_rtol0;
476  // Maximum rtol
477  double lin_rtol_max;
478  // Function norm ||F(x)|| of the previous iterate
479  mutable double fnorm_last = 0.0;
480  // Linear residual norm of the previous iterate
481  mutable double lnorm_last = 0.0;
482  // Forcing term (linear residual rtol) from the previous iterate
483  mutable double eta_last = 0.0;
484  // Eisenstat-Walker factor gamma
485  double gamma;
486  // Eisenstat-Walker factor alpha
487  double alpha;
488 
489  /** @brief Method for the adaptive linear solver rtol invoked before the
490  linear solve. */
491  void AdaptiveLinRtolPreSolve(const Vector &x,
492  const int it,
493  const double fnorm) const;
494 
495  /** @brief Method for the adaptive linear solver rtol invoked after the
496  linear solve. */
497  void AdaptiveLinRtolPostSolve(const Vector &x,
498  const Vector &b,
499  const int it,
500  const double fnorm) const;
501 
502 public:
504 
505 #ifdef MFEM_USE_MPI
506  NewtonSolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
507 #endif
508  virtual void SetOperator(const Operator &op);
509 
510  /// Set the linear solver for inverting the Jacobian.
511  /** This method is equivalent to calling SetPreconditioner(). */
512  virtual void SetSolver(Solver &solver) { prec = &solver; }
513 
514  /// Solve the nonlinear system with right-hand side @a b.
515  /** If `b.Size() != Height()`, then @a b is assumed to be zero. */
516  virtual void Mult(const Vector &b, Vector &x) const;
517 
518  /** @brief This method can be overloaded in derived classes to implement line
519  search algorithms. */
520  /** The base class implementation (NewtonSolver) simply returns 1. A return
521  value of 0 indicates a failure, interrupting the Newton iteration. */
522  virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const
523  { return 1.0; }
524 
525  /** @brief This method can be overloaded in derived classes to perform
526  computations that need knowledge of the newest Newton state. */
527  virtual void ProcessNewState(const Vector &x) const { }
528 
529  /// Enable adaptive linear solver relative tolerance algorithm.
530  /** Compute a relative tolerance for the Krylov method after each nonlinear
531  iteration, based on the algorithm presented in [1].
532 
533  The maximum linear solver relative tolerance @a rtol_max should be < 1. For
534  @a type 1 the parameters @a alpha and @a gamma are ignored. For @a type 2
535  @a alpha has to be between 0 and 1 and @a gamma between 1 and 2.
536 
537  [1] Eisenstat, Stanley C., and Homer F. Walker. "Choosing the forcing terms
538  in an inexact Newton method."
539  */
540  void SetAdaptiveLinRtol(const int type = 2,
541  const double rtol0 = 0.5,
542  const double rtol_max = 0.9,
543  const double alpha = 0.5 * (1.0 + sqrt(5.0)),
544  const double gamma = 1.0);
545 };
546 
547 /** L-BFGS method for solving F(x)=b for a given operator F, by minimizing
548  the norm of F(x) - b. Requires only the action of the operator F. */
549 class LBFGSSolver : public NewtonSolver
550 {
551 protected:
552  int m = 10;
553 
554 public:
556 
557 #ifdef MFEM_USE_MPI
558  LBFGSSolver(MPI_Comm comm_) : NewtonSolver(comm_) { }
559 #endif
560 
561  void SetHistorySize(int dim) { m = dim; }
562 
563  /// Solve the nonlinear system with right-hand side @a b.
564  /** If `b.Size() != Height()`, then @a b is assumed to be zero. */
565  virtual void Mult(const Vector &b, Vector &x) const;
566 
567  virtual void SetPreconditioner(Solver &pr)
568  { MFEM_WARNING("L-BFGS won't use the given preconditioner."); }
569  virtual void SetSolver(Solver &solver)
570  { MFEM_WARNING("L-BFGS won't use the given solver."); }
571 };
572 
573 /** Adaptive restarted GMRES.
574  m_max and m_min(=1) are the maximal and minimal restart parameters.
575  m_step(=1) is the step to use for going from m_max and m_min.
576  cf(=0.4) is a desired convergence factor. */
577 int aGMRES(const Operator &A, Vector &x, const Vector &b,
578  const Operator &M, int &max_iter,
579  int m_max, int m_min, int m_step, double cf,
580  double &tol, double &atol, int printit);
581 
582 
583 /** Defines operators and constraints for the following optimization problem:
584  *
585  * Find x that minimizes the objective function F(x), subject to
586  * C(x) = c_e,
587  * d_lo <= D(x) <= d_hi,
588  * x_lo <= x <= x_hi.
589  *
590  * The operators F, C, D must take input of the same size (same width).
591  * Gradients of F, C, D might be needed, depending on the OptimizationSolver.
592  * When used with Hiop, gradients of C and D must be DenseMatrices.
593  * F always returns a scalar value, see CalcObjective(), CalcObjectiveGrad().
594  * C and D can have arbitrary heights.
595  * C and D can be NULL, meaning that their constraints are not used.
596  *
597  * When used in parallel, all Vectors are assumed to be true dof vectors, and
598  * the operators are expected to be defined for tdof vectors. */
600 {
601 protected:
602  /// Not owned, some can remain unused (NULL).
603  const Operator *C, *D;
604  const Vector *c_e, *d_lo, *d_hi, *x_lo, *x_hi;
605 
606 public:
607  const int input_size;
608 
609  /// In parallel, insize is the number of the local true dofs.
610  OptimizationProblem(int insize, const Operator *C_, const Operator *D_);
611 
612  /// Objective F(x). In parallel, the result should be reduced over tasks.
613  virtual double CalcObjective(const Vector &x) const = 0;
614  /// The result grad is expected to enter with the correct size.
615  virtual void CalcObjectiveGrad(const Vector &x, Vector &grad) const
616  { MFEM_ABORT("The objective gradient is not implemented."); }
617 
618  void SetEqualityConstraint(const Vector &c);
619  void SetInequalityConstraint(const Vector &dl, const Vector &dh);
620  void SetSolutionBounds(const Vector &xl, const Vector &xh);
621 
622  const Operator *GetC() const { return C; }
623  const Operator *GetD() const { return D; }
624  const Vector *GetEqualityVec() const { return c_e; }
625  const Vector *GetInequalityVec_Lo() const { return d_lo; }
626  const Vector *GetInequalityVec_Hi() const { return d_hi; }
627  const Vector *GetBoundsVec_Lo() const { return x_lo; }
628  const Vector *GetBoundsVec_Hi() const { return x_hi; }
629 
630  int GetNumConstraints() const;
631 };
632 
633 /// Abstract solver for OptimizationProblems.
635 {
636 protected:
638 
639 public:
641 #ifdef MFEM_USE_MPI
642  OptimizationSolver(MPI_Comm comm_): IterativeSolver(comm_), problem(NULL) { }
643 #endif
644  virtual ~OptimizationSolver() { }
645 
646  /** This function is virtual as solvers might need to perform some initial
647  * actions (e.g. validation) with the OptimizationProblem. */
649  { problem = &prob; }
650 
651  virtual void Mult(const Vector &xt, Vector &x) const = 0;
652 
653  virtual void SetPreconditioner(Solver &pr)
654  { MFEM_ABORT("Not meaningful for this solver."); }
655  virtual void SetOperator(const Operator &op)
656  { MFEM_ABORT("Not meaningful for this solver."); }
657 };
658 
659 /** SLBQP optimizer:
660  * (S)ingle (L)inearly Constrained with (B)ounds (Q)uadratic (P)rogram
661  *
662  * Minimize || x-x_t ||, subject to
663  * sum w_i x_i = a,
664  * x_lo <= x <= x_hi.
665  */
667 {
668 protected:
670  double a;
671 
672  /// Solve QP at fixed lambda
673  inline double solve(double l, const Vector &xt, Vector &x, int &nclip) const
674  {
675  add(xt, l, w, x);
676  if (problem == NULL) { x.median(lo,hi); }
677  else
678  {
681  }
682  nclip++;
683  if (problem == NULL) { return Dot(w, x) - a; }
684  else
685  {
686  Vector c(1);
687  // Includes parallel communication.
688  problem->GetC()->Mult(x, c);
689 
690  return c(0) - (*problem->GetEqualityVec())(0);
691  }
692  }
693 
694  inline void print_iteration(int it, double r, double l) const;
695 
696 public:
698 
699 #ifdef MFEM_USE_MPI
700  SLBQPOptimizer(MPI_Comm comm_) : OptimizationSolver(comm_) { }
701 #endif
702 
703  /** Setting an OptimizationProblem will overwrite the Vectors given by
704  * SetBounds and SetLinearConstraint. The objective function remains
705  * unchanged. */
706  virtual void SetOptimizationProblem(const OptimizationProblem &prob);
707 
708  void SetBounds(const Vector &lo_, const Vector &hi_);
709  void SetLinearConstraint(const Vector &w_, double a_);
710 
711  /** We let the target values play the role of the initial vector xt, from
712  * which the operator generates the optimal vector x. */
713  virtual void Mult(const Vector &xt, Vector &x) const;
714 };
715 
716 /** Block ILU solver:
717  * Performs a block ILU(k) approximate factorization with specified block
718  * size. Currently only k=0 is supported. This is useful as a preconditioner
719  * for DG-type discretizations, where the system matrix has a natural
720  * (elemental) block structure.
721  *
722  * In the case of DG discretizations, the block size should usually be set to
723  * either ndofs_per_element or vdim*ndofs_per_element (if the finite element
724  * space has Ordering::byVDIM). The block size must evenly divide the size of
725  * the matrix.
726  *
727  * Renumbering the blocks is also supported by specifying a reordering method.
728  * Currently greedy minimum discarded fill ordering and no reordering are
729  * supported. Renumbering the blocks can lead to a much better approximate
730  * factorization.
731  */
732 class BlockILU : public Solver
733 {
734 public:
735 
736  /// The reordering method used by the BlockILU factorization.
737  enum class Reordering
738  {
740  NONE
741  };
742 
743  /** Create an "empty" BlockILU solver. SetOperator must be called later to
744  * actually form the factorization
745  */
746  BlockILU(int block_size_,
748  int k_fill_ = 0);
749 
750  /** Create a block ILU approximate factorization for the matrix @a op.
751  * @a op should be of type either SparseMatrix or HypreParMatrix. In the
752  * case that @a op is a HypreParMatrix, the ILU factorization is performed
753  * on the diagonal blocks of the parallel decomposition.
754  */
755  BlockILU(const Operator &op, int block_size_ = 1,
757  int k_fill_ = 0);
758 
759  /** Perform the block ILU factorization for the matrix @a op.
760  * As in the constructor, @a op must either be a SparseMatrix or
761  * HypreParMatrix
762  */
763  void SetOperator(const Operator &op);
764 
765  /// Solve the system `LUx = b`, where `L` and `U` are the block ILU factors.
766  void Mult(const Vector &b, Vector &x) const;
767 
768  /** Get the I array for the block CSR representation of the factorization.
769  * Similar to SparseMatrix::GetI(). Mostly used for testing.
770  */
771  int *GetBlockI() { return IB.GetData(); }
772 
773  /** Get the J array for the block CSR representation of the factorization.
774  * Similar to SparseMatrix::GetJ(). Mostly used for testing.
775  */
776  int *GetBlockJ() { return JB.GetData(); }
777 
778  /** Get the data array for the block CSR representation of the factorization.
779  * Similar to SparseMatrix::GetData(). Mostly used for testing.
780  */
781  double *GetBlockData() { return AB.Data(); }
782 
783 private:
784  /// Set up the block CSR structure corresponding to a sparse matrix @a A
785  void CreateBlockPattern(const class SparseMatrix &A);
786 
787  /// Perform the block ILU factorization
788  void Factorize();
789 
790  int block_size;
791 
792  /// Fill level for block ILU(k) factorizations. Only k=0 is supported.
793  int k_fill;
794 
795  Reordering reordering;
796 
797  /// Temporary vector used in the Mult() function.
798  mutable Vector y;
799 
800  /// Permutation and inverse permutation vectors for the block reordering.
801  Array<int> P, Pinv;
802 
803  /** Block CSR storage of the factorization. The block upper triangular part
804  * stores the U factor. The L factor implicitly has identity on the diagonal
805  * blocks, and the rest of L is given by the strictly block lower triangular
806  * part.
807  */
808  Array<int> IB, ID, JB;
809  DenseTensor AB;
810 
811  /// DB(i) stores the LU factorization of the i'th diagonal block
812  mutable DenseTensor DB;
813  /// Pivot arrays for the LU factorizations given by #DB
814  mutable Array<int> ipiv;
815 };
816 
817 
818 /// Monitor that checks whether the residual is zero at a given set of dofs.
819 /** This monitor is useful for checking if the initial guess, rhs, operator, and
820  preconditioner are properly setup for solving in the subspace with imposed
821  essential boundary conditions. */
823 {
824 protected:
825  const Array<int> *ess_dofs_list; ///< Not owned
826 
827 public:
828  ResidualBCMonitor(const Array<int> &ess_dofs_list_)
829  : ess_dofs_list(&ess_dofs_list_) { }
830 
831  void MonitorResidual(int it, double norm, const Vector &r,
832  bool final) override;
833 };
834 
835 
836 #ifdef MFEM_USE_SUITESPARSE
837 
838 /// Direct sparse solver using UMFPACK
839 class UMFPackSolver : public Solver
840 {
841 protected:
844  void *Numeric;
845  SuiteSparse_long *AI, *AJ;
846 
847  void Init();
848 
849 public:
850  double Control[UMFPACK_CONTROL];
851  mutable double Info[UMFPACK_INFO];
852 
853  /** @brief For larger matrices, if the solver fails, set the parameter @a
854  use_long_ints_ = true. */
855  UMFPackSolver(bool use_long_ints_ = false)
856  : use_long_ints(use_long_ints_) { Init(); }
857  /** @brief Factorize the given SparseMatrix using the defaults. For larger
858  matrices, if the solver fails, set the parameter @a use_long_ints_ =
859  true. */
860  UMFPackSolver(SparseMatrix &A, bool use_long_ints_ = false)
861  : use_long_ints(use_long_ints_) { Init(); SetOperator(A); }
862 
863  /** @brief Factorize the given Operator @a op which must be a SparseMatrix.
864 
865  The factorization uses the parameters set in the #Control data member.
866  @note This method calls SparseMatrix::SortColumnIndices() with @a op,
867  modifying the matrix if the column indices are not already sorted. */
868  virtual void SetOperator(const Operator &op);
869 
870  /// Set the print level field in the #Control data member.
871  void SetPrintLevel(int print_lvl) { Control[UMFPACK_PRL] = print_lvl; }
872 
873  virtual void Mult(const Vector &b, Vector &x) const;
874  virtual void MultTranspose(const Vector &b, Vector &x) const;
875 
876  virtual ~UMFPackSolver();
877 };
878 
879 /// Direct sparse solver using KLU
880 class KLUSolver : public Solver
881 {
882 protected:
884  klu_symbolic *Symbolic;
885  klu_numeric *Numeric;
886 
887  void Init();
888 
889 public:
891  : mat(0),Symbolic(0),Numeric(0)
892  { Init(); }
894  : mat(0),Symbolic(0),Numeric(0)
895  { Init(); SetOperator(A); }
896 
897  // Works on sparse matrices only; calls SparseMatrix::SortColumnIndices().
898  virtual void SetOperator(const Operator &op);
899 
900  virtual void Mult(const Vector &b, Vector &x) const;
901  virtual void MultTranspose(const Vector &b, Vector &x) const;
902 
903  virtual ~KLUSolver();
904 
905  mutable klu_common Common;
906 };
907 
908 #endif // MFEM_USE_SUITESPARSE
909 
910 /// Block diagonal solver for A, each block is inverted by direct solver
912 {
913  SparseMatrix& block_dof;
914  mutable Array<int> local_dofs;
915  mutable Vector sub_rhs;
916  mutable Vector sub_sol;
917  std::unique_ptr<DenseMatrixInverse[]> block_solvers;
918 public:
919  /// block_dof is a boolean matrix, block_dof(i, j) = 1 if j-th dof belongs to
920  /// i-th block, block_dof(i, j) = 0 otherwise.
921  DirectSubBlockSolver(const SparseMatrix& A, const SparseMatrix& block_dof);
922  virtual void Mult(const Vector &x, Vector &y) const;
923  virtual void SetOperator(const Operator &op) { }
924 };
925 
926 /// Solver S such that I - A * S = (I - A * S1) * (I - A * S0).
927 /// That is, S = S0 + S1 - S1 * A * S0.
928 class ProductSolver : public Solver
929 {
930  OperatorPtr A;
931  OperatorPtr S0;
932  OperatorPtr S1;
933 public:
935  bool ownA, bool ownS0, bool ownS1)
936  : Solver(A_->NumRows()), A(A_, ownA), S0(S0_, ownS0), S1(S1_, ownS1) { }
937  virtual void Mult(const Vector &x, Vector &y) const;
938  virtual void MultTranspose(const Vector &x, Vector &y) const;
939  virtual void SetOperator(const Operator &op) { }
940 };
941 
942 #ifdef MFEM_USE_MPI
943 /** This smoother does relaxations on an auxiliary space (determined by a map
944  from the original space to the auxiliary space provided by the user).
945  The smoother on the auxiliary space is a HypreSmoother. Its options can be
946  modified through GetSmoother.
947  For example, the space can be the nullspace of div/curl, in which case the
948  smoother can be used to construct a Hiptmair smoother. */
949 class AuxSpaceSmoother : public Solver
950 {
951  OperatorPtr aux_map_;
952  OperatorPtr aux_system_;
953  OperatorPtr aux_smoother_;
954  void Mult(const Vector &x, Vector &y, bool transpose) const;
955 public:
956  AuxSpaceSmoother(const HypreParMatrix &op, HypreParMatrix *aux_map,
957  bool op_is_symmetric = true, bool own_aux_map = false);
958  virtual void Mult(const Vector &x, Vector &y) const { Mult(x, y, false); }
959  virtual void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y, true); }
960  virtual void SetOperator(const Operator &op) { }
961  HypreSmoother& GetSmoother() { return *aux_smoother_.As<HypreSmoother>(); }
962 };
963 #endif // MFEM_USE_MPI
964 
965 }
966 
967 #endif // MFEM_SOLVERS
const Array< int > * ess_dofs_list
Not owned.
Definition: solvers.hpp:825
MINRESSolver(MPI_Comm comm_)
Definition: solvers.hpp:436
const Vector * GetInequalityVec_Hi() const
Definition: solvers.hpp:626
GMRESSolver(MPI_Comm comm_)
Definition: solvers.hpp:357
Conjugate gradient method.
Definition: solvers.hpp:316
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:99
CGSolver(MPI_Comm comm_)
Definition: solvers.hpp:327
HypreSmoother & GetSmoother()
Definition: solvers.hpp:961
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:3032
double Info[UMFPACK_INFO]
Definition: solvers.hpp:851
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:2801
Chebyshev accelerated smoothing with given vector, no matrix necessary.
Definition: solvers.hpp:202
SparseMatrix * mat
Definition: solvers.hpp:843
int GetNumIterations() const
Definition: solvers.hpp:103
const Operator * oper
Definition: solvers.hpp:75
MPI_Comm GetComm() const
Return the associated MPI communicator, or MPI_COMM_NULL if no communicator is set.
Definition: solvers.hpp:120
SuiteSparse_long * AI
Definition: solvers.hpp:845
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:616
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:3173
void SetHistorySize(int dim)
Definition: solvers.hpp:561
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:655
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:109
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:527
OperatorJacobiSmoother(const double damping=1.0)
Default constructor: the diagonal will be computed by subsequent calls to SetOperator() using the Ope...
Definition: solvers.cpp:119
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:108
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:210
FGMRES method.
Definition: solvers.hpp:367
int * GetBlockI()
Definition: solvers.hpp:771
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1071
Direct sparse solver using KLU.
Definition: solvers.hpp:880
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:512
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:3134
void SetOperator(const Operator &op)
Definition: solvers.cpp:2570
MINRES method.
Definition: solvers.hpp:426
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:860
const Operator * GetC() const
Definition: solvers.hpp:622
Reordering
The reordering method used by the BlockILU factorization.
Definition: solvers.hpp:737
virtual ~IterativeSolverMonitor()
Definition: solvers.hpp:45
ResidualBCMonitor(const Array< int > &ess_dofs_list_)
Definition: solvers.hpp:828
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.hpp:939
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:295
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:291
virtual ~OptimizationSolver()
Definition: solvers.hpp:644
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:260
double Norm(const Vector &x) const
Definition: solvers.hpp:87
const Vector * d_hi
Definition: solvers.hpp:604
double GetFinalNorm() const
Definition: solvers.hpp:105
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:1416
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:432
void SetKDim(int dim)
Definition: solvers.hpp:379
const Vector * x_lo
Definition: solvers.hpp:604
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:839
Operator * grad
Definition: solvers.hpp:468
const Vector * x_hi
Definition: solvers.hpp:604
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:884
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1453
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:439
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:71
const Vector * GetBoundsVec_Lo() const
Definition: solvers.hpp:627
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:98
double solve(double l, const Vector &xt, Vector &x, int &nclip) const
Solve QP at fixed lambda.
Definition: solvers.hpp:673
int GetNumConstraints() const
Definition: solvers.cpp:2147
Data type sparse matrix.
Definition: sparsemat.hpp:41
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:1614
virtual ~UMFPackSolver()
Definition: solvers.cpp:3067
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:130
double b
Definition: lissajous.cpp:42
void SetMaxIter(int max_it)
Definition: solvers.hpp:100
OperatorChebyshevSmoother(const Operator &oper_, const Vector &d, const Array< int > &ess_tdof_list, int order, double max_eig_estimate)
Definition: solvers.cpp:235
void SetLinearConstraint(const Vector &w_, double a_)
Definition: solvers.cpp:2174
void median(const Vector &lo, const Vector &hi)
v = median(v,lo,hi) entrywise. Implementation assumes lo &lt;= hi.
Definition: vector.cpp:501
DirectSubBlockSolver(const SparseMatrix &A, const SparseMatrix &block_dof)
Definition: solvers.cpp:3158
SuiteSparse_long * AJ
Definition: solvers.hpp:845
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition: solvers.hpp:361
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:464
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:802
void print_iteration(int it, double r, double l) const
Definition: solvers.cpp:2180
void Setup(const Vector &diag)
Definition: solvers.cpp:196
Parallel smoothers in hypre.
Definition: hypre.hpp:840
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:817
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.cpp:3089
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1439
double Dot(const Vector &x, const Vector &y) const
Definition: solvers.cpp:55
int GetConverged() const
Definition: solvers.hpp:104
prob_type prob
Definition: ex25.cpp:153
klu_symbolic * Symbolic
Definition: solvers.hpp:884
void SetIterativeSolver(const IterativeSolver &solver)
This method is invoked by ItertiveSolver::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:960
const Vector * GetBoundsVec_Hi() const
Definition: solvers.hpp:628
BlockILU(int block_size_, Reordering reordering_=Reordering::MINIMUM_DISCARDED_FILL, int k_fill_=0)
Definition: solvers.cpp:2552
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1653
void SetPrintLevel(int print_lvl)
Set the print level field in the Control data member.
Definition: solvers.hpp:871
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:1805
Stationary linear iteration: x &lt;- x + B (b - A x)
Definition: solvers.hpp:284
void SetInequalityConstraint(const Vector &dl, const Vector &dh)
Definition: solvers.cpp:2129
void SetEqualityConstraint(const Vector &c)
Definition: solvers.cpp:2121
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:653
void SetOperator(const Operator &op_)
Set/update the solver for the given operator.
Definition: solvers.hpp:262
ProductSolver(Operator *A_, Solver *S0_, Solver *S1_, bool ownA, bool ownS0, bool ownS1)
Definition: solvers.hpp:934
void SetAbsTol(double atol)
Definition: solvers.hpp:99
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:3205
void SetRelTol(double rtol)
Definition: solvers.hpp:98
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:850
Monitor that checks whether the residual is zero at a given set of dofs.
Definition: solvers.hpp:822
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:855
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1825
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1280
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: operator.hpp:69
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3120
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:1745
LBFGSSolver(MPI_Comm comm_)
Definition: solvers.hpp:558
double * GetBlockData()
Definition: solvers.hpp:781
GMRES method.
Definition: solvers.hpp:348
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: solvers.hpp:569
const Vector * GetEqualityVec() const
Definition: solvers.hpp:624
FGMRESSolver(MPI_Comm comm_)
Definition: solvers.hpp:376
virtual void SetOptimizationProblem(const OptimizationProblem &prob)
Definition: solvers.cpp:2155
void SetMonitor(IterativeSolverMonitor &m)
Set the iterative solver monitor.
Definition: solvers.hpp:114
const Operator * GetD() const
Definition: solvers.hpp:623
double a
Definition: lissajous.cpp:41
void UpdateVectors()
Definition: solvers.cpp:426
Abstract solver for OptimizationProblems.
Definition: solvers.hpp:634
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
klu_common Common
Definition: solvers.hpp:905
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:522
OptimizationProblem(int insize, const Operator *C_, const Operator *D_)
In parallel, insize is the number of the local true dofs.
Definition: solvers.cpp:2111
void SetSolutionBounds(const Vector &xl, const Vector &xh)
Definition: solvers.cpp:2139
const Vector * GetInequalityVec_Lo() const
Definition: solvers.hpp:625
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.hpp:958
AuxSpaceSmoother(const HypreParMatrix &op, HypreParMatrix *aux_map, bool op_is_symmetric=true, bool own_aux_map=false)
Definition: solvers.cpp:3223
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:386
virtual void Mult(const Vector &xt, Vector &x) const =0
Operator application: y=A(x).
BiCGSTABSolver(MPI_Comm comm_)
Definition: solvers.hpp:406
void SetOperator(const Operator &op)
Recompute the diagonal using the method AssembleDiagonal of the given new Operator, op.
Definition: solvers.cpp:160
int dim
Definition: ex24.cpp:53
SLBQPOptimizer(MPI_Comm comm_)
Definition: solvers.hpp:700
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:1758
int * GetBlockJ()
Definition: solvers.hpp:776
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:409
virtual void Mult(const Vector &xt, Vector &x) const
Definition: solvers.cpp:2187
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:167
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:330
virtual void CalcObjectiveGrad(const Vector &x, Vector &grad) const
The result grad is expected to enter with the correct size.
Definition: solvers.hpp:615
Vector data type.
Definition: vector.hpp:60
void MonitorResidual(int it, double norm, const Vector &r, bool final) override
Monitor the residual vector r.
Definition: solvers.cpp:2852
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:959
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:567
const Vector * d_lo
Definition: solvers.hpp:604
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:92
KLUSolver(SparseMatrix &A)
Definition: solvers.hpp:893
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:1960
virtual ~KLUSolver()
Definition: solvers.cpp:3148
double * Data()
Definition: densemat.hpp:849
void UpdateVectors()
Definition: solvers.cpp:607
BiCGSTAB method.
Definition: solvers.hpp:395
OptimizationSolver(MPI_Comm comm_)
Definition: solvers.hpp:642
SparseMatrix * mat
Definition: solvers.hpp:883
Base class for solvers.
Definition: operator.hpp:648
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: solvers.hpp:923
const OptimizationProblem * problem
Definition: solvers.hpp:637
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:277
Block diagonal solver for A, each block is inverted by direct solver.
Definition: solvers.hpp:911
void SetBounds(const Vector &lo_, const Vector &hi_)
Definition: solvers.cpp:2168
const Operator * D
Definition: solvers.hpp:603
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:745
virtual void SetOptimizationProblem(const OptimizationProblem &prob)
Definition: solvers.hpp:648
NewtonSolver(MPI_Comm comm_)
Definition: solvers.hpp:506
const Vector * c_e
Definition: solvers.hpp:604
IterativeSolverMonitor * monitor
Definition: solvers.hpp:77
klu_numeric * Numeric
Definition: solvers.hpp:885
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:3188
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:2999
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:298
const Operator * C
Not owned, some can remain unused (NULL).
Definition: solvers.hpp:603
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:575
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:1242
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1642
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:2904