Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3067-g24f3489fdc 2025-04-14 21:10:00+00:00
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
solver_minres.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2000 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_solver_minres_h
16#define dealii_solver_minres_h
17
18
19#include <deal.II/base/config.h>
20
25
26#include <deal.II/lac/solver.h>
28
29#include <cmath>
30
32
69template <typename VectorType = Vector<double>>
71class SolverMinRes : public SolverBase<VectorType>
72{
73public:
79 {};
80
87
94
98 virtual ~SolverMinRes() override = default;
99
103 template <typename MatrixType, typename PreconditionerType>
107 void solve(const MatrixType &A,
108 VectorType &x,
109 const VectorType &b,
110 const PreconditionerType &preconditioner);
111
120 DeclExceptionMsg(ExcPreconditionerNotDefinite,
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.");
128protected:
132 virtual double
133 criterion();
134
140 virtual void
141 print_vectors(const unsigned int step,
142 const VectorType &x,
143 const VectorType &r,
144 const VectorType &d) const;
145
152 double res2;
153};
154
156/*------------------------- Implementation ----------------------------*/
157
158#ifndef DOXYGEN
159
160template <typename VectorType>
164 const AdditionalData &)
165 : SolverBase<VectorType>(cn, mem)
166 , res2(numbers::signaling_nan<double>())
167{}
168
169
170
171template <typename VectorType>
174 const AdditionalData &)
175 : SolverBase<VectorType>(cn)
176 , res2(numbers::signaling_nan<double>())
177{}
178
179
180
181template <typename VectorType>
184{
185 return res2;
186}
187
188
189template <typename VectorType>
191void SolverMinRes<VectorType>::print_vectors(const unsigned int,
192 const VectorType &,
193 const VectorType &,
194 const VectorType &) const
195{}
196
197
198
199template <typename VectorType>
201template <typename MatrixType, typename PreconditionerType>
205void SolverMinRes<VectorType>::solve(const MatrixType &A,
206 VectorType &x,
207 const VectorType &b,
208 const PreconditionerType &preconditioner)
209{
210 LogStream::Prefix prefix("minres");
211
212 // Memory allocation
213 typename VectorMemory<VectorType>::Pointer Vu0(this->memory);
214 typename VectorMemory<VectorType>::Pointer Vu1(this->memory);
215 typename VectorMemory<VectorType>::Pointer Vu2(this->memory);
216
217 typename VectorMemory<VectorType>::Pointer Vm0(this->memory);
218 typename VectorMemory<VectorType>::Pointer Vm1(this->memory);
219 typename VectorMemory<VectorType>::Pointer Vm2(this->memory);
220
221 typename VectorMemory<VectorType>::Pointer Vv(this->memory);
222
223 // define some aliases for simpler access
224 using vecptr = VectorType *;
225 vecptr u[3] = {Vu0.get(), Vu1.get(), Vu2.get()};
226 vecptr m[3] = {Vm0.get(), Vm1.get(), Vm2.get()};
227 VectorType &v = *Vv;
228
229 // resize the vectors, but do not set the values since they'd be overwritten
230 // soon anyway.
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);
237 v.reinit(b, true);
238
239 // some values needed
240 double delta[3] = {0, 0, 0};
241 double f[2] = {0, 0};
242 double e[2] = {0, 0};
243
244 double r_l2 = 0;
245 double r0 = 0;
246 double tau = 0;
247 double c = 0;
248 double s = 0;
249 double d_ = 0;
250
251 // The iteration step.
252 unsigned int j = 1;
253
254
255 // Start of the solution process
256 A.vmult(*m[0], x);
257 *u[1] = b;
258 *u[1] -= *m[0];
259 // Precondition is applied.
260 // The preconditioner has to be
261 // positive definite and symmetric
262
263 // M v = u[1]
264 preconditioner.vmult(v, *u[1]);
265
266 delta[1] = v * (*u[1]);
267 // Preconditioner positive
268 Assert(delta[1] >= 0, ExcPreconditionerNotDefinite());
269
270 r0 = std::sqrt(delta[1]);
271 r_l2 = r0;
272
273
274 u[0]->reinit(b);
275 delta[0] = 1.;
276 m[0]->reinit(b);
277 m[1]->reinit(b);
278 m[2]->reinit(b);
279
280 SolverControl::State conv = this->iteration_status(0, r_l2, x);
281 while (conv == SolverControl::iterate)
282 {
283 if (delta[1] != 0)
284 v *= 1. / std::sqrt(delta[1]);
285 else
286 v.reinit(b);
287
288 A.vmult(*u[2], v);
289 u[2]->add(-std::sqrt(delta[1] / delta[0]), *u[0]);
290
291 const double gamma = *u[2] * v;
292 u[2]->add(-gamma / std::sqrt(delta[1]), *u[1]);
293 *m[0] = v;
294
295 // precondition: solve M v = u[2]
296 // Preconditioner has to be positive
297 // definite and symmetric.
298 preconditioner.vmult(v, *u[2]);
299
300 delta[2] = v * (*u[2]);
301
302 Assert(delta[2] >= 0, ExcPreconditionerNotDefinite());
303
304 if (j == 1)
305 {
306 d_ = gamma;
307 e[1] = std::sqrt(delta[2]);
308 }
309 if (j > 1)
310 {
311 d_ = s * e[0] - c * gamma;
312 e[0] = c * e[0] + s * gamma;
313 f[1] = s * std::sqrt(delta[2]);
314 e[1] = -c * std::sqrt(delta[2]);
315 }
316
317 const double d = std::sqrt(d_ * d_ + delta[2]);
318
319 if (j > 1)
320 tau *= s / c;
321 c = d_ / d;
322 tau *= c;
323
324 s = std::sqrt(delta[2]) / d;
325
326 if (j == 1)
327 tau = r0 * c;
328
329 m[0]->add(-e[0], *m[1]);
330 if (j > 1)
331 m[0]->add(-f[0], *m[2]);
332 *m[0] *= 1. / d;
333 x.add(tau, *m[0]);
334 r_l2 *= std::fabs(s);
335
336 conv = this->iteration_status(j, r_l2, x);
337
338 // next iteration step
339 ++j;
340 // All vectors have to be shifted
341 // one iteration step.
342 // This should be changed one time.
343 swap(*m[2], *m[1]);
344 swap(*m[1], *m[0]);
345
346 // likewise, but reverse direction:
347 // u[0] = u[1];
348 // u[1] = u[2];
349 swap(*u[0], *u[1]);
350 swap(*u[1], *u[2]);
351
352 // these are scalars, so need
353 // to bother
354 f[0] = f[1];
355 e[0] = e[1];
356 delta[0] = delta[1];
357 delta[1] = delta[2];
358 }
359
360 // in case of failure: throw exception
363
364 // otherwise exit as normal
365}
366
367#endif // DOXYGEN
368
370
371#endif
@ 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
Definition config.h:35
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:242
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define Assert(cond, exc)
#define DeclExceptionMsg(Exception, defaulttext)
#define AssertThrow(cond, exc)
std::vector< index_type > data
Definition mpi.cc:746
constexpr char A
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)