octave_packages/m/sparse/cgs.m

   1 ## Copyright (C) 2008-2012 Radek Salac
   2 ## Copyright (C) 2012 Carlo de Falco
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ##
  22 ## @deftypefn  {Function File} {@var{x} =} cgs (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
  23 ## @deftypefnx {Function File} {@var{x} =} cgs (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
  24 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} cgs (@var{A}, @var{b}, @dots{})
  25 ## Solve @code{A x = b}, where @var{A} is a square matrix, using the
  26 ## Conjugate Gradients Squared method.
  27 ##
  28 ## @itemize @minus
  29 ## @item @var{rtol} is the relative tolerance, if not given or set to []
  30 ## the default value 1e-6 is used.
  31 ##
  32 ## @item @var{maxit} the maximum number of outer iterations, if not
  33 ## given or set to [] the default value @code{min (20, numel (b))} is
  34 ## used.
  35 ##
  36 ## @item @var{x0} the initial guess, if not given or set to [] the
  37 ## default value @code{zeros (size (b))} is used.
  38 ## @end itemize
  39 ##
  40 ## @var{A} can be passed as a matrix or as a function handle or
  41 ## inline function @code{f} such that @code{f(x) = A*x}.
  42 ##
  43 ## The preconditioner @var{P} is given as @code{P = M1 * M2}.
  44 ## Both @var{M1} and @var{M2} can be passed as a matrix or as a function
  45 ## handle or inline function @code{g} such that @code{g(x) = M1 \ x} or
  46 ## @code{g(x) = M2 \ x}.
  47 ##
  48 ## If called with more than one output parameter
  49 ##
  50 ## @itemize @minus
  51 ## @item @var{flag} indicates the exit status:
  52 ## @itemize @minus
  53 ## @item 0: iteration converged to the within the chosen tolerance
  54 ##
  55 ## @item 1: the maximum number of iterations was reached before convergence
  56 ##
  57 ## @item 3: the algorithm reached stagnation
  58 ## @end itemize
  59 ## (the value 2 is unused but skipped for compatibility).
  60 ##
  61 ## @item @var{relres} is the final value of the relative residual.
  62 ##
  63 ## @item @var{iter} is the number of iterations performed.
  64 ##
  65 ## @item @var{resvec} is a vector containing the relative residual at
  66 ## each iteration.
  67 ## @end itemize
  68 ##
  69 ## @seealso{pcg, bicgstab, bicg, gmres}
  70 ## @end deftypefn
  71
  72 function [x, flag, relres, iter, resvec] = cgs (A, b, tol, maxit, M1, M2, x0)
  73
  74   if (nargin >= 2 && nargin <= 7 && isvector (full (b)))
  75
  76     if (ischar (A))
  77       A = str2func (A);
  78     elseif (ismatrix (A))
  79       Ax = @(x) A * x;
  80     elseif (isa (A, "function_handle"))
  81       Ax = @(x) feval (A, x);
  82     else
  83       error (["cgs: first argument is expected to "...
  84               "be a function or a square matrix"]);
  85     endif
  86
  87     if (nargin < 3 || isempty (tol))
  88       tol = 1e-6;
  89     endif
  90
  91     if (nargin < 4 || isempty (maxit))
  92       maxit = min (rows (b), 20);
  93     endif
  94
  95     if (nargin < 5 || isempty (M1))
  96       M1m1x = @(x) x;
  97     elseif (ischar (M1))
  98       M1m1x = str2func (M1);
  99     elseif (ismatrix (M1))
 100       M1m1x = @(x) M1 \ x;
 101     elseif (isa (M1, "function_handle"))
 102       M1m1x = @(x) feval (M1, x);
 103     else
 104       error ("cgs: preconditioner is expected to be a function or matrix");
 105     endif
 106
 107     if (nargin < 6 || isempty (M2))
 108       M2m1x = @(x) x;
 109     elseif (ischar (M2))
 110       M2m1x = str2func (M2);
 111     elseif (ismatrix (M2))
 112       M2m1x = @(x) M2 \ x;
 113     elseif (isa (M2, "function_handle"))
 114       M2m1x = @(x) feval (M2, x);
 115     else
 116       error ("cgs: preconditioner is expected to be a function or matrix");
 117     endif
 118
 119     precon = @(x) M2m1x (M1m1x (x));
 120
 121     if (nargin < 7 || isempty (x0))
 122       x0 = zeros (size (b));
 123     endif
 124
 125
 126     x = x0;
 127
 128     res = b - Ax (x);
 129     norm_b = norm (b);
 130     ## Vector of the residual norms for each iteration.
 131     resvec = norm (res) / norm_b;
 132     ro = 0;
 133     ## Default behavior we don't reach tolerance tol within maxit iterations.
 134     flag = 1;
 135     for iter = 1:maxit
 136
 137       z = precon (res);
 138
 139       ## Cache.
 140       ro_old = ro;
 141       ro = res' * z;
 142       if (iter == 1)
 143         p = z;
 144       else
 145         beta = ro / ro_old;
 146         p = z + beta * p;
 147       endif
 148       ## Cache.
 149       q = Ax (p);
 150       alpha = ro / (p' * q);
 151       x = x + alpha * p;
 152
 153       res = res - alpha * q;
 154       relres = norm (res) / norm_b;
 155       resvec = [resvec; relres];
 156
 157       if (relres <= tol)
 158         ## We reach tolerance tol within maxit iterations.
 159         flag = 0;
 160         break
 161       elseif (resvec (end) == resvec (end - 1))
 162         ## The method stagnates.
 163         flag = 3;
 164         break
 165       endif
 166     endfor
 167
 168     if (nargout < 1)
 169       if (flag == 0)
 170         printf ("cgs converged at iteration %i to a solution with relative residual %e\n",
 171                 iter, relres);
 172       elseif (flag == 3)
 173         printf (["cgs stopped at iteration %i without converging to the desired tolerance %e\n",
 174                  "because the method stagnated.\n",
 175                  "The iterate returned (number %i) has relative residual %e\n"],
 176                 iter, tol, iter, relres);
 177       else
 178         printf (["cgs stopped at iteration %i without converging to the desired tolerance %e\n",
 179                  "because the maximum number of iterations was reached.\n",
 180                  "The iterate returned (number %i) has relative residual %e\n"],
 181                 iter, tol, iter, relres);
 182       endif
 183     endif
 184
 185   else
 186     print_usage ();
 187   endif
 188
 189 endfunction
 190
 191
 192
 193 %!demo
 194 %! % Solve system of A*x=b
 195 %! A=[5 -1 3;-1 2 -2;3 -2 3]
 196 %! b=[7;-1;4]
 197 %! [a,b,c,d,e]=cgs(A,b)
 198
 199 %!shared A, b, n, M
 200 %!
 201 %!test
 202 %! n = 100;
 203 %! A = spdiags ([-ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
 204 %! b = sum (A, 2);
 205 %! tol = 1e-8;
 206 %! maxit = 1000;
 207 %! M = 4*eye (n);
 208 %! [x, flag, relres, iter, resvec] = cgs (A, b, tol, maxit, M);
 209 %! assert (x, ones (size (b)), 1e-7);
 210 %!
 211 %!test
 212 %! tol = 1e-8;
 213 %! maxit = 15;
 214 %!
 215 %! [x, flag, relres, iter, resvec] = cgs (@(x) A * x, b, tol, maxit, M);
 216 %! assert (x, ones (size (b)), 1e-7);
 217
 218 %!test
 219 %! n = 100;
 220 %! tol = 1e-8;
 221 %! a = sprand (n, n, .1);
 222 %! A = a'*a + 100 * eye (n);
 223 %! b = sum (A, 2);
 224 %! [x, flag, relres, iter, resvec] = cgs (A, b, tol, [], diag (diag (A)));
 225 %! assert (x, ones (size (b)), 1e-7);