Mercurial > hg > octave-nkf
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