MTLPreconPfc.hh 2.63 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/

void MTLPreconPfc::init(const BlockMatrix& A_, const MTLMatrix& fullMatrix_)
19
{
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
  assert(tau != NULL);
  super::init(A_, fullMatrix_);
  
  // helper-matrix MpL = M + sqrt(tau)*L
  MpL.change_dim(num_rows(getM()), num_cols(getM()));
  MpL = getM() + sqrt(*tau)*getL();
  
  // schur-komplement-Matrix S = M - 2sqrt(tau)*L + sqrt(tau)*L*M^(-1)*L*M^(-1)*L
  S = PfcSchurMatrix<MTLPreconPfc>(this);
    
  // temporary variables
  y0.change_dim(num_rows(getM()));
  y1.change_dim(num_rows(getM()));
  tmp.change_dim(num_rows(getM()));  
  
35
  // init sub-solvers
36
  Pid = new itl::pc::identity<MTLTypes::MTLMatrix, double>(getM());
37
  Pdiag = new itl::pc::diagonal<MTLTypes::MTLMatrix, double>(MpL);
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
  
  runnerM->init(A_, getM());
  runnerM2->init(A_, getM());
  runnerMpL->init(A_, MpL);
}


void MTLPreconPfc::solve(const MTLVector& b, MTLVector& x) const
{ FUNCNAME("MTLPreconPfc::solve()");

  x.change_dim(num_rows(b));
  
  const MTLVector b0(b[rows[0]]);
  const MTLVector b1(b[rows[1]]);
  const MTLVector b2(b[rows[2]]);
  
  MTLVector x0(x[rows[0]]);
  MTLVector x1(x[rows[1]]);
  MTLVector x2(x[rows[2]]);
  
  double delta = sqrt(getTau());
  
  solveM(y0, b0);
  y1 = getL() * y0;		// y1 = K*y0
  tmp = b1 - getTau()*y1;	// tmp := b1 - tau*y1

  solveMpK(y1, tmp);
  x0 = y0 + (1.0/delta)*y1;// x0 = y0 + (1/delta)*y1

  tmp = getM() * y1;		// tmp := M*y1
  solveS(S, x1, tmp);

  x0-= (1.0/delta)*x1;	// x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1)

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


template<>
template <typename VectorIn, typename VectorOut, typename Assign>
void PfcSchurMatrix<MTLPreconPfc>::mult(const VectorIn& b, VectorOut& x, Assign) const
{ FUNCNAME("PfcSchurMatrix<MTLPreconPfc>::mult()");

  int m = num_rows(super->getM());
  assert(int(size(b)) == m);
  assert(size(b) == size(x));
  
  double delta = sqrt(super->getTau());
  
  y1 = super->getL() * b;
  super->solveM(y2, y1);
  y3 = super->getL() * y2;

  y2 = super->getM() * b;
  y2-= (2.0*delta) * y1;
  y2+= delta * y3;

  Assign::apply(x, y2);
}