From 4b13578c5e6c270398e2cd33cf1f024720028960 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 17 May 1999 13:30:47 +0000 Subject: [PATCH] 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@1332 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/solver_gmres.h | 130 ------------------------- 1 file changed, 130 deletions(-) diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h index 4d47b06bb3..c5cd615dba 100644 --- a/deal.II/lac/include/lac/solver_gmres.h +++ b/deal.II/lac/include/lac/solver_gmres.h @@ -136,136 +136,6 @@ 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 ------------------- */ -- 2.39.5