ITL_Solver.hpp 14.5 KB
Newer Older
1
2
3
4
/** \file ITL_Solver.h */

#pragma once

5
6
#include <string>

7
// MTL4 includes
8
9
10
11
12
13
14
15
16
17
18
#include <boost/numeric/itl/itl.hpp>
// #include <boost/numeric/itl/krylov/bicg.hpp>
// #include <boost/numeric/itl/krylov/bicgstab_2.hpp>
// #include <boost/numeric/itl/krylov/bicgstab_ell.hpp>
// #include <boost/numeric/itl/krylov/bicgstab.hpp>
// #include <boost/numeric/itl/krylov/cg.hpp>
// #include <boost/numeric/itl/krylov/cgs.hpp>
// #include <boost/numeric/itl/krylov/gmres.hpp>
// #include <boost/numeric/itl/krylov/idr_s.hpp>
// #include <boost/numeric/itl/krylov/qmr.hpp>
// #include <boost/numeric/itl/krylov/tfqmr.hpp>
19
20

// more solvers defined in AMDiS
21
22
23
24
25
26
27
#include <amdis/linear_algebra/mtl/itl/minres.hpp>
#include <amdis/linear_algebra/mtl/itl/gcr.hpp>
#include <amdis/linear_algebra/mtl/itl/fgmres.hpp>
#include <amdis/linear_algebra/mtl/itl/fgmres_householder.hpp>
#include <amdis/linear_algebra/mtl/itl/gmres2.hpp>
#include <amdis/linear_algebra/mtl/itl/gmres_householder.hpp>
#include <amdis/linear_algebra/mtl/itl/preonly.hpp>
28

29
30
31
32
33
#include <amdis/CreatorMap.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/linear_algebra/mtl/LinearSolver.hpp>
#include <amdis/linear_algebra/mtl/KrylovRunner.hpp>
#include <amdis/linear_algebra/mtl/UmfpackRunner.hpp>
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

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:
49
    explicit cg_solver_type(std::string /*name*/) {}
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
    {
      return itl::cg(A, x, b, l, r, iter);
    }
  };


  /**
   * \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:
69
    explicit cgs_solver_type(std::string /*name*/) {}
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter)
    {
      return itl::cgs(A, x, b, l, iter);
    }
  };


  /**
   * \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:
89
    explicit bicg_solver_type(std::string /*name*/) {}
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter)
    {
      return itl::bicg(A, x, b, l, iter);
    }
  };


  /**
   * \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:
109
    explicit bicgstab_type(std::string /*name*/) {}
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter)
    {
      return itl::bicgstab(A, x, b, l, iter);
    }
  };


  /**
   * \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:
129
    explicit bicgstab2_type(std::string /*name*/) {}
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter)
    {
      return itl::bicgstab_2(A, x, b, l, iter);
    }
  };


  /**
   * \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:
148
    explicit qmr_solver_type(std::string /*name*/) {}
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
    {
      return itl::qmr(A, x, b, l, r, iter);
    }
  };


  /**
   * \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
   * (TFQMR). Does not use preconditioning currently.
   */
  class tfqmr_solver_type
  {
  public:
168
    explicit tfqmr_solver_type(std::string /*name*/) {}
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
    {
      return itl::tfqmr(A, x, b, l, r, iter);
    }
  };


  /**
   * \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:
189
    explicit bicgstab_ell_type(std::string name) : ell(3)
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
    {
      Parameters::get(name + "->ell", ell);
    }
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
    {
      return itl::bicgstab_ell(A, x, b, l, r, iter, ell);
    }
  };


  /**
   * \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:
214
    explicit gmres_type(std::string name) : restart(30), ortho(GRAM_SCHMIDT)
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
    {
      Parameters::get(name + "->restart", restart);
      Parameters::get(name + "->orthogonalization", ortho);
    }
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
    {
      switch (ortho)
      {
      default:
      case GRAM_SCHMIDT:
        return itl::gmres2(A, x, b, l, r, iter, restart);
        break;
      case HOUSEHOLDER:
        return itl::gmres_householder(A, x, b, l, iter, restart);
        break;
      }
    }
233

234
235
236
  private:
    static constexpr int GRAM_SCHMIDT = 1;
    static constexpr int HOUSEHOLDER = 2;
237

238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
    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)
   */
  class idr_s_type
  {
    int s;
  public:
259
    explicit idr_s_type(std::string name) : s(30)
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
    {
      Parameters::get(name + "->s", s);
    }
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
    {
      return itl::idr_s(A, x, b, l, r, iter, s);
    }
  };


  /**
   * \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:
282
    explicit minres_solver_type(std::string /*name*/) {}
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
    {
      return itl::minres(A, x, b, l, r, iter);
    }
  };


  /**
   * \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:
305
    explicit gcr_type(std::string name) : restart(30)
306
307
308
309
310
311
312
313
314
    {
      Parameters::get(name + "->restart", restart);
    }
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
    {
      return itl::gcr(A, x, b, l, r, iter, restart);
    }
  };
315

316
317
318
319
320
321
322
323
324
325
326
327
328
329

  /**
   * \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:
330
    explicit fgmres_type(std::string name) : restart(30), orthogonalization(GRAM_SCHMIDT)
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
    {
      Parameters::get(name + "->restart", restart);
      Parameters::get(name + "->orthogonalization", orthogonalization);
    }
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
    {
      switch (orthogonalization)
      {
      default:
      case GRAM_SCHMIDT:
        return itl::fgmres(A, x, b, l, r, iter, restart);
        break;
      case HOUSEHOLDER:
        return itl::fgmres_householder(A, x, b, r, iter, restart);
        break;
      }
    }
349

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

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

358

359
360
361
362
363
364
365
366
367
368
  /**
   * \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:
369
    explicit preonly_type(std::string /*name*/) {}
370
371
372
373
374
375
376
    template <class LinOp, class X, class B, class L, class R, class I>
    int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& /*r*/, I& iter)
    {
      return itl::preonly(A, x, b, l, iter);
    }
  };

377

378
379
  // ===========================================================================

380

381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
  /// 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
   **/
400
  template <class T, class Param, class Vector>
401
  class DefaultCreators<LinearSolverInterface<mtl::compressed2D<T, Param>, Vector, Vector>>
402
403
404
  {
    using Matrix = mtl::compressed2D<T, Param>;
    using SolverBase = LinearSolverInterface<Matrix, Vector, Vector>;
405

406
    template <class ITLSolver>
407
408
    using SolverCreator
      = typename LinearSolver<Matrix, Vector,
409
          KrylovRunner<ITLSolver, Matrix, BaseVector_t<Vector>>>::Creator;
410

411
#ifdef HAVE_UMFPACK
412
    using UmfpackSolverCreator
413
      = typename LinearSolver<Matrix, Vector,
414
          UmfpackRunner<Matrix, BaseVector_t<Vector>>>::Creator;
415
#endif
416

417
    using Map = CreatorMap<SolverBase>;
418

419
  public:
420
    static void init()
421
    {
422
      auto cg = new SolverCreator<cg_solver_type>;
423
      Map::addCreator("cg", cg);
424

425
      auto cgs = new SolverCreator<cgs_solver_type>;
426
      Map::addCreator("cgs", cgs);
427

428
      auto bicg = new SolverCreator<bicg_solver_type>;
429
      Map::addCreator("bicg", bicg);
430

431
      auto bcgs = new SolverCreator<bicgstab_type>;
432
433
      Map::addCreator("bicgstab", bcgs);
      Map::addCreator("bcgs", bcgs);
434

435
      auto bcgs2 = new SolverCreator<bicgstab2_type>;
436
437
      Map::addCreator("bicgstab2", bcgs2);
      Map::addCreator("bcgs2", bcgs2);
438

439
      auto qmr = new SolverCreator<qmr_solver_type>;
440
      Map::addCreator("qmr", qmr);
441

442
      auto tfqmr = new SolverCreator<tfqmr_solver_type>;
443
      Map::addCreator("tfqmr", tfqmr);
444

445
      auto bcgsl = new SolverCreator<bicgstab_ell_type>;
446
447
      Map::addCreator("bicgstab_ell", bcgsl);
      Map::addCreator("bcgsl", bcgsl);
448

449
      auto gmres = new SolverCreator<gmres_type>;
450
      Map::addCreator("gmres", gmres);
451

452
      auto fgmres = new SolverCreator<fgmres_type>;
453
      Map::addCreator("fgmres", fgmres);
454

455
      auto minres = new SolverCreator<minres_solver_type>;
456
      Map::addCreator("minres", minres);
457

458
      auto idrs = new SolverCreator<idr_s_type>;
459
      Map::addCreator("idrs", idrs);
460

461
      auto gcr = new SolverCreator<gcr_type>;
462
      Map::addCreator("gcr", gcr);
463

464
      auto preonly = new SolverCreator<preonly_type>;
465
      Map::addCreator("preonly", preonly);
466

467
#ifdef HAVE_UMFPACK
468
      auto umfpack = new UmfpackSolverCreator;
469
470
      Map::addCreator("umfpack", umfpack);
      Map::addCreator("direct", umfpack);
471
#endif
472

473
      Map::addCreator("default", gmres);
474
475
    }
  };
476

477
} // end namespace AMDiS