ITL_Solver.hpp 14.3 KB
Newer Older
1
2
#pragma once

3
4
#include <string>

5
// MTL4 includes
6
#include <boost/numeric/itl/itl.hpp>
7
8

// more solvers defined in AMDiS
9
10
11
12
13
14
15
#include <amdis/linearalgebra/mtl/itl/minres.hpp>
#include <amdis/linearalgebra/mtl/itl/gcr.hpp>
#include <amdis/linearalgebra/mtl/itl/fgmres.hpp>
#include <amdis/linearalgebra/mtl/itl/fgmres_householder.hpp>
#include <amdis/linearalgebra/mtl/itl/gmres2.hpp>
#include <amdis/linearalgebra/mtl/itl/gmres_householder.hpp>
#include <amdis/linearalgebra/mtl/itl/preonly.hpp>
16

17
18
#include <amdis/CreatorMap.hpp>
#include <amdis/Initfile.hpp>
19
20
#include <amdis/linearalgebra/LinearSolver.hpp>
#include <amdis/linearalgebra/mtl/KrylovRunner.hpp>
21
#include <amdis/linearalgebra/mtl/Traits.hpp>
22
#include <amdis/linearalgebra/mtl/UmfpackRunner.hpp>
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

namespace AMDiS
{
  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref cg_solver_type> implementation of conjugate gradient
   * method \implements ITL_Solver
   *
   * Solves a linear system \f$ Ax=b \f$ by the conjugate gradient method (CG)
   * and can be used for symmetric positive definite system matrices.
   * Right preconditioner is ignored.
   */
  class cg_solver_type
  {
  public:
38
    explicit cg_solver_type(std::string const& /*name*/) {}
39
40
    template <class LinOp, class X, class B, class P, class I>
    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
47
48
49
50
51
52
53
54
55
56
57
    }
  };


  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref cgs_solver_type> implementation of squared conjugate
   * gradient method \implements ITL_Solver
   *
   * Solves a linear system \f$ Ax=b \f$ by the squared conjugate gradient method
   * (CGS). Right preconditioner is ignored.
   */
  class cgs_solver_type
  {
  public:
58
    explicit cgs_solver_type(std::string const& /*name*/) {}
59
60
    template <class LinOp, class X, class B, class P, class I>
    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
67
68
69
70
71
72
73
74
75
76
77
    }
  };


  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref bicg_solver_type> implementation of bi-conjugate
   * gradient method \implements ITL_Solver
   *
   * Solves a linear system \f$ Ax=b \f$ by a BiCG method and can be used for
   * system matrices.
   */
  class bicg_solver_type
  {
  public:
78
    explicit bicg_solver_type(std::string const& /*name*/) {}
79
80
    template <class LinOp, class X, class B, class P, class I>
    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
87
88
89
90
91
92
93
94
95
96
97
    }
  };


  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref bicgstab_type> implementation of stabilized
   * bi-conjugate gradient method \implements ITL_Solver
   *
   * Solves a linear system \f$ Ax=b \f$ by a stabilized BiCG method and can be
   * used for system matrices.
   */
  class bicgstab_type
  {
  public:
98
    explicit bicgstab_type(std::string const& /*name*/) {}
99
100
    template <class LinOp, class X, class B, class P, class I>
    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
107
108
109
110
111
112
113
114
115
116
117
    }
  };


  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref bicgstab2_type> implementation of BiCGStab(l) method
   * with l=2 \implements ITL_Solver
   *
   * Solves a linear system \f$ Ax=b \f$ by a stabilized BiCG(2) method and can
   * be used for system matrices.
   */
  class bicgstab2_type
  {
  public:
118
    explicit bicgstab2_type(std::string const& /*name*/) {}
119
120
    template <class LinOp, class X, class B, class P, class I>
    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
127
128
129
130
131
132
133
134
135
136
    }
  };


  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref qmr_solver_type> implementation of Quasi-Minimal
   * Residual method \implements ITL_Solver
   *
   * Solves a linear system \f$ Ax=b \f$ by the Quasi-Minimal Residual method (QMR).
   */
  class qmr_solver_type
  {
  public:
137
    explicit qmr_solver_type(std::string const& /*name*/) {}
138
139
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
140
    {
141
142
      itl::pc::identity<LinOp> id(A);
      return itl::qmr(A, x, b, precon, id, iter);
143
144
145
146
147
148
149
150
151
152
    }
  };


  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref tfqmr_solver_type> implementation of Transposed-Free
   * Quasi-Minimal Residual method \implements ITL_Solver
   *
   * Solves a linear system by the Transposed-Free Quasi-Minimal Residual method
153
   * (TFQMR).
154
155
156
157
   */
  class tfqmr_solver_type
  {
  public:
158
    explicit tfqmr_solver_type(std::string const& /*name*/) {}
159
160
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
161
    {
162
163
      itl::pc::identity<LinOp> id(A);
      return itl::tfqmr(A, x, b, precon, id, iter);
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
    }
  };


  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref bicgstab_ell_type> implementation of stabilized
   * BiCG(ell) method \implements ITL_Solver
   *
   * Solves a linear system by a stabilized BiCG(ell) method and can be used for
   * system matrices. The parameter ell [3] can be specified.
   */
  class bicgstab_ell_type
  {
    int ell;
  public:
180
    explicit bicgstab_ell_type(std::string const& name) : ell(3)
181
182
183
    {
      Parameters::get(name + "->ell", ell);
    }
184
185
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
186
    {
187
188
      itl::pc::identity<LinOp> id(A);
      return itl::bicgstab_ell(A, x, b, precon, id, iter, ell);
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
    }
  };


  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref gmres_type> implementation of generalized minimal
   * residual method \implements ITL_Solver
   *
   * Solves a linear system by the GMRES method.
   * The parameter restart [30] is the maximal number of orthogonalized vectors.
   * The method is not preconditioned
   */
  class gmres_type
  {

  public:
206
    explicit gmres_type(std::string const& name) : restart(30), ortho(GRAM_SCHMIDT)
207
208
209
210
    {
      Parameters::get(name + "->restart", restart);
      Parameters::get(name + "->orthogonalization", ortho);
    }
211
212
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
213
    {
214
      itl::pc::identity<LinOp> id(A);
215
216
217
218
      switch (ortho)
      {
      default:
      case GRAM_SCHMIDT:
219
        return itl::gmres2(A, x, b, id, precon, iter, restart);
220
221
        break;
      case HOUSEHOLDER:
222
        return itl::gmres_householder(A, x, b, precon, iter, restart);
223
224
225
        break;
      }
    }
226

227
228
229
  private:
    static constexpr int GRAM_SCHMIDT = 1;
    static constexpr int HOUSEHOLDER = 2;
230

231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
    int restart;
    int ortho;
  };


  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref idr_s_type> implementation of Induced Dimension
   * Reduction method \implements ITL_Solver
   *
   * Solves a linear system by an Induced Dimension Reduction method and can be
   * used for system matrices.  The parameter s [30] can be specified.
   *
   * Peter Sonneveld and Martin B. van Gijzen, IDR(s): a family of simple and fast
   * algorithms for solving large nonsymmetric linear systems.
   * SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035-1062 (2008). (copyright SIAM)
247
248
   *
   * NOTE: preconditioners are ignored.
249
250
251
252
253
   */
  class idr_s_type
  {
    int s;
  public:
254
    explicit idr_s_type(std::string const& name) : s(30)
255
256
257
    {
      Parameters::get(name + "->s", s);
    }
258
259
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const&, I& iter)
260
    {
261
262
      itl::pc::identity<LinOp> id(A);
      return itl::idr_s(A, x, b, id, id, iter, s);
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
    }
  };


  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref minres_solver_type> implementation of minimal
   * residual method \implements ITL_Solver
   *
   * Solves a linear system by the Minres method. Can be used for symmetric
   * indefinite systems.
   */
  class minres_solver_type
  {
  public:
278
    explicit minres_solver_type(std::string const& /*name*/) {}
279
280
    template <class LinOp, class X, class B, class P, class I>
    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
287
288
289
290
291
292
293
294
295
296
297
298
299
300
    }
  };


  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref gcr_type> implementation of generalized conjugate
   * residual method \implements ITL_Solver
   *
   * Solves a linear system by the GCR method - generalized conjugate residual
   * method. The parameter restart [30] is the maximal number of orthogonalized
   * vectors.
   */
  class gcr_type
  {
    int restart;

  public:
301
    explicit gcr_type(std::string const& name) : restart(30)
302
303
304
    {
      Parameters::get(name + "->restart", restart);
    }
305
306
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
307
    {
308
309
      itl::pc::identity<LinOp> id(A);
      return itl::gcr(A, x, b, id, precon, iter, restart);
310
311
    }
  };
312

313
314
315
316
317
318
319
320
321
322
323
324
325
326

  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref fgmres_type> implementation of flexible GMRes method
   * \implements ITL_Solver
   *
   * Solves a linear system by the FGMRES method.
   * The parameter restart [30] is the maximal number of orthogonalized vectors.
   * See reference "A Flexible Inner-Outer Preconditiones GMRES Algorithm",
   * Youcef Saad, (1993)
   */
  class fgmres_type
  {
  public:
Praetorius, Simon's avatar
Praetorius, Simon committed
327
328
    explicit fgmres_type(std::string const& name)
      : restart(30), orthogonalization(GRAM_SCHMIDT)
329
330
331
332
    {
      Parameters::get(name + "->restart", restart);
      Parameters::get(name + "->orthogonalization", orthogonalization);
    }
333
334
    template <class LinOp, class X, class B, class P, class I>
    int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
335
    {
336
      itl::pc::identity<LinOp> id(A);
337
338
339
340
      switch (orthogonalization)
      {
      default:
      case GRAM_SCHMIDT:
341
        return itl::fgmres(A, x, b, id, precon, iter, restart);
342
343
        break;
      case HOUSEHOLDER:
344
        return itl::fgmres_householder(A, x, b, precon, iter, restart);
345
346
347
        break;
      }
    }
348

349
350
351
  private:
    static constexpr int GRAM_SCHMIDT = 1;
    static constexpr int HOUSEHOLDER = 2;
352

353
354
355
356
    int restart;
    int orthogonalization;
  };

357

358
359
360
361
362
363
364
365
366
367
  /**
   * \ingroup Solver
   * \brief ITL_Solver <\ref preonly_type> implementation of preconditioner as
   * \implements ITL_Solver
   *
   * Solves a linear system by applying a preconditioner only.
   */
  class preonly_type
  {
  public:
368
    explicit preonly_type(std::string const& /*name*/) {}
369
370
    template <class LinOp, class X, class B, class P, class I>
    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
  // ===========================================================================

379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
  /// Adds default creators for linear solvers based on `mtl::compressed2D`.
  /**
   * Adds creators for full-matrix aware solvers.
   * - *cg*: conjugate gradient method, \see cg_solver_type
   * - *cgs*: stabilized conjugate gradient mtheod, \see cgs_solver_type
   * - *bicg*: BiCG mtheod, \see bicg_solver_type
   * - *bcgs*: stabilized bi-conjugate gradient method, \see bicgstab_type
   * - *bcgs2*: stabilized bi-conjugate gradient method, \see bicgstab2_type
   * - *qmr*: Quasi-Minimal Residual method, \see qmr_solver_type
   * - *tfqmr*: Transposed-Free Quasi-Minimal Residual method, \see tfqmr_solver_type
   * - *bcgsl*: stabilized BiCG(ell) method, \see bicgstab_ell_type
   * - *gmres*: Generalized minimal residula method, \see gmres_type
   * - *fgmres*: Flexible GMRES, \see fgmres_type
   * - *minres*: Minimal residul method, \see minres_solver_type
   * - *idrs*: Induced Dimension Reduction method, \see idr_s_type
   * - *gcr*: Generalized conjugate residual method, \see gcr_type
   * - *preonly*: Just a pply a preconditioner, \see preonly_type
   * - *umfpack*: external UMFPACK solver, \see UmfpackRunner
   **/
Praetorius, Simon's avatar
Praetorius, Simon committed
398
399
  template <class Mat, class Vec>
  class DefaultCreators< LinearSolverInterface<Mat,Vec> >
400
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
401
    using SolverBase = LinearSolverInterface<Mat,Vec>;
402

403
    template <class ITLSolver>
404
    using SolverCreator
Praetorius, Simon's avatar
Praetorius, Simon committed
405
      = typename LinearSolver<Mat,Vec,KrylovRunner<Mat,Vec,ITLSolver>>::Creator;
406

407
#ifdef HAVE_UMFPACK
408
    using UmfpackSolverCreator
Praetorius, Simon's avatar
Praetorius, Simon committed
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
428
      Map::addCreator("bicgstab", bcgs);
      Map::addCreator("bcgs", bcgs);
429

430
      auto bcgs2 = new SolverCreator<bicgstab2_type>;
431
432
      Map::addCreator("bicgstab2", bcgs2);
      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
442
      Map::addCreator("bicgstab_ell", bcgsl);
      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
463
464
      // default iterative solver
      Map::addCreator("default", gmres);

Praetorius, Simon's avatar
Praetorius, Simon committed
465
      init_direct(std::is_same<typename Dune::FieldTraits<typename Mat::value_type>::real_type, double>{});
466
467
468
469
470
    }

    static void init_direct(std::false_type) {}
    static void init_direct(std::true_type)
    {
471
#ifdef HAVE_UMFPACK
472
      auto umfpack = new UmfpackSolverCreator;
473
      Map::addCreator("umfpack", umfpack);
Praetorius, Simon's avatar
Praetorius, Simon committed
474
475

      // default direct solvers
476
      Map::addCreator("direct", umfpack);
477
#endif
478
479
    }
  };
480

481
} // end namespace AMDiS