Mercurial > hg > octave-nkf
changeset 12655:52d79740a76c
maint: Remove linear-algebra/gmres.m which was moved to sparse/ directory.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Sat, 07 May 2011 14:30:20 -0700 |
parents | 25d4e163ee38 |
children | 6b2f14af2360 |
files | scripts/linear-algebra/gmres.m |
diffstat | 1 files changed, 0 insertions(+), 214 deletions(-) [+] |
line wrap: on
line diff
deleted file mode 100644 --- a/scripts/linear-algebra/gmres.m +++ /dev/null @@ -1,214 +0,0 @@ -## 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));