diff scripts/sparse/bicgstab.m @ 9279:1673a0dc019f

fix & improve preconditioning strategy in cgs & bicgstab
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 28 May 2009 07:37:13 +0200
parents 2b35cb600d50
children be55736a0783
line wrap: on
line diff
--- 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);