15#ifndef dealii_solver_minres_h
16#define dealii_solver_minres_h
69template <
typename VectorType = Vector<
double>>
103 template <
typename MatrixType,
typename PreconditionerType>
107 void solve(const MatrixType &A,
110 const PreconditionerType &preconditioner);
121 "The preconditioner for MinRes must be a symmetric and "
122 "definite operator, even though MinRes can solve linear "
123 "systems with symmetric and *indefinite* operators. "
124 "During iterations, MinRes has detected that the "
125 "preconditioner is apparently not definite.");
141 print_vectors(const
unsigned int step,
144 const VectorType &d) const;
160template <
typename VectorType>
166 , res2(
numbers::signaling_nan<double>())
171template <
typename VectorType>
174 const AdditionalData &)
176 , res2(
numbers::signaling_nan<double>())
181template <
typename VectorType>
189template <
typename VectorType>
194 const VectorType &)
const
199template <
typename VectorType>
201template <
typename MatrixType,
typename PreconditionerType>
208 const PreconditionerType &preconditioner)
225 vecptr u[3] = {Vu0.get(), Vu1.get(), Vu2.get()};
226 vecptr m[3] = {Vm0.get(), Vm1.get(), Vm2.get()};
231 u[0]->reinit(b,
true);
232 u[1]->reinit(b,
true);
233 u[2]->reinit(b,
true);
234 m[0]->reinit(b,
true);
235 m[1]->reinit(b,
true);
236 m[2]->reinit(b,
true);
240 double delta[3] = {0, 0, 0};
241 double f[2] = {0, 0};
242 double e[2] = {0, 0};
264 preconditioner.vmult(v, *u[1]);
266 delta[1] = v * (*u[1]);
268 Assert(delta[1] >= 0, ExcPreconditionerNotDefinite());
289 u[2]->add(-
std::sqrt(delta[1] / delta[0]), *u[0]);
291 const double gamma = *u[2] * v;
292 u[2]->add(-gamma /
std::sqrt(delta[1]), *u[1]);
298 preconditioner.vmult(v, *u[2]);
300 delta[2] = v * (*u[2]);
302 Assert(delta[2] >= 0, ExcPreconditionerNotDefinite());
311 d_ = s *
e[0] - c *
gamma;
317 const double d =
std::sqrt(d_ * d_ + delta[2]);
329 m[0]->add(-e[0], *m[1]);
331 m[0]->add(-f[0], *m[2]);
334 r_l2 *= std::fabs(s);
336 conv = this->iteration_status(j, r_l2, x);
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
virtual ~SolverMinRes() override=default
SolverMinRes(SolverControl &cn, const AdditionalData &data=AdditionalData())
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
virtual double criterion()
SolverMinRes(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define DeclExceptionMsg(Exception, defaulttext)
#define AssertThrow(cond, exc)
std::vector< index_type > data
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
long double gamma(const unsigned int n)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)