Mercurial > octave-nkf
changeset 12650:b7f16997df4d
maint: Move gmres.m to sparse directory rather than linear-algebra directory.
* sparse/module.mk: Update list to include gmres.m
* sparse/gmres.m: Newly moved function
| author | Rik <octave@nomad.inbox5.com> |
|---|---|
| date | Fri, 06 May 2011 09:39:08 -0700 |
| parents | 8579e7a6e6e0 |
| children | 332bb3b9600e |
| files | scripts/sparse/gmres.m scripts/sparse/module.mk |
| diffstat | 2 files changed, 215 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/sparse/gmres.m Fri May 06 09:39:08 2011 -0700 @@ -0,0 +1,214 @@ +## Copyright (C) 2009-2011 Carlo de Falco +## +## This file is part of Octave. +## +## Octave 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 3 of the License, or (at your +## option) any later version. +## +## Octave 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 Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{x} =} gmres (@var{A}, @var{b}, @var{m}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0}) +## @deftypefnx {Function File} {@var{x} =} gmres (@var{A}, @var{b}, @var{m}, @var{rtol}, @var{maxit}, @var{P}) +## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} gmres (@dots{}) +## Solve @code{A x = b} using the Preconditioned GMRES iterative method +## with restart, a.k.a. PGMRES(m). +## +## @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} is the maximum number of outer iterations, +## if not given or set to [] the default value +## @code{min (10, numel (b) / restart)} is used. +## +## @item @var{x0} is the initial guess, +## if not given or set to [] the default value @code{zeros(size (b))} is used. +## +## @item @var{m} is the restart parameter, +## if not given or set to [] the default value @code{numel (b)} is used. +## @end itemize +## +## Argument @var{A} can be passed as a matrix, 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, function handle, or +## inline function @code{g} such that @code{g(x) = M1\x} or @code{g(x) = M2\x}. +## +## Besides the vector @var{x}, additional outputs are: +## +## @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 +## @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, bigcstab} +## @end deftypefn + +function [x, flag, prec_res_norm, itcnt] = gmres (A, b, restart, rtol, maxit, M1, M2, x0) + + if (nargin < 2 || nargin > 8) + print_usage (); + end + + if (ischar (A)) + Ax = str2func (A); + elseif (ismatrix (A)) + Ax = @(x) A*x; + elseif (isa (A, "function_handle")) + Ax = A; + else + error ("gmres: A must be a function or matrix"); + endif + + if (nargin < 3 || isempty (restart)) + restart = rows (b); + endif + + if (nargin < 4 || isempty (rtol)) + rtol = 1e-6; + endif + + if (nargin < 5 || isempty (maxit)) + maxit = min (rows (b)/restart, 10); + endif + + if (nargin < 6 || isempty (M1)) + M1m1x = @(x) x; + elseif (ischar (M1)) + M1m1x = str2func (M1); + elseif (ismatrix (M1)) + M1m1x = @(x) M1 \ x; + elseif (isa (M1, "function_handle")) + M1m1x = M1; + else + error ("gmres: preconditioner M1 must be a function or matrix"); + endif + + if (nargin < 7 || isempty (M2)) + M2m1x = @(x) x; + elseif (ischar (M2)) + M2m1x = str2func (M2); + elseif (ismatrix (M2)) + M2m1x = @(x) M2 \ x; + elseif (isa (M2, "function_handle")) + M2m1x = M2; + else + error ("gmres: preconditioner M2 must be a function or matrix"); + endif + + Pm1x = @(x) M2m1x (M1m1x (x)); + + if (nargin < 8 || isempty (x0)) + x0 = zeros (size (b)); + endif + + x_old = x0; + x = x_old; + prec_res = Pm1x (b - Ax (x_old)); + prec_res_norm = norm (prec_res, 2); + + B = zeros (restart + 1, 1); + V = zeros (rows (x), restart); + H = zeros (restart + 1, restart); + + ## begin loop + iter = 1; + restart_it = restart + 1; + resids = zeros (maxit, 1); + resids(1) = prec_res_norm; + prec_b_norm = norm (Pm1x (b), 2); + flag = 1; + + while ((iter <= maxit * restart) && (prec_res_norm > 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; + H(:) = 0; + V(:, 1) = prec_res / prec_res_norm; + 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)); + + 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); + + x = x_old + V(:, 1:restart_it) * Y(1:restart_it); + + resids(iter) = prec_res_norm ; + if (norm (x - x_old, inf) <= eps) + flag = 3; + break + endif + + restart_it++ ; iter++; + endwhile + + if (prec_res_norm > rtol * prec_b_norm) + flag = 0; + endif + + resids = resids(1:iter-1); + itcnt = [floor(maxit/restart), rem(maxit, restart)]; +endfunction + + +%!shared A, b, dim +%! dim = 100; +%!test +%! A = spdiags ([-ones(dim,1) 2*ones(dim,1) ones(dim,1)], [-1:1], dim, dim); +%! b = ones(dim, 1); +%! x = gmres (A, b, 10, 1e-10, dim, @(x) x./diag(A), [], b); +%! assert(x, A\b, 1e-9*norm(x,inf)); +%! +%!test +%! x = gmres (A, b, dim, 1e-10, 1e4, @(x) diag(diag(A))\x, [], b); +%! assert(x, A\b, 1e-7*norm(x,inf)); +%! +%!test +%! A = spdiags ([[1./(2:2:2*(dim-1)) 0]; 1./(1:2:2*dim-1); [0 1./(2:2:2*(dim-1))]]', -1:1, dim, dim); +%! A = A'*A; +%! b = rand (dim, 1); +%! [x, resids] = gmres (@(x) A*x, b, dim, 1e-10, dim, @(x) x./diag(A), [], []); +%! assert(x, A\b, 1e-9*norm(x,inf)) +%! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) diag(diag(A))\x, [], []); +%! assert(x, A\b, 1e-9*norm(x,inf)); +%!test +%! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) x./diag(A), [], []); +%! assert(x, A\b, 1e-7*norm(x,inf));