Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2000 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_solver_minres_h
17#define dealii_solver_minres_h
18
19
20#include <deal.II/base/config.h>
21
25
26#include <deal.II/lac/solver.h>
28
29#include <cmath>
30
32
69template <class VectorType = Vector<double>>
70class SolverMinRes : public SolverBase<VectorType>
71{
72public:
78 {};
79
85 const AdditionalData & data = AdditionalData());
86
92 const AdditionalData &data = AdditionalData());
93
97 virtual ~SolverMinRes() override = default;
98
102 template <typename MatrixType, typename PreconditionerType>
103 void
104 solve(const MatrixType & A,
105 VectorType & x,
106 const VectorType & b,
107 const PreconditionerType &preconditioner);
108
120protected:
124 virtual double
126
132 virtual void
133 print_vectors(const unsigned int step,
134 const VectorType & x,
135 const VectorType & r,
136 const VectorType & d) const;
137
144 double res2;
145};
146
148/*------------------------- Implementation ----------------------------*/
149
150#ifndef DOXYGEN
151
152template <class VectorType>
155 const AdditionalData &)
156 : SolverBase<VectorType>(cn, mem)
157 , res2(numbers::signaling_nan<double>())
158{}
159
160
161
162template <class VectorType>
164 const AdditionalData &)
165 : SolverBase<VectorType>(cn)
166 , res2(numbers::signaling_nan<double>())
167{}
168
169
170
171template <class VectorType>
172double
174{
175 return res2;
176}
177
178
179template <class VectorType>
180void
182 const VectorType &,
183 const VectorType &,
184 const VectorType &) const
185{}
186
187
188
189template <class VectorType>
190template <typename MatrixType, typename PreconditionerType>
191void
192SolverMinRes<VectorType>::solve(const MatrixType & A,
193 VectorType & x,
194 const VectorType & b,
195 const PreconditionerType &preconditioner)
196{
197 LogStream::Prefix prefix("minres");
198
199 // Memory allocation
200 typename VectorMemory<VectorType>::Pointer Vu0(this->memory);
201 typename VectorMemory<VectorType>::Pointer Vu1(this->memory);
202 typename VectorMemory<VectorType>::Pointer Vu2(this->memory);
203
204 typename VectorMemory<VectorType>::Pointer Vm0(this->memory);
205 typename VectorMemory<VectorType>::Pointer Vm1(this->memory);
206 typename VectorMemory<VectorType>::Pointer Vm2(this->memory);
207
208 typename VectorMemory<VectorType>::Pointer Vv(this->memory);
209
210 // define some aliases for simpler access
211 using vecptr = VectorType *;
212 vecptr u[3] = {Vu0.get(), Vu1.get(), Vu2.get()};
213 vecptr m[3] = {Vm0.get(), Vm1.get(), Vm2.get()};
214 VectorType &v = *Vv;
215
216 // resize the vectors, but do not set the values since they'd be overwritten
217 // soon anyway.
218 u[0]->reinit(b, true);
219 u[1]->reinit(b, true);
220 u[2]->reinit(b, true);
221 m[0]->reinit(b, true);
222 m[1]->reinit(b, true);
223 m[2]->reinit(b, true);
224 v.reinit(b, true);
225
226 // some values needed
227 double delta[3] = {0, 0, 0};
228 double f[2] = {0, 0};
229 double e[2] = {0, 0};
230
231 double r_l2 = 0;
232 double r0 = 0;
233 double tau = 0;
234 double c = 0;
235 double s = 0;
236 double d_ = 0;
237
238 // The iteration step.
239 unsigned int j = 1;
240
241
242 // Start of the solution process
243 A.vmult(*m[0], x);
244 *u[1] = b;
245 *u[1] -= *m[0];
246 // Precondition is applied.
247 // The preconditioner has to be
248 // positive definite and symmetric
249
250 // M v = u[1]
251 preconditioner.vmult(v, *u[1]);
252
253 delta[1] = v * (*u[1]);
254 // Preconditioner positive
255 Assert(delta[1] >= 0, ExcPreconditionerNotDefinite());
256
257 r0 = std::sqrt(delta[1]);
258 r_l2 = r0;
259
260
261 u[0]->reinit(b);
262 delta[0] = 1.;
263 m[0]->reinit(b);
264 m[1]->reinit(b);
265 m[2]->reinit(b);
266
267 SolverControl::State conv = this->iteration_status(0, r_l2, x);
268 while (conv == SolverControl::iterate)
269 {
270 if (delta[1] != 0)
271 v *= 1. / std::sqrt(delta[1]);
272 else
273 v.reinit(b);
274
275 A.vmult(*u[2], v);
276 u[2]->add(-std::sqrt(delta[1] / delta[0]), *u[0]);
277
278 const double gamma = *u[2] * v;
279 u[2]->add(-gamma / std::sqrt(delta[1]), *u[1]);
280 *m[0] = v;
281
282 // precondition: solve M v = u[2]
283 // Preconditioner has to be positive
284 // definite and symmetric.
285 preconditioner.vmult(v, *u[2]);
286
287 delta[2] = v * (*u[2]);
288
289 Assert(delta[2] >= 0, ExcPreconditionerNotDefinite());
290
291 if (j == 1)
292 {
293 d_ = gamma;
294 e[1] = std::sqrt(delta[2]);
295 }
296 if (j > 1)
297 {
298 d_ = s * e[0] - c * gamma;
299 e[0] = c * e[0] + s * gamma;
300 f[1] = s * std::sqrt(delta[2]);
301 e[1] = -c * std::sqrt(delta[2]);
302 }
303
304 const double d = std::sqrt(d_ * d_ + delta[2]);
305
306 if (j > 1)
307 tau *= s / c;
308 c = d_ / d;
309 tau *= c;
310
311 s = std::sqrt(delta[2]) / d;
312
313 if (j == 1)
314 tau = r0 * c;
315
316 m[0]->add(-e[0], *m[1]);
317 if (j > 1)
318 m[0]->add(-f[0], *m[2]);
319 *m[0] *= 1. / d;
320 x.add(tau, *m[0]);
321 r_l2 *= std::fabs(s);
322
323 conv = this->iteration_status(j, r_l2, x);
324
325 // next iteration step
326 ++j;
327 // All vectors have to be shifted
328 // one iteration step.
329 // This should be changed one time.
330 swap(*m[2], *m[1]);
331 swap(*m[1], *m[0]);
332
333 // likewise, but reverse direction:
334 // u[0] = u[1];
335 // u[1] = u[2];
336 swap(*u[0], *u[1]);
337 swap(*u[1], *u[2]);
338
339 // these are scalars, so need
340 // to bother
341 f[0] = f[1];
342 e[0] = e[1];
343 delta[0] = delta[1];
344 delta[1] = delta[2];
345 }
346
347 // in case of failure: throw exception
350
351 // otherwise exit as normal
352}
353
354#endif // DOXYGEN
355
357
358#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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define DeclException0(Exception0)
Definition exceptions.h:465
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)