35 template <
class Matrix,
class Vector >
38 Update(Vector &x,
int k, Matrix &h, Vector &s, std::vector<Vector> &v)
43 for (
int i = k; i >= 0; i--) {
45 for (
int j = i - 1; j >= 0; j--)
46 y(j) -= h(j,i) * y(i);
49 for (
int j = 0; j <= k; j++) {
54 template <
class Real >
58 return (x > 0 ? x : -x);
61 double norm(
const CPPL::dcovector& V) {
64 for (i = 0; i < V.l; i++) {
70 double dot(
const CPPL::dcovector& V,
const CPPL::dcovector& W) {
73 for (i = 0; i < V.l; i++) {
81 void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
86 }
else if (abs(dy) > abs(dx)) {
88 sn = 1.0 / sqrt( 1.0 + temp*temp );
92 cs = 1.0 / sqrt( 1.0 + temp*temp );
99 void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
101 Real temp = cs * dx + sn * dy;
102 dy = -sn * dx + cs * dy;
107 template <
class Operator,
class Vector,
class Preconditioner,
108 class Matrix,
class Real >
110 GMRES(
const Operator &A, Vector &x,
const Vector &b,
111 const Preconditioner &M, Matrix &H,
int &
m,
int &max_iter,
116 Vector s(m+1), cs(m+1), sn(m+1);
120 Real normb = norm(M.solve(b));
123 Vector r = M.solve(b - A * x);
131 if (
false && (resid = norm(r) / normb) <= tol) {
138 std::vector<Vector> v(m+1, Vector());
141 while (j <= max_iter) {
142 v[0] = r * (1.0 / beta);
144 for (i = 0; i < s.l; i++) s(i) = 0;
147 for (i = 0; i < m && j <= max_iter; i++, j++) {
148 w = M.solve(A * v[i]);
150 for (k = 0; k <= i; k++) {
151 H(k, i) = dot(w, v[k]);
157 v[i+1] = w * (1.0 / H(i+1, i));
159 for (k = 0; k < i; k++)
160 ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
162 GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
163 ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
164 ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
167 if ((resid = abs(s(i+1)) / normb) < tol) {
168 Update(x, i, H, s, v);
175 Update(x, i - 1, H, s, v);
176 r = M.solve(b - A * x);
179 if ((resid = beta / normb) < tol) {
static const double m
mass of electron in grams