#include "DOFVector.h" #include "DOFMatrix.h" #include "MatVecMultiplier.h" #include "StlVector.h" #include "V3Vector.h" namespace AMDiS { template CGSolver::CGSolver(std::string name) : OEMSolver(name) {} template CGSolver::~CGSolver() {} template int CGSolver::solveSystem(MatVecMultiplier *matVec, VectorType *x, VectorType *b, bool reuseMatrix) { FUNCNAME("CGSolver::solve()"); const double TOL = 1.e-30; // std::cout << "b: " << std::endl; // print(b); // If norm of b is smaller than the tolarance, we can assume b to be zero. // Hence, x = 0 is the solution of the linear system. if (norm(b) < TOL) { INFO(this->info,2)("b == 0, x = 0 is the solution of the linear system\n"); setValue(*x, 0.0); this->residual = 0.0; return(0); } // std::cout << "x: " << std::endl; // print(x); // p is temporaly used to calculate -x. *p = *x; *p *= -1.0; // r = b - Ax matVec->matVec(NoTranspose, *p, *r); // std::cout << "-Ax: " << std::endl; // print(r); *r += *b; // std::cout << "r: " << std::endl; // print(r); *p = *r; if (this->leftPrecon) { this->leftPrecon->precon(p); } // \alpha = (r,p)_2 // without precondition, \alpha = (r,p)_2 = (r,r)_2 = \|r\|_2^2 // std::cout << "Berechne Alpha: " << std::endl; // std::cout << "r: " << std::endl; // print(r); // std::cout << "p: " << std::endl; // print(p); double alpha = *r * *p; // std::cout << "alpha=" << alpha << std::endl; double alpha1, lambda, old_res = -1.0; START_INFO(); // Test if x is already a solution of the lineas system. if (SOLVE_INFO(0, sqrt(alpha), &old_res)) { return(0); } for (int iter = 1; iter <= this->max_iter; iter++) { // std::cout << "---------------------------Ein Schritt---------------------------" << std::endl; // v = Ap matVec->matVec(NoTranspose, *p, *v); // std::cout << "p: " << std::endl; // print(p); // std::cout << "Ap: " << std::endl; // print(v); // \lambda = alpha / (v,p)_2 lambda = *v * *p; if (abs(lambda) < TOL) { BREAK_INFO("(v,p)_2 = 0", iter, sqrt(alpha), &old_res); return(iter); } lambda = alpha / lambda; // x = x + \lambda * p axpy(lambda, *p, *x); // r = r - \lambda * v axpy(-lambda, *v, *r); *z = *r; if (this->leftPrecon) { this->leftPrecon->precon(z); } // \alpha_{m + 1} = (r,z)_2 // without precondition, \alpha_{m + 1} = (r,z)_2 = (r,r)_2 = \|r\|_2^2 alpha1 = *r * *z; // Check if x is a solution of the linear system. if (SOLVE_INFO(iter, sqrt(alpha1), &old_res)) { *p = *x; *p *= -1.0; // r = b - Ax matVec->matVec(NoTranspose, *p, *r); *r += *b; return(iter); } // p = z + (alpha_{m + 1} / alpha_{m}) * p xpay(alpha1 / alpha, *z, *p); // Put the value of the alpha_{m + 1} variable to alpha_{m} variable // for the next iteration. alpha = alpha1; } ERROR_EXIT("Should not be reached\n"); return(0); } }