--- /dev/null
+## Copyright (C) 2009-2012 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{bicg, bicgstab, cgs, pcg}
+## @end deftypefn
+
+function [x, flag, presn, it, resids] = gmres (A, b, restart, rtol, maxit, M1, M2, x0)
+
+ if (nargin < 2 || nargin > 8)
+ print_usage ();
+ endif
+
+ 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));
+ presn = 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) = presn;
+ prec_b_norm = norm (Pm1x (b), 2);
+ flag = 1;
+
+ 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));
+ presn = norm (prec_res, 2);
+ B(1) = presn;
+ H(:) = 0;
+ 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));
+
+ 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);
+
+ presn = norm (little_res, 2);
+
+ x = x_old + V(:, 1:restart_it) * Y(1:restart_it);
+
+ resids(iter) = presn;
+ if (norm (x - x_old, inf) <= eps)
+ flag = 3;
+ break
+ endif
+
+ restart_it++ ;
+ iter++;
+ endwhile
+
+ if (presn > rtol * prec_b_norm)
+ flag = 0;
+ endif
+
+ resids = resids(1:iter-1);
+ it = [ceil(iter / restart), rem(iter, 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));
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff