octave_packages/m/sparse/bicg.m

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