octave_packages/m/sparse/gmres.m

   1 ## Copyright (C) 2009-2012 Carlo de Falco
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by the
   7 ## Free Software Foundation; either version 3 of the License, or (at your
   8 ## option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but WITHOUT
  11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12 ## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13 ## for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {@var{x} =} gmres (@var{A}, @var{b}, @var{m}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
  21 ## @deftypefnx {Function File} {@var{x} =} gmres (@var{A}, @var{b}, @var{m}, @var{rtol}, @var{maxit}, @var{P})
  22 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} gmres (@dots{})
  23 ## Solve @code{A x = b} using the Preconditioned GMRES iterative method
  24 ## with restart, a.k.a. PGMRES(m).
  25 ##
  26 ## @itemize @minus
  27 ## @item @var{rtol} is the relative tolerance,
  28 ## if not given or set to [] the default value 1e-6 is used.
  29 ##
  30 ## @item @var{maxit} is the maximum number of outer iterations,
  31 ## if not given or set to [] the default value
  32 ## @code{min (10, numel (b) / restart)} is used.
  33 ##
  34 ## @item @var{x0} is the initial guess,
  35 ## if not given or set to [] the default value @code{zeros(size (b))} is used.
  36 ##
  37 ## @item @var{m} is the restart parameter,
  38 ## if not given or set to [] the default value @code{numel (b)} is used.
  39 ## @end itemize
  40 ##
  41 ## Argument @var{A} can be passed as a matrix, function handle, or
  42 ## inline function @code{f} such that @code{f(x) = A*x}.
  43 ##
  44 ## The preconditioner @var{P} is given as @code{P = M1 * M2}.
  45 ## Both @var{M1} and @var{M2} can be passed as a matrix, function handle, or
  46 ## inline function @code{g} such that @code{g(x) = M1\x} or @code{g(x) = M2\x}.
  47 ##
  48 ## Besides the vector @var{x}, additional outputs are:
  49 ##
  50 ## @itemize @minus
  51 ## @item @var{flag} indicates the exit status:
  52 ## @table @asis
  53 ## @item 0 : iteration converged to within the specified tolerance
  54 ##
  55 ## @item 1 : maximum number of iterations exceeded
  56 ##
  57 ## @item 2 : unused, but skipped for compatibility
  58 ##
  59 ## @item 3 : algorithm reached stagnation
  60 ## @end table
  61 ##
  62 ## @item @var{relres} is the final value of the relative residual.
  63 ##
  64 ## @item @var{iter} is a vector containing the number of outer iterations and
  65 ## total iterations performed.
  66 ##
  67 ## @item @var{resvec} is a vector containing the relative residual at each
  68 ## iteration.
  69 ## @end itemize
  70 ##
  71 ## @seealso{bicg, bicgstab, cgs, pcg}
  72 ## @end deftypefn
  73
  74 function [x, flag, presn, it, resids] = gmres (A, b, restart, rtol, maxit, M1, M2, x0)
  75
  76   if (nargin < 2 || nargin > 8)
  77     print_usage ();
  78   endif
  79
  80   if (ischar (A))
  81     Ax = str2func (A);
  82   elseif (ismatrix (A))
  83     Ax = @(x) A*x;
  84   elseif (isa (A, "function_handle"))
  85     Ax = A;
  86   else
  87     error ("gmres: A must be a function or matrix");
  88   endif
  89
  90   if (nargin < 3 || isempty (restart))
  91     restart = rows (b);
  92   endif
  93
  94   if (nargin < 4 || isempty (rtol))
  95     rtol = 1e-6;
  96   endif
  97
  98   if (nargin < 5 || isempty (maxit))
  99     maxit = min (rows (b)/restart, 10);
 100   endif
 101
 102   if (nargin < 6 || isempty (M1))
 103     M1m1x = @(x) x;
 104   elseif (ischar (M1))
 105     M1m1x = str2func (M1);
 106   elseif (ismatrix (M1))
 107     M1m1x = @(x) M1 \ x;
 108   elseif (isa (M1, "function_handle"))
 109     M1m1x = M1;
 110   else
 111     error ("gmres: preconditioner M1 must be a function or matrix");
 112   endif
 113
 114   if (nargin < 7 || isempty (M2))
 115     M2m1x = @(x) x;
 116   elseif (ischar (M2))
 117     M2m1x = str2func (M2);
 118   elseif (ismatrix (M2))
 119     M2m1x = @(x) M2 \ x;
 120   elseif (isa (M2, "function_handle"))
 121     M2m1x = M2;
 122   else
 123     error ("gmres: preconditioner M2 must be a function or matrix");
 124   endif
 125
 126   Pm1x = @(x) M2m1x (M1m1x (x));
 127
 128   if (nargin < 8 || isempty (x0))
 129     x0 = zeros (size (b));
 130   endif
 131
 132   x_old = x0;
 133   x = x_old;
 134   prec_res = Pm1x (b - Ax (x_old));
 135   presn = norm (prec_res, 2);
 136
 137   B = zeros (restart + 1, 1);
 138   V = zeros (rows (x), restart);
 139   H = zeros (restart + 1, restart);
 140
 141   ## begin loop
 142   iter = 1;
 143   restart_it  = restart + 1;
 144   resids      = zeros (maxit, 1);
 145   resids(1)   = presn;
 146   prec_b_norm = norm (Pm1x (b), 2);
 147   flag        = 1;
 148
 149   while (iter <= maxit * restart && presn > rtol * prec_b_norm)
 150
 151     ## restart
 152     if (restart_it > restart)
 153       restart_it = 1;
 154       x_old = x;
 155       prec_res = Pm1x (b - Ax (x_old));
 156       presn = norm (prec_res, 2);
 157       B(1) = presn;
 158       H(:) = 0;
 159       V(:, 1) = prec_res / presn;
 160     endif
 161
 162     ## basic iteration
 163     tmp = Pm1x (Ax (V(:, restart_it)));
 164     [V(:,restart_it+1), H(1:restart_it+1, restart_it)] = ...
 165         mgorth (tmp, V(:,1:restart_it));
 166
 167     Y = (H(1:restart_it+1, 1:restart_it) \ B (1:restart_it+1));
 168
 169     little_res = B(1:restart_it+1) - ...
 170         H(1:restart_it+1, 1:restart_it) * Y(1:restart_it);
 171
 172     presn = norm (little_res, 2);
 173
 174     x = x_old + V(:, 1:restart_it) * Y(1:restart_it);
 175
 176     resids(iter) = presn;
 177     if (norm (x - x_old, inf) <= eps)
 178       flag = 3;
 179       break
 180     endif
 181
 182     restart_it++ ;
 183     iter++;
 184   endwhile
 185
 186   if (presn > rtol * prec_b_norm)
 187     flag = 0;
 188   endif
 189
 190   resids = resids(1:iter-1);
 191   it = [ceil(iter / restart), rem(iter, restart)];
 192
 193 endfunction
 194
 195
 196 %!shared A, b, dim
 197 %! dim = 100;
 198 %!test
 199 %! A = spdiags ([-ones(dim,1) 2*ones(dim,1) ones(dim,1)], [-1:1], dim, dim);
 200 %! b = ones(dim, 1);
 201 %! x = gmres (A, b, 10, 1e-10, dim, @(x) x./diag(A), [],  b);
 202 %! assert(x, A\b, 1e-9*norm(x,inf));
 203 %!
 204 %!test
 205 %! x = gmres (A, b, dim, 1e-10, 1e4, @(x) diag(diag(A))\x, [],  b);
 206 %! assert(x, A\b, 1e-7*norm(x,inf));
 207 %!
 208 %!test
 209 %! 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);
 210 %! A = A'*A;
 211 %! b = rand (dim, 1);
 212 %! [x, resids] = gmres (@(x) A*x, b, dim, 1e-10, dim, @(x) x./diag(A), [],  []);
 213 %! assert(x, A\b, 1e-9*norm(x,inf))
 214 %! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) diag(diag(A))\x, [],  []);
 215 %! assert(x, A\b, 1e-9*norm(x,inf));
 216 %!test
 217 %! x =  gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) x./diag(A), [],  []);
 218 %! assert(x, A\b, 1e-7*norm(x,inf));