MTLPfcPrecon.inc.hpp 1.77 KB
Newer Older
1
namespace AMDiS
2
{
3

4
5
6
  template <class Matrix, class Vector>
  void MTLPfcPrecon<Matrix, Vector>::init(Matrix const& A)
  {
7
    test_exit(tau, "Preconditioner not initialized!");
8
    Super::init(A);
9

10
11
    using size_type = typename MTLMatrix::size_type;
    size_type n = num_rows(getM());
12

13
    double delta = std::sqrt(M0 * getTau());
14
15
    msg("delta = ", delta);

16
17
    // helper-matrix MpL = M + delta*L
    MpL.change_dim(n, n);
18
    MpL = getM() + (1.0/delta) * getL0();
19
20
21
22

    // helper-matrix MpL = M + sqrt(delta)*L
    MpL2.change_dim(n, n);
    MpL2 = getM() + std::sqrt(delta) * getL();
23

24
25
26
27
28
29
30
31
32
    // temporary variables
    y0.change_dim(num_rows(getM()));
    y1.change_dim(num_rows(getM()));
    tmp.change_dim(num_rows(getM()));
  }


  template <class Matrix, class Vector>
  void MTLPfcPrecon<Matrix, Vector>::solve(Vector const& b, Vector& x) const
33
  {
34
    x.change_dim(num_rows(b));
35

36
37
38
    Vector x0(x[Super::getColRange(0_c)]);
    Vector x1(x[Super::getColRange(1_c)]);
    Vector x2(x[Super::getColRange(2_c)]);
39

40
41
    const Vector b0(b[Super::getRowRange(0_c)]);
    const Vector b1(b[Super::getRowRange(1_c)]);
42
43
    const Vector b2(b[Super::getRowRange(2_c)]);

44
    double delta = std::sqrt(M0 * getTau());
45

46
47
48
49
50
51
52
53
54
55
56
57
    solveM(y0, b0);		// M*y0 = b0
    y1 = getL0() * y0;		// y1 = K*y0
    tmp = b1 - y1;		// tmp := b1 - tau*y1

    solveMpL(y1, tmp);		// (M + delta*K) * y1 = tmp
    x0 = y0 + (1.0/delta)*y1;	// x0 = y0 + (1/delta)*y1

    tmp = getM() * y1;		// tmp := M*y1
    solveMpL2(y1, tmp);		// (M+eps*sqrt(delta)K) * y1 = tmp
    tmp = getM() * y1;		// tmp := M*y1
    solveMpL2(x1, tmp);		// (M+eps*sqrt(delta)K) * x1 = tmp

58
    x0-= (1.0/delta)*x1;	// x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1)
59
60
61
62
63
64

    tmp = b2 - getL() * x1;	// tmp := b2 - K*x1
    solveM(x2, tmp);
  }

} // end namespace AMDiS