CGSolver.hh 3.21 KB
 Peter Gottschling committed Feb 15, 2008 1 2 3 4 5 6 7 8 9 #include "DOFVector.h" #include "DOFMatrix.h" #include "MatVecMultiplier.h" #include "StlVector.h" #include "V3Vector.h" namespace AMDiS { template  Thomas Witkowski committed Aug 29, 2008 10  CGSolver::CGSolver(std::string name)  Peter Gottschling committed Feb 15, 2008 11 12 13 14 15 16 17 18  : OEMSolver(name) {} template CGSolver::~CGSolver() {} template int CGSolver::solveSystem(MatVecMultiplier *matVec,  19  VectorType *x, VectorType *b, bool reuseMatrix)  Peter Gottschling committed Feb 15, 2008 20 21 22 23 24  { FUNCNAME("CGSolver::solve()"); const double TOL = 1.e-30;  Thomas Witkowski committed Aug 29, 2008 25 // std::cout << "b: " << std::endl;  Peter Gottschling committed Feb 15, 2008 26 27 28 29 30 31 32 33 34 35 36 37 // 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); }  Thomas Witkowski committed Aug 29, 2008 38 // std::cout << "x: " << std::endl;  Peter Gottschling committed Feb 15, 2008 39 40 41 42 43 44 45 46 // print(x); // p is temporaly used to calculate -x. *p = *x; *p *= -1.0; // r = b - Ax matVec->matVec(NoTranspose, *p, *r);  Thomas Witkowski committed Aug 29, 2008 47 // std::cout << "-Ax: " << std::endl;  Peter Gottschling committed Feb 15, 2008 48 49 50 // print(r); *r += *b;  Thomas Witkowski committed Aug 29, 2008 51 // std::cout << "r: " << std::endl;  Peter Gottschling committed Feb 15, 2008 52 53 54 55 56 57 58 59 60 // 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  Thomas Witkowski committed Aug 29, 2008 61 62 // std::cout << "Berechne Alpha: " << std::endl; // std::cout << "r: " << std::endl;  Peter Gottschling committed Feb 15, 2008 63 // print(r);  Thomas Witkowski committed Aug 29, 2008 64 // std::cout << "p: " << std::endl;  Peter Gottschling committed Feb 15, 2008 65 66 // print(p); double alpha = *r * *p;  Thomas Witkowski committed Aug 29, 2008 67 // std::cout << "alpha=" << alpha << std::endl;  Peter Gottschling committed Feb 15, 2008 68 69 70 71 72 73 74 75 76 77  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++) {  Thomas Witkowski committed Aug 29, 2008 78  // std::cout << "---------------------------Ein Schritt---------------------------" << std::endl;  Peter Gottschling committed Feb 15, 2008 79 80  // v = Ap matVec->matVec(NoTranspose, *p, *v);  Thomas Witkowski committed Aug 29, 2008 81 // std::cout << "p: " << std::endl;  Peter Gottschling committed Feb 15, 2008 82 // print(p);  Thomas Witkowski committed Aug 29, 2008 83 // std::cout << "Ap: " << std::endl;  Peter Gottschling committed Feb 15, 2008 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 // 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)) {  Thomas Witkowski committed May 23, 2008 111 112 113 114 115 116 117  *p = *x; *p *= -1.0; // r = b - Ax matVec->matVec(NoTranspose, *p, *r); *r += *b;  Peter Gottschling committed Feb 15, 2008 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132  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); } }