changeset 13091:e5aaba072d2b

maint: style fixes in iterative linear solvers * bicg.m, bicgstab.m, cgs.m, gmres.m: Style fixes.
author Carlo de Falco <kingcrimson@tiscali.it>
date Sun, 04 Sep 2011 01:21:43 +0200
parents 7f127e079a7c
children 186c3b80ba54
files scripts/sparse/bicg.m scripts/sparse/bicgstab.m scripts/sparse/gmres.m
diffstat 3 files changed, 104 insertions(+), 107 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/sparse/bicg.m
+++ b/scripts/sparse/bicg.m
@@ -19,46 +19,51 @@
 ## @deftypefn {Function File} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
 ## @deftypefnx {Function File} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicg (@var{A}, @var{b}, ...)
-##
-##   Solve @code{A x = b} using the Bi-conjugate gradient iterative method.
+## Solve @code{A x = b} using the Bi-conjugate gradient iterative method.
 ##
-##   @itemize @minus
-##   @item @var{rtol} is the relative tolerance, if not given or set to [] the default value 1e-6 is used.
-##   @item @var{maxit} the maximum number of outer iterations, if not given or set to [] the default value @code{min (20, numel (b))} is used.
-##   @item @var{x0} the initial guess, if not given or set to [] the default value @code{zeros (size (b))} is used. 
-##   @end itemize
+## @itemize @minus
+## @item @var{rtol} is the relative tolerance, if not given 
+## or set to [] the default value 1e-6 is used.
+## @item @var{maxit} the maximum number of outer iterations, 
+## if not given or set to [] the default value 
+## @code{min (20, numel (b))} is used.
+## @item @var{x0} the initial guess, if not given or set to [] 
+## the default value @code{zeros (size (b))} is used. 
+## @end itemize
 ##
-##   @var{A} can be passed as a matrix or as a function handle or 
-##   inline function @code{f} such that @code{f(x, "notransp") = A*x} and @code{f(x, "transp") = A'*x}.
-##
-##   The preconditioner @var{P} is given as @code{P = M1 * M2}. 
-##   Both @var{M1} and @var{M2} can be passed as a matrix or as a function handle or 
-##   inline function @code{g} such that @code{g(x, 'notransp') = M1 \ x} or @code{g(x, 'notransp') = M2 \ x} and 
-##   @code{g(x, 'transp') = M1' \ x} or @code{g(x, 'transp') = M2' \ x}.
+## @var{A} can be passed as a matrix or as a function handle or 
+## inline function @code{f} such that @code{f(x, "notransp") = A*x}
+## and @code{f(x, "transp") = A'*x}.
 ##
-##   If colled with more than one output parameter
+## The preconditioner @var{P} is given as @code{P = M1 * M2}. 
+## Both @var{M1} and @var{M2} can be passed as a matrix or as 
+## a function handle or inline function @code{g} such that 
+## @code{g(x, 'notransp') = M1 \ x} or @code{g(x, 'notransp') = M2 \ x} and 
+## @code{g(x, 'transp') = M1' \ x} or @code{g(x, 'transp') = M2' \ x}.
+##
+## If colled with more than one output parameter
 ##
-##   @itemize @minus
-##   @item @var{flag} indicates the exit status:
-##   @itemize @minus
-##     @item 0: iteration converged to the within the chosen tolerance
-##     @item 1: the maximum number of iterations was reached before convergence
-##     @item 3: the algorithm reached stagnation
-##   @end itemize
-##   (the value 2 is unused but skipped for compatibility).
-##   @item @var{relres} is the final value of the relative residual.
-##   @item @var{iter} is the number of iterations performed. 
-##   @item @var{resvec} is a vector containing the relative residual at each iteration.
-##   @end itemize
+## @itemize @minus
+## @item @var{flag} indicates the exit status:
+## @itemize @minus
+## @item 0: iteration converged to the within the chosen tolerance
+## @item 1: the maximum number of iterations was reached before convergence
+## @item 3: the algorithm reached stagnation
+## @end itemize
+## (the value 2 is unused but skipped for compatibility).
+## @item @var{relres} is the final value of the relative residual.
+## @item @var{iter} is the number of iterations performed. 
+## @item @var{resvec} is a vector containing the relative residual at each iteration.
+## @end itemize
 ##
-##   @seealso{pcg,cgs,bicgstab,gmres}
+## @seealso{bicgstab,cgs,gmres,pcg}
 ##
 ## @end deftypefn
 
 
 function [x, flag, res1, k, resvec] = bicg (A, b, tol, maxit, M1, M2, x0)
 
-  if ((nargin >= 2) && isvector (full (b)))
+  if (nargin >= 2 && isvector (full (b)))
     
     if (ischar (A))
       fun = str2func (A);
@@ -75,15 +80,15 @@
               "be a function or a square matrix"]);
     endif
 
-    if ((nargin < 3) || (isempty (tol)))
+    if (nargin < 3 || isempty (tol))
       tol = 1e-6;
     endif
 
-    if ((nargin < 4) || (isempty (maxit)))
+    if (nargin < 4 || isempty (maxit))
       maxit = min (rows (b), 20);
     endif
 
-    if ((nargin < 5) || isempty (M1))
+    if (nargin < 5 || isempty (M1))
       M1m1x = @(x, ignore) x;
       M1tm1x = M1m1x;
     elseif (ischar (M1))
@@ -101,7 +106,7 @@
               "be a function or matrix"]);
     endif
     
-    if ((nargin < 6) || isempty (M2))
+    if (nargin < 6 || isempty (M2))
       M2m1x = @(x, ignore) x;
       M2tm1x = M2m1x;
     elseif (ischar (M2))
@@ -122,7 +127,7 @@
     Pm1x  = @(x) M2m1x  (M1m1x (x));
     Ptm1x = @(x) M1tm1x (M2tm1x (x));
 
-    if ((nargin < 7) || (isempty (x0)))
+    if (nargin < 7 || isempty (x0))
       x0 = zeros (size (b));
     endif
 
--- a/scripts/sparse/bicgstab.m
+++ b/scripts/sparse/bicgstab.m
@@ -22,53 +22,48 @@
 ## @deftypefn {Function File} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
 ## @deftypefnx {Function File} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicgstab (@var{A}, @var{b}, ...)
-##
-##   Solve @code{A x = b} using the stabilizied Bi-conjugate gradient iterative method.
-##
-##   @itemize @minus
-##
-##   @item @var{rtol} is the relative tolerance, if not given or set to
-##   [] the default value 1e-6 is used.
+## Solve @code{A x = b} using the stabilizied Bi-conjugate gradient iterative method.
 ##
-##   @item @var{maxit} the maximum number of outer iterations, if not
-##   given or set to [] the default value @code{min (20, numel (b))} is
-##   used.
+## @itemize @minus
+## @item @var{rtol} is the relative tolerance, if not given or set to
+## [] the default value 1e-6 is used.
+## @item @var{maxit} the maximum number of outer iterations, if not
+## given or set to [] the default value @code{min (20, numel (b))} is
+## used.
+## @item @var{x0} the initial guess, if not given or set to [] the
+## default value @code{zeros (size (b))} is used.
+## @end itemize
 ##
-##   @item @var{x0} the initial guess, if not given or set to [] the
-##   default value @code{zeros (size (b))} is used.
-##
-##   @end itemize
-##
-##   @var{A} can be passed as a matrix or as a function handle or 
-##   inline function @code{f} such that @code{f(x) = A*x}.
+## @var{A} can be passed as a matrix or as a function handle or 
+## inline function @code{f} such that @code{f(x) = A*x}.
 ##
-##   The preconditioner @var{P} is given as @code{P = M1 * M2}. 
-##   Both @var{M1} and @var{M2} can be passed as a matrix or as a function handle or 
-##   inline function @code{g} such that @code{g(x) = M1 \ x} or @code{g(x) = M2 \ x}.
+## The preconditioner @var{P} is given as @code{P = M1 * M2}. 
+## Both @var{M1} and @var{M2} can be passed as a matrix or as a function handle or 
+## inline function @code{g} such that @code{g(x) = M1 \ x} or @code{g(x) = M2 \ x}.
 ##
-##   If called with more than one output parameter
+## If called with more than one output parameter
 ##
-##   @itemize @minus
-##   @item @var{flag} indicates the exit status:
-##   @itemize @minus
-##     @item 0: iteration converged to the within the chosen tolerance
-##     @item 1: the maximum number of iterations was reached before convergence
-##     @item 3: the algorithm reached stagnation
-##   @end itemize
-##   (the value 2 is unused but skipped for compatibility).
-##   @item @var{relres} is the final value of the relative residual.
-##   @item @var{iter} is the number of iterations performed. 
-##   @item @var{resvec} is a vector containing the relative residual at each iteration.
-##   @end itemize
+## @itemize @minus
+## @item @var{flag} indicates the exit status:
+## @itemize @minus
+## @item 0: iteration converged to the within the chosen tolerance
+## @item 1: the maximum number of iterations was reached before convergence
+## @item 3: the algorithm reached stagnation
+## @end itemize
+## (the value 2 is unused but skipped for compatibility).
+## @item @var{relres} is the final value of the relative residual.
+## @item @var{iter} is the number of iterations performed. 
+## @item @var{resvec} is a vector containing the relative residual at each iteration.
+## @end itemize
 ##
-##   @seealso{pcg,cgs,bicg,gmres}
+## @seealso{bicg,cgs,gmres,pcg}
 ##
 ## @end deftypefn
 
 function [x, flag, relres, iter, resvec] = bicgstab (A, b, tol, maxit, 
                                                      M1, M2, x0)
 
-  if ((nargin >= 2) && (nargin <= 7) && isvector (full (b)))
+  if (nargin >= 2 && nargin <= 7 && isvector (full (b)))
     
     if (ischar (A))
       A = str2func (A);
@@ -81,47 +76,46 @@
               "to be a function or a square matrix"]);
     endif
     
-    if ((nargin < 3) || (isempty (tol)))
+    if (nargin < 3 || isempty (tol))
       tol = 1e-6;
     endif
 
-    if ((nargin < 4) || (isempty (maxit)))
+    if (nargin < 4 || isempty (maxit))
       maxit = min (rows (b), 20);
     endif
 
-    if ((nargin < 5) || isempty (M1))
+    if (nargin < 5 || isempty (M1))
       M1m1x = @(x) x;
     elseif (ischar (M1))
       M1m1x = str2func (M1);
     elseif (ismatrix (M1))
-      M1m1x  = @(x) M1  \ x;
+      M1m1x = @(x) M1  \ x;
     elseif (isa (M1, "function_handle"))
-      M1m1x  = @(x) feval (M1, x);
+      M1m1x = @(x) feval (M1, x);
     else
       error (["bicgstab: preconditioner is " ...
               "expected to be a function or matrix"]);
     endif
     
-    if ((nargin < 6) || isempty (M2))
+    if (nargin < 6 || isempty (M2))
       M2m1x = @(x) x;
     elseif (ischar (M2))
       M2m1x = str2func (M2);
     elseif (ismatrix (M2))
-      M2m1x  = @(x) M2  \ x;
+      M2m1x = @(x) M2  \ x;
     elseif (isa (M2, "function_handle"))
-      M2m1x  = @(x) feval (M2, x);
+      M2m1x = @(x) feval (M2, x);
     else
       error (["bicgstab: preconditioner is "...
               "expected to be a function or matrix"]);
     endif
 
-    precon  = @(x) M2m1x (M1m1x (x));
+    precon = @(x) M2m1x (M1m1x (x));
 
-    if ((nargin < 7) || (isempty (x0)))
+    if (nargin < 7 || isempty (x0))
       x0 = zeros (size (b));
     endif
-
-
+    
     ## specifies initial estimate x0
     if (nargin < 7)
       x = zeros (rows (b), 1);
--- a/scripts/sparse/gmres.m
+++ b/scripts/sparse/gmres.m
@@ -49,30 +49,23 @@
 ##
 ## @itemize @minus
 ## @item @var{flag} indicates the exit status:
-##
 ## @table @asis
-##   @item 0 : iteration converged to within the specified tolerance
-##
-##   @item 1 : maximum number of iterations exceeded
-##
-##   @item 2 : unused, but skipped for compatibility
-##
-##   @item 3 : algorithm reached stagnation
+## @item 0 : iteration converged to within the specified tolerance
+## @item 1 : maximum number of iterations exceeded
+## @item 2 : unused, but skipped for compatibility
+## @item 3 : algorithm reached stagnation
 ## @end table
-##
 ## @item @var{relres} is the final value of the relative residual.
-##
 ## @item @var{iter} is a vector containing the number of outer iterations and
 ## total iterations performed.
-##
 ## @item @var{resvec} is a vector containing the relative residual at each
 ## iteration.
 ## @end itemize
 ##
-## @seealso{pcg, cgs, bicgstab}
+## @seealso{bicg, bicgstab, cgs, pcg}
 ## @end deftypefn
 
-function [x, flag, prec_res_norm, itcnt] = gmres (A, b, restart, rtol, maxit, M1, M2, x0)
+function [x, flag, presn, it] = gmres (A, b, restart, rtol, maxit, M1, M2, x0)
 
   if (nargin < 2 || nargin > 8)
     print_usage ();    
@@ -133,7 +126,7 @@
   x_old = x0; 
   x = x_old;
   prec_res = Pm1x (b - Ax (x_old));
-  prec_res_norm = norm (prec_res, 2);
+  presn = norm (prec_res, 2);
   
   B = zeros (restart + 1, 1);
   V = zeros (rows (x), restart);
@@ -143,49 +136,54 @@
   iter = 1;
   restart_it  = restart + 1; 
   resids      = zeros (maxit, 1);
-  resids(1)   = prec_res_norm;
+  resids(1)   = presn;
   prec_b_norm = norm (Pm1x (b), 2);
   flag        = 1;
 
-  while ((iter <= maxit * restart) && (prec_res_norm > rtol * prec_b_norm))
+  while (iter <= maxit * restart && presn > rtol * prec_b_norm)
   
     ## restart
     if (restart_it > restart)
       restart_it = 1;
       x_old = x;	      
       prec_res = Pm1x (b - Ax (x_old));
-      prec_res_norm = norm (prec_res, 2);
-      B(1) = prec_res_norm;
+      presn = norm (prec_res, 2);
+      B(1) = presn;
       H(:) = 0;
-      V(:, 1) = prec_res / prec_res_norm;
+      V(:, 1) = prec_res / presn;
     endif  
     
     ## basic iteration
     tmp = Pm1x (Ax (V(:, restart_it)));
-    [V(:,restart_it+1), H(1:restart_it+1, restart_it)] = mgorth (tmp, V(:,1:restart_it));
+    [V(:,restart_it+1), H(1:restart_it+1, restart_it)] = ...
+        mgorth (tmp, V(:,1:restart_it));
     
     Y = (H(1:restart_it+1, 1:restart_it) \ B (1:restart_it+1));
 	      
-    little_res = B(1:restart_it+1) - H(1:restart_it+1, 1:restart_it) * Y(1:restart_it);
-    prec_res_norm = norm (little_res, 2);
+    little_res = B(1:restart_it+1) - ...
+        H(1:restart_it+1, 1:restart_it) * Y(1:restart_it);
+
+    presn = norm (little_res, 2);
 	      
     x = x_old + V(:, 1:restart_it) * Y(1:restart_it);
     
-    resids(iter) = prec_res_norm ;
+    resids(iter) = presn;
     if (norm (x - x_old, inf) <= eps)
       flag = 3;
       break
     endif
 
-    restart_it++ ; iter++;
+    restart_it++ ; 
+    iter++;
   endwhile
 
-  if (prec_res_norm > rtol * prec_b_norm)
+  if (presn > rtol * prec_b_norm)
     flag = 0;
   endif
 
   resids = resids(1:iter-1);
-  itcnt = [floor(maxit/restart), rem(maxit, restart)];
+  it = [floor(maxit/restart), rem(maxit, restart)];
+
 endfunction