Reference documentation for deal.II version GIT relicensing-16-gf4f69b3b42 2024-02-27 08: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 - 2023 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
24
25#include <deal.II/lac/solver.h>
27
28#include <cmath>
29
31
68template <typename VectorType = Vector<double>>
69class SolverMinRes : public SolverBase<VectorType>
70{
71public:
77 {};
78
84 const AdditionalData &data = AdditionalData());
85
91 const AdditionalData &data = AdditionalData());
92
96 virtual ~SolverMinRes() override = default;
97
101 template <typename MatrixType, typename PreconditionerType>
102 void
103 solve(const MatrixType &A,
104 VectorType &x,
105 const VectorType &b,
106 const PreconditionerType &preconditioner);
107
119protected:
123 virtual double
125
131 virtual void
132 print_vectors(const unsigned int step,
133 const VectorType &x,
134 const VectorType &r,
135 const VectorType &d) const;
136
143 double res2;
144};
145
147/*------------------------- Implementation ----------------------------*/
148
149#ifndef DOXYGEN
150
151template <typename VectorType>
154 const AdditionalData &)
155 : SolverBase<VectorType>(cn, mem)
156 , res2(numbers::signaling_nan<double>())
157{}
158
159
160
161template <typename VectorType>
163 const AdditionalData &)
164 : SolverBase<VectorType>(cn)
165 , res2(numbers::signaling_nan<double>())
166{}
167
168
169
170template <typename VectorType>
171double
173{
174 return res2;
175}
176
177
178template <typename VectorType>
179void
181 const VectorType &,
182 const VectorType &,
183 const VectorType &) const
184{}
185
186
187
188template <typename VectorType>
189template <typename MatrixType, typename PreconditionerType>
190void
191SolverMinRes<VectorType>::solve(const MatrixType &A,
192 VectorType &x,
193 const VectorType &b,
194 const PreconditionerType &preconditioner)
195{
196 LogStream::Prefix prefix("minres");
197
198 // Memory allocation
199 typename VectorMemory<VectorType>::Pointer Vu0(this->memory);
200 typename VectorMemory<VectorType>::Pointer Vu1(this->memory);
201 typename VectorMemory<VectorType>::Pointer Vu2(this->memory);
202
203 typename VectorMemory<VectorType>::Pointer Vm0(this->memory);
204 typename VectorMemory<VectorType>::Pointer Vm1(this->memory);
205 typename VectorMemory<VectorType>::Pointer Vm2(this->memory);
206
207 typename VectorMemory<VectorType>::Pointer Vv(this->memory);
208
209 // define some aliases for simpler access
210 using vecptr = VectorType *;
211 vecptr u[3] = {Vu0.get(), Vu1.get(), Vu2.get()};
212 vecptr m[3] = {Vm0.get(), Vm1.get(), Vm2.get()};
213 VectorType &v = *Vv;
214
215 // resize the vectors, but do not set the values since they'd be overwritten
216 // soon anyway.
217 u[0]->reinit(b, true);
218 u[1]->reinit(b, true);
219 u[2]->reinit(b, true);
220 m[0]->reinit(b, true);
221 m[1]->reinit(b, true);
222 m[2]->reinit(b, true);
223 v.reinit(b, true);
224
225 // some values needed
226 double delta[3] = {0, 0, 0};
227 double f[2] = {0, 0};
228 double e[2] = {0, 0};
229
230 double r_l2 = 0;
231 double r0 = 0;
232 double tau = 0;
233 double c = 0;
234 double s = 0;
235 double d_ = 0;
236
237 // The iteration step.
238 unsigned int j = 1;
239
240
241 // Start of the solution process
242 A.vmult(*m[0], x);
243 *u[1] = b;
244 *u[1] -= *m[0];
245 // Precondition is applied.
246 // The preconditioner has to be
247 // positive definite and symmetric
248
249 // M v = u[1]
250 preconditioner.vmult(v, *u[1]);
251
252 delta[1] = v * (*u[1]);
253 // Preconditioner positive
254 Assert(delta[1] >= 0, ExcPreconditionerNotDefinite());
255
256 r0 = std::sqrt(delta[1]);
257 r_l2 = r0;
258
259
260 u[0]->reinit(b);
261 delta[0] = 1.;
262 m[0]->reinit(b);
263 m[1]->reinit(b);
264 m[2]->reinit(b);
265
266 SolverControl::State conv = this->iteration_status(0, r_l2, x);
267 while (conv == SolverControl::iterate)
268 {
269 if (delta[1] != 0)
270 v *= 1. / std::sqrt(delta[1]);
271 else
272 v.reinit(b);
273
274 A.vmult(*u[2], v);
275 u[2]->add(-std::sqrt(delta[1] / delta[0]), *u[0]);
276
277 const double gamma = *u[2] * v;
278 u[2]->add(-gamma / std::sqrt(delta[1]), *u[1]);
279 *m[0] = v;
280
281 // precondition: solve M v = u[2]
282 // Preconditioner has to be positive
283 // definite and symmetric.
284 preconditioner.vmult(v, *u[2]);
285
286 delta[2] = v * (*u[2]);
287
288 Assert(delta[2] >= 0, ExcPreconditionerNotDefinite());
289
290 if (j == 1)
291 {
292 d_ = gamma;
293 e[1] = std::sqrt(delta[2]);
294 }
295 if (j > 1)
296 {
297 d_ = s * e[0] - c * gamma;
298 e[0] = c * e[0] + s * gamma;
299 f[1] = s * std::sqrt(delta[2]);
300 e[1] = -c * std::sqrt(delta[2]);
301 }
302
303 const double d = std::sqrt(d_ * d_ + delta[2]);
304
305 if (j > 1)
306 tau *= s / c;
307 c = d_ / d;
308 tau *= c;
309
310 s = std::sqrt(delta[2]) / d;
311
312 if (j == 1)
313 tau = r0 * c;
314
315 m[0]->add(-e[0], *m[1]);
316 if (j > 1)
317 m[0]->add(-f[0], *m[2]);
318 *m[0] *= 1. / d;
319 x.add(tau, *m[0]);
320 r_l2 *= std::fabs(s);
321
322 conv = this->iteration_status(j, r_l2, x);
323
324 // next iteration step
325 ++j;
326 // All vectors have to be shifted
327 // one iteration step.
328 // This should be changed one time.
329 swap(*m[2], *m[1]);
330 swap(*m[1], *m[0]);
331
332 // likewise, but reverse direction:
333 // u[0] = u[1];
334 // u[1] = u[2];
335 swap(*u[0], *u[1]);
336 swap(*u[1], *u[2]);
337
338 // these are scalars, so need
339 // to bother
340 f[0] = f[1];
341 e[0] = e[1];
342 delta[0] = delta[1];
343 delta[1] = delta[2];
344 }
345
346 // in case of failure: throw exception
349
350 // otherwise exit as normal
351}
352
353#endif // DOXYGEN
354
356
357#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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define DeclException0(Exception0)
Definition exceptions.h:471
static ::ExceptionBase & ExcPreconditionerNotDefinite()
#define Assert(cond, exc)
#define AssertThrow(cond, exc)
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)
T signaling_nan()
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
void swap(SmartPointer< T, P > &t1, SmartPointer< T, Q > &t2)