AMDiS  0.3
The Adaptive Multi-Dimensional Simulation Toolbox
ITL_Solver.hpp
1 #pragma once
2 
3 #include <string>
4 
5 // MTL4 includes
6 #include <boost/numeric/itl/itl.hpp>
7 
8 // more solvers defined in AMDiS
9 #include <amdis/linearalgebra/mtl/itl/minres.hpp>
10 #include <amdis/linearalgebra/mtl/itl/gcr.hpp>
11 #include <amdis/linearalgebra/mtl/itl/fgmres.hpp>
12 #include <amdis/linearalgebra/mtl/itl/fgmres_householder.hpp>
13 #include <amdis/linearalgebra/mtl/itl/gmres2.hpp>
14 #include <amdis/linearalgebra/mtl/itl/gmres_householder.hpp>
15 #include <amdis/linearalgebra/mtl/itl/preonly.hpp>
16 
17 #include <amdis/CreatorMap.hpp>
18 #include <amdis/Initfile.hpp>
19 #include <amdis/linearalgebra/LinearSolver.hpp>
20 #include <amdis/linearalgebra/mtl/KrylovRunner.hpp>
21 #include <amdis/linearalgebra/mtl/Traits.hpp>
22 #include <amdis/linearalgebra/mtl/UmfpackRunner.hpp>
23 
24 namespace AMDiS
25 {
36  {
37  public:
38  explicit cg_solver_type(std::string const& /*name*/) {}
39  template <class LinOp, class X, class B, class P, class I>
40  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
41  {
42  return itl::cg(A, x, b, precon, iter);
43  }
44  };
45 
46 
56  {
57  public:
58  explicit cgs_solver_type(std::string const& /*name*/) {}
59  template <class LinOp, class X, class B, class P, class I>
60  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
61  {
62  return itl::cgs(A, x, b, precon, iter);
63  }
64  };
65 
66 
76  {
77  public:
78  explicit bicg_solver_type(std::string const& /*name*/) {}
79  template <class LinOp, class X, class B, class P, class I>
80  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
81  {
82  return itl::bicg(A, x, b, precon, iter);
83  }
84  };
85 
86 
96  {
97  public:
98  explicit bicgstab_type(std::string const& /*name*/) {}
99  template <class LinOp, class X, class B, class P, class I>
100  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
101  {
102  return itl::bicgstab(A, x, b, precon, iter);
103  }
104  };
105 
106 
116  {
117  public:
118  explicit bicgstab2_type(std::string const& /*name*/) {}
119  template <class LinOp, class X, class B, class P, class I>
120  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
121  {
122  return itl::bicgstab_2(A, x, b, precon, iter);
123  }
124  };
125 
126 
135  {
136  public:
137  explicit qmr_solver_type(std::string const& /*name*/) {}
138  template <class LinOp, class X, class B, class P, class I>
139  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
140  {
141  itl::pc::identity<LinOp> id(A);
142  return itl::qmr(A, x, b, precon, id, iter);
143  }
144  };
145 
146 
156  {
157  public:
158  explicit tfqmr_solver_type(std::string const& /*name*/) {}
159  template <class LinOp, class X, class B, class P, class I>
160  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
161  {
162  itl::pc::identity<LinOp> id(A);
163  return itl::tfqmr(A, x, b, precon, id, iter);
164  }
165  };
166 
167 
177  {
178  int ell;
179  public:
180  explicit bicgstab_ell_type(std::string const& name) : ell(3)
181  {
182  Parameters::get(name + "->ell", ell);
183  }
184  template <class LinOp, class X, class B, class P, class I>
185  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
186  {
187  itl::pc::identity<LinOp> id(A);
188  return itl::bicgstab_ell(A, x, b, precon, id, iter, ell);
189  }
190  };
191 
192 
203  {
204 
205  public:
206  explicit gmres_type(std::string const& name) : restart(30), ortho(GRAM_SCHMIDT)
207  {
208  Parameters::get(name + "->restart", restart);
209  Parameters::get(name + "->orthogonalization", ortho);
210  }
211  template <class LinOp, class X, class B, class P, class I>
212  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
213  {
214  itl::pc::identity<LinOp> id(A);
215  switch (ortho)
216  {
217  default:
218  case GRAM_SCHMIDT:
219  return itl::gmres2(A, x, b, id, precon, iter, restart);
220  break;
221  case HOUSEHOLDER:
222  return itl::gmres_householder(A, x, b, precon, iter, restart);
223  break;
224  }
225  }
226 
227  private:
228  static constexpr int GRAM_SCHMIDT = 1;
229  static constexpr int HOUSEHOLDER = 2;
230 
231  int restart;
232  int ortho;
233  };
234 
235 
251  {
252  int s;
253  public:
254  explicit idr_s_type(std::string const& name) : s(30)
255  {
256  Parameters::get(name + "->s", s);
257  }
258  template <class LinOp, class X, class B, class P, class I>
259  int operator()(LinOp const& A, X& x, B const& b, P const&, I& iter)
260  {
261  itl::pc::identity<LinOp> id(A);
262  return itl::idr_s(A, x, b, id, id, iter, s);
263  }
264  };
265 
266 
276  {
277  public:
278  explicit minres_solver_type(std::string const& /*name*/) {}
279  template <class LinOp, class X, class B, class P, class I>
280  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
281  {
282  return itl::minres(A, x, b, precon, iter);
283  }
284  };
285 
286 
296  class gcr_type
297  {
298  int restart;
299 
300  public:
301  explicit gcr_type(std::string const& name) : restart(30)
302  {
303  Parameters::get(name + "->restart", restart);
304  }
305  template <class LinOp, class X, class B, class P, class I>
306  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
307  {
308  itl::pc::identity<LinOp> id(A);
309  return itl::gcr(A, x, b, id, precon, iter, restart);
310  }
311  };
312 
313 
325  {
326  public:
327  explicit fgmres_type(std::string const& name)
328  : restart(30), orthogonalization(GRAM_SCHMIDT)
329  {
330  Parameters::get(name + "->restart", restart);
331  Parameters::get(name + "->orthogonalization", orthogonalization);
332  }
333  template <class LinOp, class X, class B, class P, class I>
334  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
335  {
336  itl::pc::identity<LinOp> id(A);
337  switch (orthogonalization)
338  {
339  default:
340  case GRAM_SCHMIDT:
341  return itl::fgmres(A, x, b, id, precon, iter, restart);
342  break;
343  case HOUSEHOLDER:
344  return itl::fgmres_householder(A, x, b, precon, iter, restart);
345  break;
346  }
347  }
348 
349  private:
350  static constexpr int GRAM_SCHMIDT = 1;
351  static constexpr int HOUSEHOLDER = 2;
352 
353  int restart;
354  int orthogonalization;
355  };
356 
357 
366  {
367  public:
368  explicit preonly_type(std::string const& /*name*/) {}
369  template <class LinOp, class X, class B, class P, class I>
370  int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
371  {
372  return itl::preonly(A, x, b, precon, iter);
373  }
374  };
375 
376 
377  // ===========================================================================
378 
380 
398  template <class Mat, class Vec>
399  class DefaultCreators< LinearSolverInterface<Mat,Vec> >
400  {
401  using SolverBase = LinearSolverInterface<Mat,Vec>;
402 
403  template <class ITLSolver>
404  using SolverCreator
406 
407 #ifdef HAVE_UMFPACK
408  using UmfpackSolverCreator
409  = typename LinearSolver<Mat,Vec,UmfpackRunner<Mat,Vec>>::Creator;
410 #endif
411 
412  using Map = CreatorMap<SolverBase>;
413 
414  public:
415  static void init()
416  {
417  auto cg = new SolverCreator<cg_solver_type>;
418  Map::addCreator("cg", cg);
419 
420  auto cgs = new SolverCreator<cgs_solver_type>;
421  Map::addCreator("cgs", cgs);
422 
423  auto bicg = new SolverCreator<bicg_solver_type>;
424  Map::addCreator("bicg", bicg);
425 
426  auto bcgs = new SolverCreator<bicgstab_type>;
427  Map::addCreator("bicgstab", bcgs);
428  Map::addCreator("bcgs", bcgs);
429 
430  auto bcgs2 = new SolverCreator<bicgstab2_type>;
431  Map::addCreator("bicgstab2", bcgs2);
432  Map::addCreator("bcgs2", bcgs2);
433 
434  auto qmr = new SolverCreator<qmr_solver_type>;
435  Map::addCreator("qmr", qmr);
436 
437  auto tfqmr = new SolverCreator<tfqmr_solver_type>;
438  Map::addCreator("tfqmr", tfqmr);
439 
440  auto bcgsl = new SolverCreator<bicgstab_ell_type>;
441  Map::addCreator("bicgstab_ell", bcgsl);
442  Map::addCreator("bcgsl", bcgsl);
443 
444  auto gmres = new SolverCreator<gmres_type>;
445  Map::addCreator("gmres", gmres);
446 
447  auto fgmres = new SolverCreator<fgmres_type>;
448  Map::addCreator("fgmres", fgmres);
449 
450  auto minres = new SolverCreator<minres_solver_type>;
451  Map::addCreator("minres", minres);
452 
453  auto idrs = new SolverCreator<idr_s_type>;
454  Map::addCreator("idrs", idrs);
455 
456  auto gcr = new SolverCreator<gcr_type>;
457  Map::addCreator("gcr", gcr);
458 
459  auto preonly = new SolverCreator<preonly_type>;
460  Map::addCreator("preonly", preonly);
461 
462  // default iterative solver
463  Map::addCreator("default", gmres);
464 
465  init_direct(std::is_same<typename Dune::FieldTraits<typename Mat::value_type>::real_type, double>{});
466  }
467 
468  static void init_direct(std::false_type) {}
469  static void init_direct(std::true_type)
470  {
471 #ifdef HAVE_UMFPACK
472  auto umfpack = new UmfpackSolverCreator;
473  Map::addCreator("umfpack", umfpack);
474 
475  // default direct solvers
476  Map::addCreator("direct", umfpack);
477 #endif
478  }
479  };
480 
481 } // end namespace AMDiS
ITL_Solver <fgmres_type> implementation of flexible GMRes method.
Definition: ITL_Solver.hpp:324
ITL_Solver <qmr_solver_type> implementation of Quasi-Minimal Residual method.
Definition: ITL_Solver.hpp:134
ITL_Solver <bicgstab_type> implementation of stabilized bi-conjugate gradient method.
Definition: ITL_Solver.hpp:95
ITL_Solver <tfqmr_solver_type> implementation of Transposed-Free Quasi-Minimal Residual method...
Definition: ITL_Solver.hpp:155
ITL_Solver <gmres_type> implementation of generalized minimal residual method.
Definition: ITL_Solver.hpp:202
ITL_Solver <bicgstab2_type> implementation of BiCGStab(l) method with l=2.
Definition: ITL_Solver.hpp:115
A CreatorMap is used to construct objects, which types depends on key words determined at run time...
Definition: CreatorMap.hpp:29
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
static std::optional< T > get(std::string const &key)
Get parameter-values from parameter-tree.
Definition: Initfile.hpp:25
Abstract base class for linear solvers.
Definition: LinearSolverInterface.hpp:27
ITL_Solver <gcr_type> implementation of generalized conjugate residual method.
Definition: ITL_Solver.hpp:296
Wrapper class for various solvers. These solvers are parametrized by MatrixType and VectorType...
Definition: LinearSolver.hpp:25
Definition: CreatorMap.hpp:16
ITL_Solver <cg_solver_type> implementation of conjugate gradient method.
Definition: ITL_Solver.hpp:35
ITL_Solver <minres_solver_type> implementation of minimal residual method.
Definition: ITL_Solver.hpp:275
ITL_Solver <bicgstab_ell_type> implementation of stabilized BiCG(ell) method.
Definition: ITL_Solver.hpp:176
ITL_Solver <idr_s_type> implementation of Induced Dimension Reduction method.
Definition: ITL_Solver.hpp:250
ITL_Solver <bicg_solver_type> implementation of bi-conjugate gradient method.
Definition: ITL_Solver.hpp:75
void init(int &argc, char **&argv, std::string const &initFileName="")
Initialized the Environment for MPI.
Definition: AMDiS.hpp:29
ITL_Solver <preonly_type> implementation of preconditioner as.
Definition: ITL_Solver.hpp:365
ITL_Solver <cgs_solver_type> implementation of squared conjugate gradient method. ...
Definition: ITL_Solver.hpp:55
auto X()
Generator for CoordsFunction.
Definition: CoordsGridFunction.hpp:123