# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1243489033 -7200
# Node ID 1673a0dc019f9889e3cf76dd58527dcec4059b5e
# Parent  2b35cb600d5038a9ec5e6229eb1f8b70cd4e8f12
fix & improve preconditioning strategy in cgs & bicgstab

diff -r 2b35cb600d50 -r 1673a0dc019f scripts/ChangeLog
--- a/scripts/ChangeLog	Thu May 28 07:03:11 2009 +0200
+++ b/scripts/ChangeLog	Thu May 28 07:37:13 2009 +0200
@@ -1,3 +1,8 @@
+2009-05-28  Jaroslav Hajek  <highegg@gmail.com>
+
+	* sparse/bicgstab.m: Improve preconditioning; avoid explicit inverse.
+	* sparse/cgs.m: Improve preconditioning; avoid explicit inverse.
+
 2009-05-28  Radek Salac  <salac.r@gmail.com>
 
 	* sparse/bicgstab.m: New output when calling without arguments.
diff -r 2b35cb600d50 -r 1673a0dc019f scripts/sparse/bicgstab.m
--- a/scripts/sparse/bicgstab.m	Thu May 28 07:03:11 2009 +0200
+++ b/scripts/sparse/bicgstab.m	Thu May 28 07:37:13 2009 +0200
@@ -75,20 +75,20 @@
 
   ## Left preconditioner.
   if (nargin == 5)
-    precon = M1;
+    if (isnumeric (M1))
+      precon = @(x) M1 \ x;
+    endif
   elseif (nargin > 5)
-    if (isparse(M1) && issparse(M2))
-      precon = @(x) M1 * (M2 * x);
+    if (issparse (M1) && issparse (M2))
+      precon = @(x) M2 \ (M1 \ x);
     else
-      precon = M1 * M2;
-    endif 
+      M = M1*M2;
+      precon = @(x) M \ x;
+    endif
+  else
+    precon = @(x) x;
   endif
 
-  if (nargin > 4 && isnumeric(precon) )
-    precon = inv(precon);
-  endif
-
-
   ## specifies initial estimate x0
   if (nargin < 7)
     x = zeros (rows (b), 1);
@@ -117,27 +117,13 @@
       p = res + beta * (p - omega * v);
     endif
 
-    if (nargin > 4 && isnumeric (precon))
-      phat = precon * p;
-    elseif (nargin > 4)
-      ## Our preconditioner is a function.
-      phat = feval (precon, p);
-    else
-      phat = p;
-    endif
+    phat = precon (p);
 
     v = A * phat;
     alpha = rho_1 / (rr' * v);
     s = res - alpha * v;
 
-    if (nargin > 4 && isnumeric (precon))
-      shat = precon * s;
-    elseif (nargin > 4)
-      ## Our preconditioner is a function.
-      shat = feval (precon, s);
-    else
-      shat = s;
-    endif
+    shat = precon (s);
 
     t = A * shat; 
     omega = (t' * s) / (t' * t);
diff -r 2b35cb600d50 -r 1673a0dc019f scripts/sparse/cgs.m
--- a/scripts/sparse/cgs.m	Thu May 28 07:03:11 2009 +0200
+++ b/scripts/sparse/cgs.m	Thu May 28 07:37:13 2009 +0200
@@ -65,19 +65,19 @@
   endif
 
   ## Left preconditioner.
-  precon = [];
   if (nargin == 5)
-    precon = M1;
+    if (isnumeric (M1))
+      precon = @(x) M1 \ x;
+    endif
   elseif (nargin > 5)
-    if (isparse(M1) && issparse(M2))
-      precon = @(x) M1 * (M2 * x);
+    if (issparse (M1) && issparse (M2))
+      precon = @(x) M2 \ (M1 \ x);
     else
-      precon = M1 * M2;
+      M = M1*M2;
+      precon = @(x) M \ x;
     endif
-  endif
-
-  if (nargin > 4 && isnumeric(precon) )
-    precon = inv(precon);
+  else
+    precon = @(x) x;
   endif
 
   ## Specifies initial estimate x0.
@@ -96,15 +96,7 @@
   flag = 1;
   for iter = 1 : maxit
 
-    if (nargin > 4 && isnumeric (precon))
-      z = precon * res;
-    elseif (nargin > 4)
-      ## Our preconditioner is a function.
-      z = feval (precon, res);
-    else
-      ## We don't use preconditioning.
-      z = res;
-    endif
+    z = precon (res);
 
     ## Cache.
     ro_old = ro;