Mercurial > octave-nkf
changeset 13073:51bc892d5cf8
Move bicg.m to scripts/sparse for uniformity
* bicg.m : Move to scripts/sparse where the other iterative
solvers are.
| author | Carlo de Falco <kingcrimson@tiscali.it> |
|---|---|
| date | Sat, 03 Sep 2011 21:51:26 +0200 |
| parents | 39814ef3278b |
| children | 0b789a03bde1 |
| files | scripts/linear-algebra/bicg.m scripts/linear-algebra/module.mk scripts/sparse/bicg.m scripts/sparse/module.mk |
| diffstat | 4 files changed, 250 insertions(+), 250 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/linear-algebra/bicg.m Sat Sep 03 12:44:15 2011 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,249 +0,0 @@ -## Copyright (C) 2006 Sylvain Pelissier <sylvain.pelissier@gmail.com> -## Copyright (C) 2011 Carlo de Falco -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; If not, see <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## -## @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. -## -## @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}. -## -## 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 -## -## @seealso{pcg,cgs,bicgstab,gmres} -## -## @end deftypefn - - -function [x, flag, res1, k, resvec] = bicg (A, b, tol, maxit, M1, M2, x0) - - if ((nargin >= 2) && isvector (full (b))) - - if (ischar (A)) - fun = str2func (A); - Ax = @(x) feval (fun, x, "notransp"); - Atx = @(x) feval (fun, x, "transp"); - elseif (ismatrix (A)) - Ax = @(x) A * x; - Atx = @(x) A' * x; - elseif (isa (A, "function_handle")) - Ax = @(x) feval (A, x, "notransp"); - Atx = @(x) feval (A, x, "transp"); - else - error (["bicg: first argument is expected to " ... - "be a function or a square matrix"]); - endif - - if ((nargin < 3) || (isempty (tol))) - tol = 1e-6; - endif - - if ((nargin < 4) || (isempty (maxit))) - maxit = min (rows (b), 20); - endif - - if ((nargin < 5) || isempty (M1)) - M1m1x = @(x, ignore) x; - M1tm1x = M1m1x; - elseif (ischar (M1)) - fun = str2func (M1); - M1m1x = @(x) feval (fun, x, "notransp"); - M1tm1x = @(x) feval (fun, x, "transp"); - elseif (ismatrix (M1)) - M1m1x = @(x) M1 \ x; - M1tm1x = @(x) M1' \ x; - elseif (isa (M1, "function_handle")) - M1m1x = @(x) feval (M1, x, "notransp"); - M1tm1x = @(x) feval (M1, x, "transp"); - else - error (["bicg: preconditioner is expected to " ... - "be a function or matrix"]); - endif - - if ((nargin < 6) || isempty (M2)) - M2m1x = @(x, ignore) x; - M2tm1x = M2m1x; - elseif (ischar (M2)) - fun = str2func (M2); - M2m1x = @(x) feval (fun, x, "notransp"); - M2tm1x = @(x) feval (fun, x, "transp"); - elseif (ismatrix (M2)) - M2m1x = @(x) M2 \ x; - M2tm1x = @(x) M2' \ x; - elseif (isa (M2, "function_handle")) - M2m1x = @(x) feval (M2, x, "notransp"); - M2tm1x = @(x) feval (M2, x, "transp"); - else - error (["bicg: preconditioner is expected to " ... - "be a function or matrix"]); - endif - - Pm1x = @(x) M2m1x (M1m1x (x)); - Ptm1x = @(x) M1tm1x (M2tm1x (x)); - - if ((nargin < 7) || (isempty (x0))) - x0 = zeros (size (b)); - endif - - y = x = x0; - c = b; - - r0 = b - Ax (x); - s0 = c - Atx (y); - - d = Pm1x (r0); - f = Ptm1x (s0); - - bnorm = norm (b); - res0 = Inf; - - if (any (r0 != 0)) - - for k = 1:maxit - - a = (s0' * Pm1x (r0)) ./ (f' * Ax (d)); - - x += a * d; - y += conj (a) * f; - - r1 = r0 - a * Ax (d); - s1 = s0 - conj (a) * Atx (f); - - beta = (s1' * Pm1x (r1)) ./ (s0' * Pm1x (r0)); - - d = Pm1x (r1) + beta * d; - f = Ptm1x (s1) + conj (beta) * f; - - r0 = r1; - s0 = s1; - - res1 = norm (b - Ax (x)) / bnorm; - if (res1 < tol) - flag = 0; - if (nargout < 2) - printf ("bicg converged at iteration %i ", k); - printf ("to a solution with relative residual %e\n", res1); - endif - break; - endif - - if (res0 <= res1) - flag = 3; - printf ("bicg stopped at iteration %i ", k); - printf ("without converging to the desired tolerance %e\n", tol); - printf ("because the method stagnated.\n"); - printf ("The iterate returned (number %i) ", k-1); - printf ("has relative residual %e\n", res0); - break - endif - res0 = res1; - if (nargout > 4) - resvec(k) = res0; - endif - endfor - - if (k == maxit) - flag = 1; - printf ("bicg stopped at iteration %i ", maxit); - printf ("without converging to the desired tolerance %e\n", tol); - printf ("because the maximum number of iterations was reached. "); - printf ("The iterate returned (number %i) has ", maxit); - printf ("relative residual %e\n", res1); - endif - - else - flag = 0; - if (nargout < 2) - printf ("bicg converged after 0 interations\n"); - endif - endif - - else - print_usage (); - endif - -endfunction; - - -%!test -%! n = 100; -%! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n); -%! b = sum (A, 2); -%! tol = 1e-8; -%! maxit = 15; -%! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n); -%! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n); -%! [x, flag, relres, iter, resvec] = bicg (A, b, tol, maxit, M1, M2); -%! assert (x, ones (size (b)), 1e-7); -%! -%!test -%! n = 100; -%! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n); -%! b = sum (A, 2); -%! tol = 1e-8; -%! maxit = 15; -%! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n); -%! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n); -%! -%! function y = afun (x, t, a) -%! switch t -%! case "notransp" -%! y = a * x; -%! case "transp" -%! y = a' * x; -%! endswitch -%! endfunction -%! -%! [x, flag, relres, iter, resvec] = bicg (@(x, t) afun (x, t, A), -%! b, tol, maxit, M1, M2); -%! assert (x, ones (size (b)), 1e-7); - -%!test -%! n = 100; -%! tol = 1e-8; -%! a = sprand (n, n, .1); -%! A = a' * a + 100 * eye (n); -%! b = sum (A, 2); -%! [x, flag, relres, iter, resvec] = bicg (A, b, tol, [], diag (diag (A))); -%! assert (x, ones (size (b)), 1e-7);
--- a/scripts/linear-algebra/module.mk Sat Sep 03 12:44:15 2011 -0700 +++ b/scripts/linear-algebra/module.mk Sat Sep 03 21:51:26 2011 +0200 @@ -1,7 +1,6 @@ FCN_FILE_DIRS += linear-algebra linear_algebra_FCN_FILES = \ - linear-algebra/bicg.m \ linear-algebra/commutation_matrix.m \ linear-algebra/cond.m \ linear-algebra/condest.m \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/sparse/bicg.m Sat Sep 03 21:51:26 2011 +0200 @@ -0,0 +1,249 @@ +## Copyright (C) 2006 Sylvain Pelissier <sylvain.pelissier@gmail.com> +## Copyright (C) 2011 Carlo de Falco +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## +## @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. +## +## @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}. +## +## 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 +## +## @seealso{pcg,cgs,bicgstab,gmres} +## +## @end deftypefn + + +function [x, flag, res1, k, resvec] = bicg (A, b, tol, maxit, M1, M2, x0) + + if ((nargin >= 2) && isvector (full (b))) + + if (ischar (A)) + fun = str2func (A); + Ax = @(x) feval (fun, x, "notransp"); + Atx = @(x) feval (fun, x, "transp"); + elseif (ismatrix (A)) + Ax = @(x) A * x; + Atx = @(x) A' * x; + elseif (isa (A, "function_handle")) + Ax = @(x) feval (A, x, "notransp"); + Atx = @(x) feval (A, x, "transp"); + else + error (["bicg: first argument is expected to " ... + "be a function or a square matrix"]); + endif + + if ((nargin < 3) || (isempty (tol))) + tol = 1e-6; + endif + + if ((nargin < 4) || (isempty (maxit))) + maxit = min (rows (b), 20); + endif + + if ((nargin < 5) || isempty (M1)) + M1m1x = @(x, ignore) x; + M1tm1x = M1m1x; + elseif (ischar (M1)) + fun = str2func (M1); + M1m1x = @(x) feval (fun, x, "notransp"); + M1tm1x = @(x) feval (fun, x, "transp"); + elseif (ismatrix (M1)) + M1m1x = @(x) M1 \ x; + M1tm1x = @(x) M1' \ x; + elseif (isa (M1, "function_handle")) + M1m1x = @(x) feval (M1, x, "notransp"); + M1tm1x = @(x) feval (M1, x, "transp"); + else + error (["bicg: preconditioner is expected to " ... + "be a function or matrix"]); + endif + + if ((nargin < 6) || isempty (M2)) + M2m1x = @(x, ignore) x; + M2tm1x = M2m1x; + elseif (ischar (M2)) + fun = str2func (M2); + M2m1x = @(x) feval (fun, x, "notransp"); + M2tm1x = @(x) feval (fun, x, "transp"); + elseif (ismatrix (M2)) + M2m1x = @(x) M2 \ x; + M2tm1x = @(x) M2' \ x; + elseif (isa (M2, "function_handle")) + M2m1x = @(x) feval (M2, x, "notransp"); + M2tm1x = @(x) feval (M2, x, "transp"); + else + error (["bicg: preconditioner is expected to " ... + "be a function or matrix"]); + endif + + Pm1x = @(x) M2m1x (M1m1x (x)); + Ptm1x = @(x) M1tm1x (M2tm1x (x)); + + if ((nargin < 7) || (isempty (x0))) + x0 = zeros (size (b)); + endif + + y = x = x0; + c = b; + + r0 = b - Ax (x); + s0 = c - Atx (y); + + d = Pm1x (r0); + f = Ptm1x (s0); + + bnorm = norm (b); + res0 = Inf; + + if (any (r0 != 0)) + + for k = 1:maxit + + a = (s0' * Pm1x (r0)) ./ (f' * Ax (d)); + + x += a * d; + y += conj (a) * f; + + r1 = r0 - a * Ax (d); + s1 = s0 - conj (a) * Atx (f); + + beta = (s1' * Pm1x (r1)) ./ (s0' * Pm1x (r0)); + + d = Pm1x (r1) + beta * d; + f = Ptm1x (s1) + conj (beta) * f; + + r0 = r1; + s0 = s1; + + res1 = norm (b - Ax (x)) / bnorm; + if (res1 < tol) + flag = 0; + if (nargout < 2) + printf ("bicg converged at iteration %i ", k); + printf ("to a solution with relative residual %e\n", res1); + endif + break; + endif + + if (res0 <= res1) + flag = 3; + printf ("bicg stopped at iteration %i ", k); + printf ("without converging to the desired tolerance %e\n", tol); + printf ("because the method stagnated.\n"); + printf ("The iterate returned (number %i) ", k-1); + printf ("has relative residual %e\n", res0); + break + endif + res0 = res1; + if (nargout > 4) + resvec(k) = res0; + endif + endfor + + if (k == maxit) + flag = 1; + printf ("bicg stopped at iteration %i ", maxit); + printf ("without converging to the desired tolerance %e\n", tol); + printf ("because the maximum number of iterations was reached. "); + printf ("The iterate returned (number %i) has ", maxit); + printf ("relative residual %e\n", res1); + endif + + else + flag = 0; + if (nargout < 2) + printf ("bicg converged after 0 interations\n"); + endif + endif + + else + print_usage (); + endif + +endfunction; + + +%!test +%! n = 100; +%! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n); +%! b = sum (A, 2); +%! tol = 1e-8; +%! maxit = 15; +%! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n); +%! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n); +%! [x, flag, relres, iter, resvec] = bicg (A, b, tol, maxit, M1, M2); +%! assert (x, ones (size (b)), 1e-7); +%! +%!test +%! n = 100; +%! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n); +%! b = sum (A, 2); +%! tol = 1e-8; +%! maxit = 15; +%! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n); +%! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n); +%! +%! function y = afun (x, t, a) +%! switch t +%! case "notransp" +%! y = a * x; +%! case "transp" +%! y = a' * x; +%! endswitch +%! endfunction +%! +%! [x, flag, relres, iter, resvec] = bicg (@(x, t) afun (x, t, A), +%! b, tol, maxit, M1, M2); +%! assert (x, ones (size (b)), 1e-7); + +%!test +%! n = 100; +%! tol = 1e-8; +%! a = sprand (n, n, .1); +%! A = a' * a + 100 * eye (n); +%! b = sum (A, 2); +%! [x, flag, relres, iter, resvec] = bicg (A, b, tol, [], diag (diag (A))); +%! assert (x, ones (size (b)), 1e-7);