NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
Functions
gmres.h File Reference
+ Include dependency graph for gmres.h:

Go to the source code of this file.

Functions

template<class Matrix , class Vector >
void Update (Vector &x, int k, Matrix &h, Vector &s, Vector v[])
 
template<class Real >
Real abs (Real x)
 
template<class Operator , class Vector , class Preconditioner , class Matrix , class Real >
int GMRES (const Operator &A, Vector &x, const Vector &b, const Preconditioner &M, Matrix &H, int &m, int &max_iter, Real &tol)
 
template<class Real >
void GeneratePlaneRotation (Real &dx, Real &dy, Real &cs, Real &sn)
 
template<class Real >
void ApplyPlaneRotation (Real &dx, Real &dy, Real &cs, Real &sn)
 

Function Documentation

template<class Matrix , class Vector >
void Update ( Vector &  x,
int  k,
Matrix &  h,
Vector &  s,
Vector  v[] 
)

Definition at line 25 of file gmres.h.

26 {
27  Vector y(s);
28 
29  // Backsolve:
30  for (int i = k; i >= 0; i--) {
31  y(i) /= h(i,i);
32  for (int j = i - 1; j >= 0; j--)
33  y(j) -= h(j,i) * y(i);
34  }
35 
36  for (int j = 0; j <= k; j++)
37  x += v[j] * y(j);
38 }
template<class Real >
Real abs ( Real  x)

Definition at line 43 of file gmres.h.

44 {
45  return (x > 0 ? x : -x);
46 }
template<class Operator , class Vector , class Preconditioner , class Matrix , class Real >
int GMRES ( const Operator &  A,
Vector &  x,
const Vector &  b,
const Preconditioner &  M,
Matrix &  H,
int &  m,
int &  max_iter,
Real &  tol 
)

Definition at line 52 of file gmres.h.

55 {
56  Real resid;
57  int i, j = 1, k;
58  Vector s(m+1), cs(m+1), sn(m+1), w;
59 
60  Real normb = norm(M.solve(b));
61  Vector r = M.solve(b - A * x);
62  Real beta = norm(r);
63 
64  if (normb == 0.0)
65  normb = 1;
66 
67  if ((resid = norm(r) / normb) <= tol) {
68  tol = resid;
69  max_iter = 0;
70  return 0;
71  }
72 
73  Vector *v = new Vector[m+1];
74 
75  while (j <= max_iter) {
76  v[0] = r * (1.0 / beta); // ??? r / beta
77  s = 0.0;
78  s(0) = beta;
79 
80  for (i = 0; i < m && j <= max_iter; i++, j++) {
81  w = M.solve(A * v[i]);
82  for (k = 0; k <= i; k++) {
83  H(k, i) = dot(w, v[k]);
84  w -= H(k, i) * v[k];
85  }
86  H(i+1, i) = norm(w);
87  v[i+1] = w * (1.0 / H(i+1, i)); // ??? w / H(i+1, i)
88 
89  for (k = 0; k < i; k++)
90  ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
91 
92  GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
93  ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
94  ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
95 
96  if ((resid = abs(s(i+1)) / normb) < tol) {
97  Update(x, i, H, s, v);
98  tol = resid;
99  max_iter = j;
100  delete [] v;
101  return 0;
102  }
103  }
104  Update(x, m - 1, H, s, v);
105  r = M.solve(b - A * x);
106  beta = norm(r);
107  if ((resid = beta / normb) < tol) {
108  tol = resid;
109  max_iter = j;
110  delete [] v;
111  return 0;
112  }
113  }
114 
115  tol = resid;
116  delete [] v;
117  return 1;
118 }
void Update(Vector &x, int k, Matrix &h, Vector &s, Vector v[])
Definition: gmres.h:25
void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
Definition: gmres.h:143
Real abs(Real x)
Definition: gmres.h:43
void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
Definition: gmres.h:125
A::value_type dot(const A &a, const B &b)
Dot product of array or expression.
Definition: util.ixx:91
A::value_type norm(const A &a)
Euklidean norm of array or expression.
Definition: util.hxx:54

+ Here is the call graph for this function:

template<class Real >
void GeneratePlaneRotation ( Real &  dx,
Real &  dy,
Real &  cs,
Real &  sn 
)

Definition at line 125 of file gmres.h.

126 {
127  if (dy == 0.0) {
128  cs = 1.0;
129  sn = 0.0;
130  } else if (abs(dy) > abs(dx)) {
131  Real temp = dx / dy;
132  sn = 1.0 / sqrt( 1.0 + temp*temp );
133  cs = temp * sn;
134  } else {
135  Real temp = dy / dx;
136  cs = 1.0 / sqrt( 1.0 + temp*temp );
137  sn = temp * cs;
138  }
139 }
Real abs(Real x)
Definition: gmres.h:43

+ Here is the call graph for this function:

template<class Real >
void ApplyPlaneRotation ( Real &  dx,
Real &  dy,
Real &  cs,
Real &  sn 
)

Definition at line 143 of file gmres.h.

144 {
145  Real temp = cs * dx + sn * dy;
146  dy = -sn * dx + cs * dy;
147  dx = temp;
148 }