# HG changeset patch
# User Radek Salac <salac.r@gmail.com>
# Date 1243486991 -7200
# Node ID 2b35cb600d5038a9ec5e6229eb1f8b70cd4e8f12
# Parent  be84e9654febd145973d2951c93641de5b0bc176
improve bicgstab and cgs

diff -r be84e9654feb -r 2b35cb600d50 scripts/ChangeLog
--- a/scripts/ChangeLog	Wed May 27 13:06:38 2009 -0700
+++ b/scripts/ChangeLog	Thu May 28 07:03:11 2009 +0200
@@ -1,3 +1,12 @@
+2009-05-28  Radek Salac  <salac.r@gmail.com>
+
+	* sparse/bicgstab.m: New output when calling without arguments.
+	Time optimization - remove certain checks linked to preconditioner which
+	are unacceptably slow.
+	* sparse/cgs.m: New output when calling without arguments.
+	Time optimization - remove certain checks linked to preconditioner which
+	are unacceptably slow. Rename internal variable to match bicgstab.
+
 2009-05-27  Rik Wehbring  <rdrider0-list@yahoo.com>
 
 	* plot/axis.m: Update documentation to reflect addition of "tight" option.
diff -r be84e9654feb -r 2b35cb600d50 scripts/sparse/bicgstab.m
--- a/scripts/sparse/bicgstab.m	Wed May 27 13:06:38 2009 -0700
+++ b/scripts/sparse/bicgstab.m	Thu May 28 07:03:11 2009 +0200
@@ -1,4 +1,4 @@
-## Copyright (C) 2008, 2009 Radek Salac
+## Copyright (C) 2008 Radek Salac
 ##
 ## This file is part of Octave.
 ##
@@ -84,20 +84,10 @@
     endif 
   endif
 
-  if (nargin > 4 && isnumeric (precon))
-    ## Precon can by also function.
-    if (det (precon) != 0) 
-      ## We can compute inverse preconditioner and use quicker algorithm.
-      precon = inv (precon);
-    else
-      error ("bicgstab: preconditioner is ill conditioned");
-    endif
+  if (nargin > 4 && isnumeric(precon) )
+    precon = inv(precon);
+  endif
 
-    if (isinf (cond (precon))); 
-      ## We must make test if preconditioner isn't ill conditioned.
-      error ("bicgstab: preconditioner is ill conditioned");
-    endif
-  endif
 
   ## specifies initial estimate x0
   if (nargin < 7)
@@ -112,7 +102,7 @@
   rr = res;
 
   ## Vector of the residual norms for each iteration.
-  resvec = [norm(res)];
+  resvec = [norm(res)/norm_b];
 
   ## Default behaviour we don't reach tolerance tol within maxit iterations.
   flag = 1;
@@ -169,6 +159,23 @@
     endif
   endfor
 
+  if (nargout < 2)
+    if (flag == 0) 
+      printf (["bicgstab converged at iteration %i ",
+      "to a solution with relative residual %e\n"],iter,relres);
+    elseif (flag == 3)
+      printf (["bicgstab 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 (["bicgstab stopped at iteration %i ",
+      "without converging to the desired toleranc %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
 
 %!demo
diff -r be84e9654feb -r 2b35cb600d50 scripts/sparse/cgs.m
--- a/scripts/sparse/cgs.m	Wed May 27 13:06:38 2009 -0700
+++ b/scripts/sparse/cgs.m	Thu May 28 07:03:11 2009 +0200
@@ -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