CGSolver.hh 3.19 KB
Newer Older
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<typename VectorType>
10
  CGSolver<VectorType>::CGSolver(std::string name) 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
    : OEMSolver<VectorType>(name) 
  {}

  template<typename VectorType>
  CGSolver<VectorType>::~CGSolver() {}

  template<typename VectorType>
  int CGSolver<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
					VectorType *x, VectorType *b)
  {
    FUNCNAME("CGSolver::solve()");

    const double TOL = 1.e-30;

25
//    std::cout << "b: " << std::endl;
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);
    }

38
//    std::cout << "x: " << std::endl;
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);
47
//    std::cout << "-Ax: " << std::endl;
48
49
50
//    print(r);
    *r += *b;

51
//    std::cout << "r: " << std::endl;
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
61
62
//     std::cout << "Berechne Alpha: " << std::endl;
//     std::cout << "r: " << std::endl;
63
//     print(r);
64
//     std::cout << "p: " << std::endl;
65
66
//     print(p);
    double alpha = *r * *p;
67
//    std::cout << "alpha=" << alpha << std::endl;
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++) {
78
	//    std::cout << "---------------------------Ein Schritt---------------------------" << std::endl;
79
80
      // v = Ap
      matVec->matVec(NoTranspose, *p, *v);
81
//       std::cout << "p: " << std::endl;
82
//       print(p);
83
//       std::cout << "Ap: " << std::endl;
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)) {
111
112
113
114
115
116
117
	*p = *x;
	*p *= -1.0;
	
	// r = b - Ax
	matVec->matVec(NoTranspose, *p, *r);
	*r += *b;
	
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);
  }
}