MTLPfcPrecon.hpp 3.49 KB
Newer Older
1
2
3
#pragma once

#include "BlockPreconditioner.hpp"
4
#include <amdis/common/Literals.hpp>
5

6
7
namespace AMDiS
{
8
  template <class Matrix, class Vector>
9
10
  class MTLPfcPrecon
      : public BlockPreconditioner<Matrix, Vector>
11
12
13
14
15
  {
  public:
    using Self = MTLPfcPrecon;
    using Super = BlockPreconditioner<Matrix, Vector>;
    using PreconBase = typename Super::PreconBase;
16

17
18
19
    using MTLMatrix = BaseMatrix_t<Matrix>;
    using MTLVector = Vector;
    using LinearSolverType = LinearSolverInterface<MTLMatrix, MTLVector>;
20

21
22
23
24
25
26
  public:
    struct Creator : CreatorInterface<PreconBase>
    {
      virtual std::shared_ptr<PreconBase> create() override
      {
        precon = std::make_shared<Self>();
27
        return precon;
28
      }
29

30
      std::shared_ptr<Self> precon = NULL;
31
32
    };

33
34
35
36
  public:
    MTLPfcPrecon()
      : Super()
    {
37
38
      std::string solverNameM = "cg",
                  solverNameMpL = "cg",
39
                  solverNameMpL2 = "cg";
40
41
42
      Parameters::get("precon_pfc->M->solver", solverNameM);
      Parameters::get("precon_pfc->MpL->solver", solverNameMpL);
      Parameters::get("precon_pfc->MpL2->solver", solverNameMpL2);
43

44
      CreatorInterfaceName<LinearSolverType>* creator;
45

46
      creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameM, "precon_pfc->M->solver") );
47
      solverM = creator->create("precon_pfc->M");
48

49
      creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameMpL, "precon_pfc->MpL->solver") );
50
      solverMpL = creator->create("precon_pfc->MpL");
51

52
      creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameMpL2, "precon_pfc->MpL2->solver") );
53
      solverMpL2 = creator->create("precon_pfc->MpL2");
54
    }
55

56
57
58
59
    virtual ~MTLPfcPrecon()
    {
      exit();
    }
60

61
62
63
64
65
    void setData(double* tau_, double M0_ = 1.0)
    {
      tau = tau_;
      M0 = M0_;
    }
66

67
68
    // implementation in MTLPfcPrecon.inc.hpp
    virtual void init(Matrix const& A) override;
69

70
    virtual void exit() override {}
71

72
    // implementation in MTLPfcPrecon.inc.hpp
73
74
75
    virtual void solve(MTLVector const& b, MTLVector& x) const override;


76
77
78
    template <typename VectorX, typename VectorB>
    void solveM(VectorX& x, const VectorB& b) const
    {
79
      SolverInfo solverInfo("precon_pfc->M");
80
81
      solverM->solve(getM(), x, b, solverInfo);
    }
82

83
84
85
    template <typename VectorX, typename VectorB>
    void solveMpL(VectorX& x, const VectorB& b) const
    {
86
      SolverInfo solverInfo("precon_pfc->MpL");
87
      solverMpL->solve(MpL, x, b, solverInfo);
88
    }
89

90
91
92
    template <typename VectorX, typename VectorB>
    void solveMpL2(VectorX& x, const VectorB& b) const
    {
93
      SolverInfo solverInfo("precon_pfc->MpL2");
94
      solverMpL2->solve(MpL2, x, b, solverInfo);
95
    }
96

Praetorius, Simon's avatar
Praetorius, Simon committed
97
98
99
    MTLMatrix const& getM()  const { return matM  ? *matM  : (*Super::A)(_2, _2); } // < u, v >
    MTLMatrix const& getL0() const { return matL0 ? *matL0 : (*Super::A)(_1, _0); } // < M0*tau*grad(u), grad(v) >
    MTLMatrix const& getL()  const { return matL  ? *matL  : (*Super::A)(_2, _1); } // < grad(u), grad(v) >
100

101
    double getTau() const { return *tau; }
102
103

  protected:
104
105
    double* tau = NULL;
    double M0 = 1.0;
106

107
108
109
    MTLMatrix* matM = NULL;
    MTLMatrix* matL = NULL;
    MTLMatrix* matL0 = NULL;
110

111
112
    MTLMatrix MpL;
    MTLMatrix MpL2;
113

114
115
116
117
118
119
120
121
    mutable MTLVector y0;
    mutable MTLVector y1;
    mutable MTLVector tmp;

    std::shared_ptr<LinearSolverType> solverM;
    std::shared_ptr<LinearSolverType> solverMpL;
    std::shared_ptr<LinearSolverType> solverMpL2;
  };
122

123
124
125
} // end namespace AMDiS

#include "MTLPfcPrecon.inc.hpp"