Reference documentation for deal.II version Git e038590ef8 2020-10-28 14:38:53 +0100
\(\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 - 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>
25 
26 #include <deal.II/lac/solver.h>
28 
29 #include <cmath>
30 
32 
35 
68 template <class VectorType = Vector<double>>
69 class SolverMinRes : public SolverBase<VectorType>
70 {
71 public:
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 
118 
119 protected:
123  virtual double
124  criterion();
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 
151 template <class VectorType>
154  const AdditionalData &)
155  : SolverBase<VectorType>(cn, mem)
157 {}
158 
159 
160 
161 template <class VectorType>
163  const AdditionalData &)
165  , res2(numbers::signaling_nan<double>())
166 {}
167 
168 
169 
170 template <class VectorType>
171 double
173 {
174  return res2;
175 }
176 
177 
178 template <class VectorType>
179 void
180 SolverMinRes<VectorType>::print_vectors(const unsigned int,
181  const VectorType &,
182  const VectorType &,
183  const VectorType &) const
184 {}
185 
186 
187 
188 template <class VectorType>
189 template <typename MatrixType, typename PreconditionerType>
190 void
191 SolverMinRes<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
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)
boost::signals2::signal< SolverControl::State(const unsigned int iteration, const double check_value, const VectorType &current_iterate), StateCombiner > iteration_status
Definition: solver.h:471
#define AssertThrow(cond, exc)
Definition: exceptions.h:1533
Stop iteration, goal reached.
#define Assert(cond, exc)
Definition: exceptions.h:1423
#define DeclException0(Exception0)
Definition: exceptions.h:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:369
virtual ~SolverMinRes() override=default
Expression fabs(const Expression &x)
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)
static const char A
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
T signaling_nan()
virtual double criterion()
long double gamma(const unsigned int n)
static ::ExceptionBase & ExcPreconditionerNotDefinite()
VectorMemory< VectorType > & memory
Definition: solver.h:420
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const