1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2020 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_precondition_selector_h
17 #define dealii_precondition_selector_h
18
19
20 #include <deal.II/base/config.h>
21
22 #include <deal.II/base/smartpointer.h>
23 #include <deal.II/base/subscriptor.h>
24
25 #include <string>
26
27 DEAL_II_NAMESPACE_OPEN
28
29 // Forward declarations
30 #ifndef DOXYGEN
31 template <class number>
32 class Vector;
33 template <class number>
34 class SparseMatrix;
35 #endif
36
37
38 /*! @addtogroup Preconditioners
39 *@{
40 */
41
42 /**
43 * Selects the preconditioner. The constructor of this class takes the name of
44 * the preconditioning and the damping parameter @p omega of the
45 * preconditioning and the @p use_matrix function takes the matrix that is
46 * used by the matrix-builtin precondition functions. Each time, the
47 * <tt>operator()</tt> function is called, this preselected preconditioner,
48 * this matrix and this @p omega is used for the preconditioning. This class
49 * is designed for being used as argument of the @p solve function of a @p
50 * Solver and it covers the selection of all matrix-builtin precondition
51 * functions. The selection of other preconditioners, like BlockSOR or ILU
52 * should be handled in derived classes by the user.
53 *
54 * <h3>Usage</h3> The simplest use of this class is the following:
55 * @code
56 * // generate a @p SolverControl and a @p VectorMemory
57 * SolverControl control;
58 * VectorMemory<Vector<double> > memory;
59 *
60 * // generate a solver
61 * SolverCG<SparseMatrix<double>, Vector<double> > solver(control, memory);
62 *
63 * // generate a @p PreconditionSelector
64 * PreconditionSelector<SparseMatrix<double>, Vector<double> >
65 * preconditioning("jacobi", 1.);
66 *
67 * // give a matrix whose diagonal entries are to be used for the
68 * // preconditioning. Generally the matrix of the linear equation system Ax=b.
69 * preconditioning.use_matrix(A);
70 *
71 * // call the @p solve function with this preconditioning as last argument
72 * solver.solve(A,x,b,preconditioning);
73 * @endcode
74 * The same example where also the @p SolverSelector class is used reads
75 * @code
76 * // generate a @p SolverControl and a @p VectorMemory
77 * SolverControl control;
78 * VectorMemory<Vector<double> > memory;
79 *
80 * // generate a @p SolverSelector that calls the @p SolverCG
81 * SolverSelector<SparseMatrix<double>, Vector<double> >
82 * solver_selector("cg", control, memory);
83 *
84 * // generate a @p PreconditionSelector
85 * PreconditionSelector<SparseMatrix<double>, Vector<double> >
86 * preconditioning("jacobi", 1.);
87 *
88 * preconditioning.use_matrix(A);
89 *
90 * solver_selector.solve(A,x,b,preconditioning);
91 * @endcode
92 * Now the use of the @p SolverSelector in combination with the @p
93 * PreconditionSelector allows the user to select both, the solver and the
94 * preconditioner, at the beginning of their program and each time the solver is
95 * started (that is several times e.g. in a nonlinear iteration) this
96 * preselected solver and preconditioner is called.
97 */
98 template <typename MatrixType = SparseMatrix<double>,
99 typename VectorType = dealii::Vector<double>>
100 class PreconditionSelector : public Subscriptor
101 {
102 public:
103 /**
104 * Declare type for container size.
105 */
106 using size_type = typename MatrixType::size_type;
107
108 /**
109 * Constructor. @p omega denotes the damping parameter of the
110 * preconditioning.
111 */
112 PreconditionSelector(const std::string & preconditioning,
113 const typename VectorType::value_type &omega = 1.);
114
115 /**
116 * Destructor.
117 */
118 virtual ~PreconditionSelector() override;
119
120 /**
121 * Takes the matrix that is needed for preconditionings that involves a
122 * matrix. e.g. for @p precondition_jacobi, <tt>~_sor</tt>, <tt>~_ssor</tt>.
123 */
124 void
125 use_matrix(const MatrixType &M);
126
127 /**
128 * Return the dimension of the codomain (or range) space. Note that the
129 * matrix is of dimension $m \times n$.
130 */
131 size_type
132 m() const;
133
134 /**
135 * Return the dimension of the domain space. Note that the matrix is of
136 * dimension $m \times n$.
137 */
138 size_type
139 n() const;
140
141 /**
142 * Precondition procedure. Calls the preconditioning that was specified in
143 * the constructor.
144 */
145 virtual void
146 vmult(VectorType &dst, const VectorType &src) const;
147
148 /**
149 * Transpose precondition procedure. Calls the preconditioning that was
150 * specified in the constructor.
151 */
152 virtual void
153 Tvmult(VectorType &dst, const VectorType &src) const;
154
155 /**
156 * Get the names of all implemented preconditionings. The list of possible
157 * options includes:
158 * <ul>
159 * <li> "none" </li>
160 * <li> "jacobi" </li>
161 * <li> "sor" </li>
162 * <li> "ssor" </li>
163 * </ul>
164 */
165 static std::string
166 get_precondition_names();
167
168 /**
169 * @addtogroup Exceptions
170 * @{
171 */
172
173
174 /**
175 * Exception.
176 */
177 DeclException0(ExcNoMatrixGivenToUse);
178
179 //@}
180 protected:
181 /**
182 * Stores the name of the preconditioning.
183 */
184 std::string preconditioning;
185
186 private:
187 /**
188 * Matrix that is used for the matrix-builtin preconditioning function. cf.
189 * also @p PreconditionUseMatrix.
190 */
191 SmartPointer<const MatrixType, PreconditionSelector<MatrixType, VectorType>>
192 A;
193
194 /**
195 * Stores the damping parameter of the preconditioner.
196 */
197 const typename VectorType::value_type omega;
198 };
199
200 /*@}*/
201 /* --------------------- Inline and template functions ------------------- */
202
203
204 template <typename MatrixType, typename VectorType>
PreconditionSelector(const std::string & preconditioning,const typename VectorType::value_type & omega)205 PreconditionSelector<MatrixType, VectorType>::PreconditionSelector(
206 const std::string & preconditioning,
207 const typename VectorType::value_type &omega)
208 : preconditioning(preconditioning)
209 , omega(omega)
210 {}
211
212
213 template <typename MatrixType, typename VectorType>
~PreconditionSelector()214 PreconditionSelector<MatrixType, VectorType>::~PreconditionSelector()
215 {
216 // release the matrix A
217 A = nullptr;
218 }
219
220
221 template <typename MatrixType, typename VectorType>
222 void
use_matrix(const MatrixType & M)223 PreconditionSelector<MatrixType, VectorType>::use_matrix(const MatrixType &M)
224 {
225 A = &M;
226 }
227
228
229 template <typename MatrixType, typename VectorType>
230 inline typename PreconditionSelector<MatrixType, VectorType>::size_type
m()231 PreconditionSelector<MatrixType, VectorType>::m() const
232 {
233 Assert(A != nullptr, ExcNoMatrixGivenToUse());
234 return A->m();
235 }
236
237
238 template <typename MatrixType, typename VectorType>
239 inline typename PreconditionSelector<MatrixType, VectorType>::size_type
n()240 PreconditionSelector<MatrixType, VectorType>::n() const
241 {
242 Assert(A != nullptr, ExcNoMatrixGivenToUse());
243 return A->n();
244 }
245
246
247
248 template <typename MatrixType, typename VectorType>
249 void
vmult(VectorType & dst,const VectorType & src)250 PreconditionSelector<MatrixType, VectorType>::vmult(VectorType & dst,
251 const VectorType &src) const
252 {
253 if (preconditioning == "none")
254 {
255 dst = src;
256 }
257 else
258 {
259 Assert(A != nullptr, ExcNoMatrixGivenToUse());
260
261 if (preconditioning == "jacobi")
262 {
263 A->precondition_Jacobi(dst, src, omega);
264 }
265 else if (preconditioning == "sor")
266 {
267 A->precondition_SOR(dst, src, omega);
268 }
269 else if (preconditioning == "ssor")
270 {
271 A->precondition_SSOR(dst, src, omega);
272 }
273 else
274 Assert(false, ExcNotImplemented());
275 }
276 }
277
278
279 template <typename MatrixType, typename VectorType>
280 void
Tvmult(VectorType & dst,const VectorType & src)281 PreconditionSelector<MatrixType, VectorType>::Tvmult(
282 VectorType & dst,
283 const VectorType &src) const
284 {
285 if (preconditioning == "none")
286 {
287 dst = src;
288 }
289 else
290 {
291 Assert(A != nullptr, ExcNoMatrixGivenToUse());
292
293 if (preconditioning == "jacobi")
294 {
295 A->precondition_Jacobi(dst, src, omega); // Symmetric operation
296 }
297 else if (preconditioning == "sor")
298 {
299 A->precondition_TSOR(dst, src, omega);
300 }
301 else if (preconditioning == "ssor")
302 {
303 A->precondition_SSOR(dst, src, omega); // Symmetric operation
304 }
305 else
306 Assert(false, ExcNotImplemented());
307 }
308 }
309
310
311 template <typename MatrixType, typename VectorType>
312 std::string
get_precondition_names()313 PreconditionSelector<MatrixType, VectorType>::get_precondition_names()
314 {
315 return "none|jacobi|sor|ssor";
316 }
317
318
319 DEAL_II_NAMESPACE_CLOSE
320
321 #endif
322