Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
solver_minres.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2019 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>
23 #include <deal.II/base/signaling_nan.h>
24 #include <deal.II/base/subscriptor.h>
25 
26 #include <deal.II/lac/solver.h>
27 #include <deal.II/lac/solver_control.h>
28 
29 #include <cmath>
30 
31 DEAL_II_NAMESPACE_OPEN
32 
35 
71 template <class VectorType = Vector<double>>
72 class SolverMinRes : public SolverBase<VectorType>
73 {
74 public:
80  {};
81 
87  const AdditionalData & data = AdditionalData());
88 
94  const AdditionalData &data = AdditionalData());
95 
99  virtual ~SolverMinRes() override = default;
100 
104  template <typename MatrixType, typename PreconditionerType>
105  void
106  solve(const MatrixType & A,
107  VectorType & x,
108  const VectorType & b,
109  const PreconditionerType &preconditioner);
110 
121 
122 protected:
126  virtual double
127  criterion();
128 
134  virtual void
135  print_vectors(const unsigned int step,
136  const VectorType & x,
137  const VectorType & r,
138  const VectorType & d) const;
139 
146  double res2;
147 };
148 
150 /*------------------------- Implementation ----------------------------*/
151 
152 #ifndef DOXYGEN
153 
154 template <class VectorType>
157  const AdditionalData &)
158  : SolverBase<VectorType>(cn, mem)
159  , res2(numbers::signaling_nan<double>())
160 {}
161 
162 
163 
164 template <class VectorType>
166  const AdditionalData &)
167  : SolverBase<VectorType>(cn)
168  , res2(numbers::signaling_nan<double>())
169 {}
170 
171 
172 
173 template <class VectorType>
174 double
176 {
177  return res2;
178 }
179 
180 
181 template <class VectorType>
182 void
183 SolverMinRes<VectorType>::print_vectors(const unsigned int,
184  const VectorType &,
185  const VectorType &,
186  const VectorType &) const
187 {}
188 
189 
190 
191 template <class VectorType>
192 template <typename MatrixType, typename PreconditionerType>
193 void
194 SolverMinRes<VectorType>::solve(const MatrixType & A,
195  VectorType & x,
196  const VectorType & b,
197  const PreconditionerType &preconditioner)
198 {
199  LogStream::Prefix prefix("minres");
200 
201  // Memory allocation
202  typename VectorMemory<VectorType>::Pointer Vu0(this->memory);
203  typename VectorMemory<VectorType>::Pointer Vu1(this->memory);
204  typename VectorMemory<VectorType>::Pointer Vu2(this->memory);
205 
206  typename VectorMemory<VectorType>::Pointer Vm0(this->memory);
207  typename VectorMemory<VectorType>::Pointer Vm1(this->memory);
208  typename VectorMemory<VectorType>::Pointer Vm2(this->memory);
209 
210  typename VectorMemory<VectorType>::Pointer Vv(this->memory);
211 
212  // define some aliases for simpler access
213  using vecptr = VectorType *;
214  vecptr u[3] = {Vu0.get(), Vu1.get(), Vu2.get()};
215  vecptr m[3] = {Vm0.get(), Vm1.get(), Vm2.get()};
216  VectorType &v = *Vv;
217 
218  // resize the vectors, but do not set the values since they'd be overwritten
219  // soon anyway.
220  u[0]->reinit(b, true);
221  u[1]->reinit(b, true);
222  u[2]->reinit(b, true);
223  m[0]->reinit(b, true);
224  m[1]->reinit(b, true);
225  m[2]->reinit(b, true);
226  v.reinit(b, true);
227 
228  // some values needed
229  double delta[3] = {0, 0, 0};
230  double f[2] = {0, 0};
231  double e[2] = {0, 0};
232 
233  double r_l2 = 0;
234  double r0 = 0;
235  double tau = 0;
236  double c = 0;
237  double s = 0;
238  double d_ = 0;
239 
240  // The iteration step.
241  unsigned int j = 1;
242 
243 
244  // Start of the solution process
245  A.vmult(*m[0], x);
246  *u[1] = b;
247  *u[1] -= *m[0];
248  // Precondition is applied.
249  // The preconditioner has to be
250  // positive definite and symmetric
251 
252  // M v = u[1]
253  preconditioner.vmult(v, *u[1]);
254 
255  delta[1] = v * (*u[1]);
256  // Preconditioner positive
257  Assert(delta[1] >= 0, ExcPreconditionerNotDefinite());
258 
259  r0 = std::sqrt(delta[1]);
260  r_l2 = r0;
261 
262 
263  u[0]->reinit(b);
264  delta[0] = 1.;
265  m[0]->reinit(b);
266  m[1]->reinit(b);
267  m[2]->reinit(b);
268 
269  SolverControl::State conv = this->iteration_status(0, r_l2, x);
270  while (conv == SolverControl::iterate)
271  {
272  if (delta[1] != 0)
273  v *= 1. / std::sqrt(delta[1]);
274  else
275  v.reinit(b);
276 
277  A.vmult(*u[2], v);
278  u[2]->add(-std::sqrt(delta[1] / delta[0]), *u[0]);
279 
280  const double gamma = *u[2] * v;
281  u[2]->add(-gamma / std::sqrt(delta[1]), *u[1]);
282  *m[0] = v;
283 
284  // precondition: solve M v = u[2]
285  // Preconditioner has to be positive
286  // definite and symmetric.
287  preconditioner.vmult(v, *u[2]);
288 
289  delta[2] = v * (*u[2]);
290 
291  Assert(delta[2] >= 0, ExcPreconditionerNotDefinite());
292 
293  if (j == 1)
294  {
295  d_ = gamma;
296  e[1] = std::sqrt(delta[2]);
297  }
298  if (j > 1)
299  {
300  d_ = s * e[0] - c * gamma;
301  e[0] = c * e[0] + s * gamma;
302  f[1] = s * std::sqrt(delta[2]);
303  e[1] = -c * std::sqrt(delta[2]);
304  }
305 
306  const double d = std::sqrt(d_ * d_ + delta[2]);
307 
308  if (j > 1)
309  tau *= s / c;
310  c = d_ / d;
311  tau *= c;
312 
313  s = std::sqrt(delta[2]) / d;
314 
315  if (j == 1)
316  tau = r0 * c;
317 
318  m[0]->add(-e[0], *m[1]);
319  if (j > 1)
320  m[0]->add(-f[0], *m[2]);
321  *m[0] *= 1. / d;
322  x.add(tau, *m[0]);
323  r_l2 *= std::fabs(s);
324 
325  conv = this->iteration_status(j, r_l2, x);
326 
327  // next iteration step
328  ++j;
329  // All vectors have to be shifted
330  // one iteration step.
331  // This should be changed one time.
332  swap(*m[2], *m[1]);
333  swap(*m[1], *m[0]);
334 
335  // likewise, but reverse direction:
336  // u[0] = u[1];
337  // u[1] = u[2];
338  swap(*u[0], *u[1]);
339  swap(*u[1], *u[2]);
340 
341  // these are scalars, so need
342  // to bother
343  f[0] = f[1];
344  e[0] = e[1];
345  delta[0] = delta[1];
346  delta[1] = delta[2];
347  }
348 
349  // in case of failure: throw exception
352 
353  // otherwise exit as normal
354 }
355 
356 #endif // DOXYGEN
357 
358 DEAL_II_NAMESPACE_CLOSE
359 
360 #endif
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
Continue iteration.
SolverMinRes(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
void swap(BlockIndices &u, BlockIndices &v)
Stop iteration, goal reached.
#define Assert(cond, exc)
Definition: exceptions.h:1407
#define DeclException0(Exception0)
Definition: exceptions.h:473
virtual ~SolverMinRes() override=default
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
T signaling_nan()
virtual double criterion()
static ::ExceptionBase & ExcPreconditionerNotDefinite()
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const