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_solver_selector_h
17 #define dealii_solver_selector_h
18
19
20 #include <deal.II/base/config.h>
21
22 #include <deal.II/base/smartpointer.h>
23
24 #include <deal.II/lac/precondition.h>
25 #include <deal.II/lac/solver.h>
26 #include <deal.II/lac/solver_bicgstab.h>
27 #include <deal.II/lac/solver_cg.h>
28 #include <deal.II/lac/solver_control.h>
29 #include <deal.II/lac/solver_gmres.h>
30 #include <deal.II/lac/solver_minres.h>
31 #include <deal.II/lac/solver_richardson.h>
32 #include <deal.II/lac/vector.h>
33
34 DEAL_II_NAMESPACE_OPEN
35
36
37 /*!@addtogroup Solvers */
38 /*@{*/
39
40 /**
41 * Selects a solver by changing a parameter.
42 *
43 * By calling the @p solve function of this @p SolverSelector, it selects the
44 * @p solve function of that @p Solver that was specified in the constructor
45 * of this class.
46 *
47 * <h3>Usage</h3> The simplest use of this class is the following:
48 * @code
49 * // generate a @p SolverSelector that calls the @p SolverCG
50 * SolverControl control;
51 * SolverSelector<Vector<double> > solver_selector ("cg", control);
52 *
53 * // generate e.g. a @p PreconditionRelaxation
54 * PreconditionRelaxation<SparseMatrix<double>, Vector<double> >
55 * preconditioning(
56 * A, &SparseMatrix<double>::template precondition_SSOR<double>,0.8);
57 *
58 * // call the @p solve function with this preconditioning as last argument
59 * solver_selector.solve(A,x,b,preconditioning);
60 * @endcode
61 * But the full usefulness of the @p SolverSelector class is not clear until
62 * the presentation of the following example that assumes the user using the
63 * @p ParameterHandler class and having declared a "solver" entry, e.g. with
64 * @code
65 * Parameter_Handler prm;
66 * prm.declare_entry ("solver", "none",
67 * Patterns::Selection(
68 * SolverSelector<>::get_solver_names()));
69 * ...
70 * @endcode
71 * Assuming that in the users parameter file there exists the line
72 * @code
73 * set solver = cg
74 * @endcode
75 * then the constructor call in the above example can be written as
76 * @code
77 * SolverSelector<SparseMatrix<double>, Vector<double> >
78 * solver_selector(prm.get("solver"), control);
79 * @endcode
80 *
81 *
82 * If at some time there exists a new solver "xyz" then the user does not need
83 * to change their program. Only in the implementation of the @p SolverSelector
84 * the calling of this solver has to be added and each user with program lines
85 * quoted above only needs to 'set solver = xyz' in their parameter file to get
86 * access to that new solver.
87 */
88 template <typename VectorType = Vector<double>>
89 class SolverSelector : public Subscriptor
90 {
91 public:
92 /**
93 * An alias for the underlying vector type
94 */
95 using vector_type = VectorType;
96
97 /**
98 * Constructor, filling in default values
99 */
100 SolverSelector() = default;
101
102 /**
103 * Constructor, selecting the solver @p name
104 * and the SolverControl object @p control already.
105 */
106 SolverSelector(const std::string &name, SolverControl &control);
107
108 /**
109 * Destructor
110 */
111 virtual ~SolverSelector() override = default;
112
113 /**
114 * Solver procedure. Calls the @p solve function of the @p solver whose @p
115 * SolverName was specified in the constructor.
116 */
117 template <class Matrix, class Preconditioner>
118 void
119 solve(const Matrix & A,
120 VectorType & x,
121 const VectorType & b,
122 const Preconditioner &precond) const;
123
124 /**
125 * Select a new solver. Note that all solver names used in this class are
126 * all lower case.
127 */
128 void
129 select(const std::string &name);
130
131 /**
132 * Set a new SolverControl. This needs to be set before solving.
133 */
134 void
135 set_control(SolverControl &ctrl);
136
137 /**
138 * Set the additional data. For more information see the @p Solver class.
139 */
140 void
141 set_data(const typename SolverRichardson<VectorType>::AdditionalData &data);
142
143 /**
144 * Set the additional data. For more information see the @p Solver class.
145 */
146 void
147 set_data(const typename SolverCG<VectorType>::AdditionalData &data);
148
149 /**
150 * Set the additional data. For more information see the @p Solver class.
151 */
152 void
153 set_data(const typename SolverMinRes<VectorType>::AdditionalData &data);
154
155 /**
156 * Set the additional data. For more information see the @p Solver class.
157 */
158 void
159 set_data(const typename SolverBicgstab<VectorType>::AdditionalData &data);
160
161 /**
162 * Set the additional data. For more information see the @p Solver class.
163 */
164 void
165 set_data(const typename SolverGMRES<VectorType>::AdditionalData &data);
166
167 /**
168 * Set the additional data. For more information see the @p Solver class.
169 */
170 void
171 set_data(const typename SolverFGMRES<VectorType>::AdditionalData &data);
172
173 /**
174 * Get the names of all implemented solvers. The list of possible
175 * options includes:
176 * <ul>
177 * <li> "richardson" </li>
178 * <li> "cg" </li>
179 * <li> "bicgstab" </li>
180 * <li> "gmres" </li>
181 * <li> "fgmres" </li>
182 * <li> "minres" </li>
183 * </ul>
184 */
185 static std::string
186 get_solver_names();
187
188 /**
189 * Exception.
190 */
191 DeclException1(ExcSolverDoesNotExist,
192 std::string,
193 << "Solver " << arg1 << " does not exist. Use one of "
194 << std::endl
195 << get_solver_names());
196
197
198
199 protected:
200 /**
201 * Stores the @p SolverControl that is needed in the constructor of each @p
202 * Solver class. This can be changed with @p set_control().
203 */
204 SmartPointer<SolverControl, SolverSelector<VectorType>> control;
205
206 /**
207 * Stores the name of the solver.
208 */
209 std::string solver_name;
210
211 private:
212 /**
213 * Stores the additional data.
214 */
215 typename SolverRichardson<VectorType>::AdditionalData richardson_data;
216
217 /**
218 * Stores the additional data.
219 */
220 typename SolverCG<VectorType>::AdditionalData cg_data;
221
222 /**
223 * Stores the additional data.
224 */
225 typename SolverMinRes<VectorType>::AdditionalData minres_data;
226
227 /**
228 * Stores the additional data.
229 */
230 typename SolverBicgstab<VectorType>::AdditionalData bicgstab_data;
231
232 /**
233 * Stores the additional data.
234 */
235 typename SolverGMRES<VectorType>::AdditionalData gmres_data;
236
237 /**
238 * Stores the additional data.
239 */
240 typename SolverFGMRES<VectorType>::AdditionalData fgmres_data;
241 };
242
243 /*@}*/
244 /* --------------------- Inline and template functions ------------------- */
245
246
247 template <typename VectorType>
SolverSelector(const std::string & name,SolverControl & solver_control)248 SolverSelector<VectorType>::SolverSelector(const std::string &name,
249 SolverControl & solver_control)
250 : solver_name(name)
251 , control(&solver_control)
252 {}
253
254
255
256 template <typename VectorType>
257 void
select(const std::string & name)258 SolverSelector<VectorType>::select(const std::string &name)
259 {
260 solver_name = name;
261 }
262
263
264
265 template <typename VectorType>
266 template <class Matrix, class Preconditioner>
267 void
solve(const Matrix & A,VectorType & x,const VectorType & b,const Preconditioner & precond)268 SolverSelector<VectorType>::solve(const Matrix & A,
269 VectorType & x,
270 const VectorType & b,
271 const Preconditioner &precond) const
272 {
273 if (solver_name == "richardson")
274 {
275 SolverRichardson<VectorType> solver(*control, richardson_data);
276 solver.solve(A, x, b, precond);
277 }
278 else if (solver_name == "cg")
279 {
280 SolverCG<VectorType> solver(*control, cg_data);
281 solver.solve(A, x, b, precond);
282 }
283 else if (solver_name == "minres")
284 {
285 SolverMinRes<VectorType> solver(*control, minres_data);
286 solver.solve(A, x, b, precond);
287 }
288 else if (solver_name == "bicgstab")
289 {
290 SolverBicgstab<VectorType> solver(*control, bicgstab_data);
291 solver.solve(A, x, b, precond);
292 }
293 else if (solver_name == "gmres")
294 {
295 SolverGMRES<VectorType> solver(*control, gmres_data);
296 solver.solve(A, x, b, precond);
297 }
298 else if (solver_name == "fgmres")
299 {
300 SolverFGMRES<VectorType> solver(*control, fgmres_data);
301 solver.solve(A, x, b, precond);
302 }
303 else
304 Assert(false, ExcSolverDoesNotExist(solver_name));
305 }
306
307
308
309 template <typename VectorType>
310 void
set_control(SolverControl & ctrl)311 SolverSelector<VectorType>::set_control(SolverControl &ctrl)
312 {
313 control = &ctrl;
314 }
315
316
317
318 template <typename VectorType>
319 std::string
get_solver_names()320 SolverSelector<VectorType>::get_solver_names()
321 {
322 return "richardson|cg|bicgstab|gmres|fgmres|minres";
323 }
324
325
326
327 template <typename VectorType>
328 void
set_data(const typename SolverGMRES<VectorType>::AdditionalData & data)329 SolverSelector<VectorType>::set_data(
330 const typename SolverGMRES<VectorType>::AdditionalData &data)
331 {
332 gmres_data = data;
333 }
334
335
336
337 template <typename VectorType>
338 void
set_data(const typename SolverFGMRES<VectorType>::AdditionalData & data)339 SolverSelector<VectorType>::set_data(
340 const typename SolverFGMRES<VectorType>::AdditionalData &data)
341 {
342 fgmres_data = data;
343 }
344
345
346
347 template <typename VectorType>
348 void
set_data(const typename SolverRichardson<VectorType>::AdditionalData & data)349 SolverSelector<VectorType>::set_data(
350 const typename SolverRichardson<VectorType>::AdditionalData &data)
351 {
352 richardson_data = data;
353 }
354
355
356
357 template <typename VectorType>
358 void
set_data(const typename SolverCG<VectorType>::AdditionalData & data)359 SolverSelector<VectorType>::set_data(
360 const typename SolverCG<VectorType>::AdditionalData &data)
361 {
362 cg_data = data;
363 }
364
365
366
367 template <typename VectorType>
368 void
set_data(const typename SolverMinRes<VectorType>::AdditionalData & data)369 SolverSelector<VectorType>::set_data(
370 const typename SolverMinRes<VectorType>::AdditionalData &data)
371 {
372 minres_data = data;
373 }
374
375
376
377 template <typename VectorType>
378 void
set_data(const typename SolverBicgstab<VectorType>::AdditionalData & data)379 SolverSelector<VectorType>::set_data(
380 const typename SolverBicgstab<VectorType>::AdditionalData &data)
381 {
382 bicgstab_data = data;
383 }
384
385 DEAL_II_NAMESPACE_CLOSE
386
387 #endif
388