diff scripts/sparse/cgs.m @ 9278:2b35cb600d50

improve bicgstab and cgs
author Radek Salac <salac.r@gmail.com>
date Thu, 28 May 2009 07:03:11 +0200
parents eb63fbe60fab
children 1673a0dc019f
line wrap: on
line diff
--- a/scripts/sparse/cgs.m
+++ b/scripts/sparse/cgs.m
@@ -1,4 +1,4 @@
-## Copyright (C) 2008, 2009 Radek Salac
+## Copyright (C) 2008 Radek Salac
 ##
 ## This file is part of Octave.
 ##
@@ -76,20 +76,8 @@
     endif
   endif
 
-  ## Precon can also be a function.
-  if (nargin > 4 && isnumeric (precon))
-
-    ## We can compute inverse preconditioner and use quicker algorithm.
-    if (det (precon) != 0)
-     precon = inv (precon);
-    else
-     error ("cgs: preconditioner is ill conditioned");
-    endif
-
-    ## We must make test if preconditioner isn't ill conditioned.
-    if (isinf (cond (precon))); 
-      error ("cgs: preconditioner is ill conditioned");
-    endif
+  if (nargin > 4 && isnumeric(precon) )
+    precon = inv(precon);
   endif
 
   ## Specifies initial estimate x0.
@@ -99,29 +87,28 @@
     x = x0;
   endif
 
-  relres = b - A * x;
+  res = b - A * x;
+  norm_b = norm (b);
   ## Vector of the residual norms for each iteration.
-  resvec = [norm(relres)];
+  resvec = [ norm(res)/norm_b ];
   ro = 0;
-  norm_b = norm (b);
   ## Default behavior we don't reach tolerance tol within maxit iterations.
   flag = 1;
   for iter = 1 : maxit
 
     if (nargin > 4 && isnumeric (precon))
-      ## We have computed inverse matrix so we can use quick algorithm.
-      z = precon * relres;
+      z = precon * res;
     elseif (nargin > 4)
       ## Our preconditioner is a function.
-      z = feval (precon, relres);
+      z = feval (precon, res);
     else
       ## We don't use preconditioning.
-      z = relres;
+      z = res;
     endif
 
     ## Cache.
     ro_old = ro;
-    ro = relres' * z;
+    ro = res' * z;
     if (iter == 1)
       p = z;
     else
@@ -133,11 +120,11 @@
     alpha = ro / (p' * q);
     x = x + alpha * p;
 
-    relres = relres - alpha * q;
-    resvec = [resvec; norm(relres)];
+    res = res - alpha * q;
+    relres = norm(res) / norm_b;
+    resvec = [resvec;relres];
 
-    relres_distance = resvec (end) / norm_b;
-    if (relres_distance <= tol)
+    if (relres <= tol)
       ## We reach tolerance tol within maxit iterations.
       flag = 0;
       break;
@@ -148,7 +135,23 @@
     endif
   endfor;
 
-  relres = relres_distance;
+  if (nargout < 1)
+    if ( flag == 0 ) 
+      printf (["cgs converged at iteration %i ",
+      "to a solution with relative residual %e\n"],iter,relres);
+    elseif (flag == 3)
+      printf (["cgs stopped at iteration %i ",
+      "without converging to the desired tolerance %e\n",
+      "because the method stagnated.\n",
+      "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres);
+    else
+      printf (["cgs stopped at iteration %i ",
+      "without converging to the desired tolerance %e\n",
+      "because the maximum number of iterations was reached.\n",
+      "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres);
+    endif
+  endif
+
 endfunction