From: wolf Date: Mon, 17 May 1999 13:30:01 +0000 (+0000) Subject: Commit this to have it in the CVS, but I will delete it immediately after again since... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2eeea1002b99a5fe7d1456a7396afe4ee0de37ab;p=dealii-svn.git Commit this to have it in the CVS, but I will delete it immediately after again since it did not really work from the beginning. git-svn-id: https://svn.dealii.org/trunk@1331 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h index c5cd615dba..4d47b06bb3 100644 --- a/deal.II/lac/include/lac/solver_gmres.h +++ b/deal.II/lac/include/lac/solver_gmres.h @@ -136,6 +136,136 @@ class SolverGMRES : public Solver +template +void +Update(::Vector &x, int k, FullMatrix &h, + ::Vector &s, vector &v) +{ + ::Vector y(s); + + // Backsolve: + for (int i = k; i >= 0; i--) { + y(i) /= h(i,i); + for (int j = i - 1; j >= 0; j--) + y(j) -= h(j,i) * y(i); + } + + for (int j = 0; j <= k; j++) + x.add (y(j), v[j]); +} + + + + +template < class Matrix, class Vector, class Preconditioner> +int +solve2(const Matrix &A, Vector &x, const Vector &b, + const Preconditioner &M, int &m, int &max_iter, + double &tol) +{ + FullMatrix H(m+1,m); + + double resid; + int i, j = 1, k; + ::Vector s(m+1), cs(m+1), sn(m+1); + Vector w(x.size()); + + Vector tmp (b.size()); + M (tmp, b); + double normb = tmp.l2_norm(); + Vector r(b.size()); + A.residual (tmp, x, b); + M(r, tmp); + double beta = r.l2_norm(); + + if (normb == 0.0) + normb = 1; + + if ((resid = r.l2_norm() / normb) <= tol) { + tol = resid; + max_iter = 0; + return 0; + } + + vector v(m+1); + + while (j <= max_iter) { + v[0] = r; + v[0].scale (1.0 / beta); // ??? r / beta + s = 0.0; + s(0) = beta; + + for (i = 0; i < m && j <= max_iter; i++, j++) { + A.vmult (tmp, v[i]); + M(w, tmp); + for (k = 0; k <= i; k++) { + H(k, i) = w * v[k]; + w.sadd (-H(k, i), v[k]); + } + H(i+1, i) = w.l2_norm(); + v[i+1] = w; + v[i+1].scale (1.0 / H(i+1, i)); // ??? w / H(i+1, i) + + for (k = 0; k < i; k++) + ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k)); + + GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i)); + ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i)); + ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i)); + + if ((resid = fabs(s(i+1)) / normb) < tol) { + Update(x, i, H, s, v); + tol = resid; + max_iter = j; + return 0; + } + } + Update(x, i - 1, H, s, v); + A.residual (tmp, x, b); + M(r, tmp); + beta = r.l2_norm(); + if ((resid = beta / normb) < tol) { + tol = resid; + max_iter = j; + return 0; + } + } + + tol = resid; + return 1; +} + + +template +void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn) +{ + if (dy == 0.0) { + cs = 1.0; + sn = 0.0; + } else if (fabs(dy) > fabs(dx)) { + Real temp = dx / dy; + sn = 1.0 / sqrt( 1.0 + temp*temp ); + cs = temp * sn; + } else { + Real temp = dy / dx; + cs = 1.0 / sqrt( 1.0 + temp*temp ); + sn = temp * cs; + } +} + + +template +void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn) +{ + Real temp = cs * dx + sn * dy; + dy = -sn * dx + cs * dy; + dx = temp; +} + + + + + /* --------------------- Inline and template functions ------------------- */