Reference documentation for deal.II version GIT 9d46a3bd8c 2023-10-03 12:10: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\}}\)
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 
22 #include <deal.II/base/logstream.h>
25 
26 #include <deal.II/lac/solver.h>
28 
29 #include <cmath>
30 
32 
69 template <typename VectorType = Vector<double>>
70 class SolverMinRes : public SolverBase<VectorType>
71 {
72 public:
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 
120 protected:
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 
152 template <typename VectorType>
155  const AdditionalData &)
156  : SolverBase<VectorType>(cn, mem)
157  , res2(numbers::signaling_nan<double>())
158 {}
159 
160 
161 
162 template <typename VectorType>
164  const AdditionalData &)
165  : SolverBase<VectorType>(cn)
166  , res2(numbers::signaling_nan<double>())
167 {}
168 
169 
170 
171 template <typename VectorType>
172 double
174 {
175  return res2;
176 }
177 
178 
179 template <typename VectorType>
180 void
181 SolverMinRes<VectorType>::print_vectors(const unsigned int,
182  const VectorType &,
183  const VectorType &,
184  const VectorType &) const
185 {}
186 
187 
188 
189 template <typename VectorType>
190 template <typename MatrixType, typename PreconditionerType>
191 void
192 SolverMinRes<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
void swap(BlockIndices &u, BlockIndices &v)
@ 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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define DeclException0(Exception0)
Definition: exceptions.h:467
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcPreconditionerNotDefinite()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
Expression fabs(const Expression &x)
static const char A
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
long double gamma(const unsigned int n)
T signaling_nan()
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)