Reference documentation for deal.II version GIT relicensing-1291-gf3a46e514e 2024-07-24 11:20:02+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\}}\)
Loading...
Searching...
No Matches
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
86 const AdditionalData &data = AdditionalData());
87
93 const AdditionalData &data = AdditionalData());
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 DeclException0(ExcPreconditionerNotDefinite);
123protected:
127 virtual double
128 criterion();
129
135 virtual void
136 print_vectors(const unsigned int step,
137 const VectorType &x,
138 const VectorType &r,
139 const VectorType &d) const;
140
147 double res2;
148};
149
151/*------------------------- Implementation ----------------------------*/
152
153#ifndef DOXYGEN
154
155template <typename VectorType>
159 const AdditionalData &)
160 : SolverBase<VectorType>(cn, mem)
161 , res2(numbers::signaling_nan<double>())
162{}
163
164
165
166template <typename VectorType>
169 const AdditionalData &)
170 : SolverBase<VectorType>(cn)
171 , res2(numbers::signaling_nan<double>())
172{}
173
174
175
176template <typename VectorType>
179{
180 return res2;
181}
182
183
184template <typename VectorType>
186void SolverMinRes<VectorType>::print_vectors(const unsigned int,
187 const VectorType &,
188 const VectorType &,
189 const VectorType &) const
190{}
191
192
193
194template <typename VectorType>
196template <typename MatrixType, typename PreconditionerType>
200void SolverMinRes<VectorType>::solve(const MatrixType &A,
201 VectorType &x,
202 const VectorType &b,
203 const PreconditionerType &preconditioner)
204{
205 LogStream::Prefix prefix("minres");
206
207 // Memory allocation
208 typename VectorMemory<VectorType>::Pointer Vu0(this->memory);
209 typename VectorMemory<VectorType>::Pointer Vu1(this->memory);
210 typename VectorMemory<VectorType>::Pointer Vu2(this->memory);
211
212 typename VectorMemory<VectorType>::Pointer Vm0(this->memory);
213 typename VectorMemory<VectorType>::Pointer Vm1(this->memory);
214 typename VectorMemory<VectorType>::Pointer Vm2(this->memory);
215
216 typename VectorMemory<VectorType>::Pointer Vv(this->memory);
217
218 // define some aliases for simpler access
219 using vecptr = VectorType *;
220 vecptr u[3] = {Vu0.get(), Vu1.get(), Vu2.get()};
221 vecptr m[3] = {Vm0.get(), Vm1.get(), Vm2.get()};
222 VectorType &v = *Vv;
223
224 // resize the vectors, but do not set the values since they'd be overwritten
225 // soon anyway.
226 u[0]->reinit(b, true);
227 u[1]->reinit(b, true);
228 u[2]->reinit(b, true);
229 m[0]->reinit(b, true);
230 m[1]->reinit(b, true);
231 m[2]->reinit(b, true);
232 v.reinit(b, true);
233
234 // some values needed
235 double delta[3] = {0, 0, 0};
236 double f[2] = {0, 0};
237 double e[2] = {0, 0};
238
239 double r_l2 = 0;
240 double r0 = 0;
241 double tau = 0;
242 double c = 0;
243 double s = 0;
244 double d_ = 0;
245
246 // The iteration step.
247 unsigned int j = 1;
248
249
250 // Start of the solution process
251 A.vmult(*m[0], x);
252 *u[1] = b;
253 *u[1] -= *m[0];
254 // Precondition is applied.
255 // The preconditioner has to be
256 // positive definite and symmetric
257
258 // M v = u[1]
259 preconditioner.vmult(v, *u[1]);
260
261 delta[1] = v * (*u[1]);
262 // Preconditioner positive
263 Assert(delta[1] >= 0, ExcPreconditionerNotDefinite());
264
265 r0 = std::sqrt(delta[1]);
266 r_l2 = r0;
267
268
269 u[0]->reinit(b);
270 delta[0] = 1.;
271 m[0]->reinit(b);
272 m[1]->reinit(b);
273 m[2]->reinit(b);
274
275 SolverControl::State conv = this->iteration_status(0, r_l2, x);
276 while (conv == SolverControl::iterate)
277 {
278 if (delta[1] != 0)
279 v *= 1. / std::sqrt(delta[1]);
280 else
281 v.reinit(b);
282
283 A.vmult(*u[2], v);
284 u[2]->add(-std::sqrt(delta[1] / delta[0]), *u[0]);
285
286 const double gamma = *u[2] * v;
287 u[2]->add(-gamma / std::sqrt(delta[1]), *u[1]);
288 *m[0] = v;
289
290 // precondition: solve M v = u[2]
291 // Preconditioner has to be positive
292 // definite and symmetric.
293 preconditioner.vmult(v, *u[2]);
294
295 delta[2] = v * (*u[2]);
296
297 Assert(delta[2] >= 0, ExcPreconditionerNotDefinite());
298
299 if (j == 1)
300 {
301 d_ = gamma;
302 e[1] = std::sqrt(delta[2]);
303 }
304 if (j > 1)
305 {
306 d_ = s * e[0] - c * gamma;
307 e[0] = c * e[0] + s * gamma;
308 f[1] = s * std::sqrt(delta[2]);
309 e[1] = -c * std::sqrt(delta[2]);
310 }
311
312 const double d = std::sqrt(d_ * d_ + delta[2]);
313
314 if (j > 1)
315 tau *= s / c;
316 c = d_ / d;
317 tau *= c;
318
319 s = std::sqrt(delta[2]) / d;
320
321 if (j == 1)
322 tau = r0 * c;
323
324 m[0]->add(-e[0], *m[1]);
325 if (j > 1)
326 m[0]->add(-f[0], *m[2]);
327 *m[0] *= 1. / d;
328 x.add(tau, *m[0]);
329 r_l2 *= std::fabs(s);
330
331 conv = this->iteration_status(j, r_l2, x);
332
333 // next iteration step
334 ++j;
335 // All vectors have to be shifted
336 // one iteration step.
337 // This should be changed one time.
338 swap(*m[2], *m[1]);
339 swap(*m[1], *m[0]);
340
341 // likewise, but reverse direction:
342 // u[0] = u[1];
343 // u[1] = u[2];
344 swap(*u[0], *u[1]);
345 swap(*u[1], *u[2]);
346
347 // these are scalars, so need
348 // to bother
349 f[0] = f[1];
350 e[0] = e[1];
351 delta[0] = delta[1];
352 delta[1] = delta[2];
353 }
354
355 // in case of failure: throw exception
358
359 // otherwise exit as normal
360}
361
362#endif // DOXYGEN
363
365
366#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:503
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:177
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#define DeclException0(Exception0)
Definition exceptions.h:471
#define Assert(cond, exc)
#define AssertThrow(cond, exc)
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(SmartPointer< T, P > &t1, SmartPointer< T, Q > &t2)