namespace AMDiS { template void MTLPfcPrecon::init(Matrix const& A) { test_exit(tau, "Preconditioner not initialized!"); Super::init(A); using size_type = typename MTLMatrix::size_type; size_type n = num_rows(getM()); double delta = std::sqrt(M0 * getTau()); msg("delta = ", delta); // helper-matrix MpL = M + delta*L MpL.change_dim(n, n); MpL = getM() + (1.0/delta) * getL0(); // helper-matrix MpL = M + sqrt(delta)*L MpL2.change_dim(n, n); MpL2 = getM() + std::sqrt(delta) * getL(); // temporary variables y0.change_dim(num_rows(getM())); y1.change_dim(num_rows(getM())); tmp.change_dim(num_rows(getM())); } template void MTLPfcPrecon::solve(Vector const& b, Vector& x) const { x.change_dim(num_rows(b)); Vector x0(x[Super::getColRange(0_c)]); Vector x1(x[Super::getColRange(1_c)]); Vector x2(x[Super::getColRange(2_c)]); const Vector b0(b[Super::getRowRange(0_c)]); const Vector b1(b[Super::getRowRange(1_c)]); const Vector b2(b[Super::getRowRange(2_c)]); double delta = std::sqrt(M0 * getTau()); 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 x0-= (1.0/delta)*x1; // x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1) tmp = b2 - getL() * x1; // tmp := b2 - K*x1 solveM(x2, tmp); } } // end namespace AMDiS