GMResSolver.hh 4.77 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
#include "GMResSolver.h"
2
3
4
5
6
7
8
#include "Preconditioner.h"

namespace AMDiS {

  template<typename VectorType>
  GMResSolver<VectorType>::GMResSolver(::std::string name)
    : OEMSolver<VectorType>(name),
Thomas Witkowski's avatar
Thomas Witkowski committed
9
10
11
      TOL_(1.e-25),
      restart_(0),
      w_(NULL)
12
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
13
14
15
    FUNCNAME("GMResSolver::GMResSolver()");
    GET_PARAMETER(0, name + "->restart", "%d", &restart_);
    TEST_EXIT(restart_)("restart not set\n");
16
17
18
19
20
21
22
23
24
  }

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

  template<typename VectorType>
  void GMResSolver<VectorType>::init()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
    // Create a new vector w.
    w_ = this->vectorCreator->create();
    z_ = this->vectorCreator->create();
    r_ = this->vectorCreator->create();

    // Get memory for the vector pointers and create the pointers afterwards.
    v = GET_MEMORY(VectorType*, restart_ + 1);
    for (int i = 0; i <= restart_; i++) {
      v[i] = this->vectorCreator->create();
    }

    // Resize all fields to the corresponding size.
    gamma.resize(restart_ + 1);
    c.resize(restart_);
    s.resize(restart_);
40

Thomas Witkowski's avatar
Thomas Witkowski committed
41
    y_.resize(restart_);
42

Thomas Witkowski's avatar
Thomas Witkowski committed
43
44
45
    h.resize(restart_ + 1);
    for (int i = 0; i <= restart_; i++) {
      h[i].resize(restart_);
46
47
48
49
50
51
    }
  }

  template<typename VectorType>
  void GMResSolver<VectorType>::exit()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
52
53
54
55
56
57
58
59
    // Empty used memory.
    if (w_) {      
      this->vectorCreator->free(w_);
      this->vectorCreator->free(z_);
      this->vectorCreator->free(r_);

      for (int i = 0; i <= restart_; i++) {
	this->vectorCreator->free(v[i]);
60
61
62
63
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
64

65
  template<typename VectorType>
Thomas Witkowski's avatar
Thomas Witkowski committed
66
67
68
  int GMResSolver<VectorType>::gmres(MatVecMultiplier<VectorType> *matVec,
				      VectorType *x, 
				      VectorType *b)
69
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
70
    FUNCNAME("GMResSolver::gmres()");
71

Thomas Witkowski's avatar
Thomas Witkowski committed
72
73
74
75
    // r_0 = b - Ax, where r_0 is already stored as the first vector in the
    // matrix V.
    matVec->matVec(NoTranspose, *x, *v[0]);
    xpay(-1.0, *b, *v[0]);
76

Thomas Witkowski's avatar
Thomas Witkowski committed
77
78
79
    // Left preconditioning, if required.
    if (this->leftPrecon) {
      this->leftPrecon->precon(v[0]);
80
81
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
82
83
    // Compute the norm of r_0 = v_0.
    gamma[0] = norm(v[0]);
84

Thomas Witkowski's avatar
Thomas Witkowski committed
85
86
    // Normalize v_0.
    *v[0] *= (1 / gamma[0]);
87

Thomas Witkowski's avatar
Thomas Witkowski committed
88
89
90
91
92
93
94
95
96
97
98
    // If the norm of gamma_0 is less than the tolerance bounday, x is already
    // a good solution for Ax = b.
    if (gamma[0] < this->tolerance) {
      this->residual = gamma[0];
      return(0);
    }    
    
    // Main loop of the GMRES algorithm.
    for (int j = 0; j < restart_; j++) {
      // w_j = A * v_j;
      matVec->matVec(NoTranspose, *v[j], *w_);
99
      
Thomas Witkowski's avatar
Thomas Witkowski committed
100
101
102
103
      // Preconditioning of w_j
      if (this->leftPrecon) {
	this->leftPrecon->precon(w_);
      }
104

Thomas Witkowski's avatar
Thomas Witkowski committed
105
106
107
108
      // w_j = A * v_j - \sum(h_ij * v_i)
      for (int i = 0; i <= j; i++) {
	// h_ij = (v_i, A * v_j)_2
	h[i][j] = *w_ * (*v[i]);
109

Thomas Witkowski's avatar
Thomas Witkowski committed
110
111
	// Update w_j
	axpy(-h[i][j], *v[i], *w_);
112
113
      }
      
Thomas Witkowski's avatar
Thomas Witkowski committed
114
115
116
117
118
119
120
121
122
      // h_{j + 1}{j} = \| w_j \|_2
      h[j + 1][j] = norm(w_);
     
      // Compute h_ij and h_{i + 1}{j} using parameters of the Givens rotation.
      for (int i = 0; i < j; i++) {
	double t1 = c[i] * h[i][j] + s[i] * h[i + 1][j];
	double t2 = -s[i] * h[i][j] + c[i] * h[i + 1][j];
	h[i][j] = t1;
	h[i + 1][j] = t2;
123
124
      }
      
Thomas Witkowski's avatar
Thomas Witkowski committed
125
126
127
128
129
130
131
132
133
134
135
      // Compute new beta
      double beta = sqrt(h[j][j] * h[j][j] + h[j + 1][j] * h[j + 1][j]);
      // Update s_j and c_j
      s[j] = h[j + 1][j] / beta;
      c[j] = h[j][j] / beta;
      // h_jj = beta
      h[j][j] = beta;

      // Update gammas
      gamma[j + 1] = -s[j] * gamma[j];
      gamma[j] *= c[j];
136
      
Thomas Witkowski's avatar
Thomas Witkowski committed
137
138
139
      // Set v_{j + 1} to w_j and normalize it.
      *v[j + 1] = *w_;
      *v[j + 1] *= 1 / h[j + 1][j];      
140
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
141
142
143
144
145
    
    // Compute the solution.
    for (int i = restart_ - 1; i >= 0; i--) {
      for (int k = i + 1; k < restart_; k++) {
 	gamma[i] -= h[i][k] * gamma[k];
146
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
147
148
      gamma[i] /= h[i][i];
      axpy(gamma[i], *v[i], *x);
149
150
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
153
154
    // Update the new residual.
    this->residual = abs(gamma[restart_]);
       
    return restart_;
155
156
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
157

158
159
  template<typename VectorType>
  int GMResSolver<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
Thomas Witkowski's avatar
Thomas Witkowski committed
160
161
 					    VectorType *x, 
 					    VectorType *b)
162
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
163
    FUNCNAME("GMResSolver::solveSystem()");    
164
    
Thomas Witkowski's avatar
Thomas Witkowski committed
165
    double old_res = -1.0;
166
    
Thomas Witkowski's avatar
Thomas Witkowski committed
167
168
169
    // 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_) {
170
171
172
173
174
      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's avatar
Thomas Witkowski committed
175
    
176
177
    START_INFO();
    for (int iter = 0; iter <= this->max_iter; iter++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
180
      // Solve Ax=b using GMRES(m).
      int k = gmres(matVec, x, b);
      // Check the new residual.
181
182
183
184
185
186
187
188
189
190
      if (SOLVE_INFO(iter, this->residual, &old_res))
	return(iter);
      TEST_EXIT(k != 0)("this must not happen\n");
    }
  
    ERROR_EXIT("Should not be reached\n");
    return(0);
  }

}