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
33 /*!@addtogroup Solvers */
34 /*@{*/
35
36 /**
37 * Minimal residual method for symmetric matrices.
38 *
39 * For the requirements on matrices and vectors in order to work with this
40 * class, see the documentation of the Solver base class.
41 *
42 * Like all other solver classes, this class has a local structure called @p
43 * AdditionalData which is used to pass additional parameters to the solver,
44 * like damping parameters or the number of temporary vectors. We use this
45 * additional structure instead of passing these values directly to the
46 * constructor because this makes the use of the @p SolverSelector and other
47 * classes much easier and guarantees that these will continue to work even if
48 * number or type of the additional parameters for a certain solver changes.
49 *
50 * However, since the MinRes method does not need additional data, the
51 * respective structure is empty and does not offer any functionality. The
52 * constructor has a default argument, so you may call it without the
53 * additional parameter.
54 *
55 * The preconditioner has to be positive definite and symmetric
56 *
57 * The algorithm is taken from the Master thesis of Astrid Battermann with
58 * some changes. The full text can be found at
59 * http://scholar.lib.vt.edu/theses/public/etd-12164379662151/etd-title.html
60 *
61 *
62 * <h3>Observing the progress of linear solver iterations</h3>
63 *
64 * The solve() function of this class uses the mechanism described in the
65 * Solver base class to determine convergence. This mechanism can also be used
66 * to observe the progress of the iteration.
67 */
68 template <class VectorType = Vector<double>>
69 class SolverMinRes : public SolverBase<VectorType>
70 {
71 public:
72 /**
73 * Standardized data struct to pipe additional data to the solver. This
74 * solver does not need additional data yet.
75 */
76 struct AdditionalData
77 {};
78
79 /**
80 * Constructor.
81 */
82 SolverMinRes(SolverControl & cn,
83 VectorMemory<VectorType> &mem,
84 const AdditionalData & data = AdditionalData());
85
86 /**
87 * Constructor. Use an object of type GrowingVectorMemory as a default to
88 * allocate memory.
89 */
90 SolverMinRes(SolverControl & cn,
91 const AdditionalData &data = AdditionalData());
92
93 /**
94 * Virtual destructor.
95 */
96 virtual ~SolverMinRes() override = default;
97
98 /**
99 * Solve the linear system $Ax=b$ for x.
100 */
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
108 /**
109 * @addtogroup Exceptions
110 * @{
111 */
112
113 /**
114 * Exception
115 */
116 DeclException0(ExcPreconditionerNotDefinite);
117 //@}
118
119 protected:
120 /**
121 * Implementation of the computation of the norm of the residual.
122 */
123 virtual double
124 criterion();
125
126 /**
127 * Interface for derived class. This function gets the current iteration
128 * vector, the residual and the update vector in each step. It can be used
129 * for graphical output of the convergence history.
130 */
131 virtual void
132 print_vectors(const unsigned int step,
133 const VectorType & x,
134 const VectorType & r,
135 const VectorType & d) const;
136
137 /**
138 * Within the iteration loop, the square of the residual vector is stored in
139 * this variable. The function @p criterion uses this variable to compute
140 * the convergence value, which in this class is the norm of the residual
141 * vector and thus the square root of the @p res2 value.
142 */
143 double res2;
144 };
145
146 /*@}*/
147 /*------------------------- Implementation ----------------------------*/
148
149 #ifndef DOXYGEN
150
151 template <class VectorType>
SolverMinRes(SolverControl & cn,VectorMemory<VectorType> & mem,const AdditionalData &)152 SolverMinRes<VectorType>::SolverMinRes(SolverControl & cn,
153 VectorMemory<VectorType> &mem,
154 const AdditionalData &)
155 : SolverBase<VectorType>(cn, mem)
156 , res2(numbers::signaling_nan<double>())
157 {}
158
159
160
161 template <class VectorType>
SolverMinRes(SolverControl & cn,const AdditionalData &)162 SolverMinRes<VectorType>::SolverMinRes(SolverControl &cn,
163 const AdditionalData &)
164 : SolverBase<VectorType>(cn)
165 , res2(numbers::signaling_nan<double>())
166 {}
167
168
169
170 template <class VectorType>
171 double
criterion()172 SolverMinRes<VectorType>::criterion()
173 {
174 return res2;
175 }
176
177
178 template <class VectorType>
179 void
print_vectors(const unsigned int,const VectorType &,const VectorType &,const VectorType &)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
solve(const MatrixType & A,VectorType & x,const VectorType & b,const PreconditionerType & preconditioner)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
347 AssertThrow(conv == SolverControl::success,
348 SolverControl::NoConvergence(j, r_l2));
349
350 // otherwise exit as normal
351 }
352
353 #endif // DOXYGEN
354
355 DEAL_II_NAMESPACE_CLOSE
356
357 #endif
358