]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Commit this to have it in the CVS, but I will delete it immediately after again since...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 17 May 1999 13:30:01 +0000 (13:30 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 17 May 1999 13:30:01 +0000 (13:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@1331 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/solver_gmres.h

index c5cd615dba3aa141b23163500e72f9c053c7d475..4d47b06bb37c1758652ea8b874a67e2e62054ffd 100644 (file)
@@ -136,6 +136,136 @@ class SolverGMRES : public Solver<Matrix, Vector>
 
 
 
+template <class Vector >
+void 
+Update(::Vector<double> &x, int k, FullMatrix<double> &h,
+       ::Vector<double> &s, vector<Vector> &v)
+{
+  ::Vector<double> 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<double> H(m+1,m);
+  
+  double resid;
+  int i, j = 1, k;
+  ::Vector<double> 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<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<class Real> 
+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<class Real> 
+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 ------------------- */
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.