octave-lyh: scripts/sparse/bicgstab.m annotate

annotate scripts/sparse/bicgstab.m @ 13070:a0d854f079d2

codesprint: Fix building of docs for new bicg functions

author	Jordi Gutiérrez Hermoso <jordigh@gmail.com>
date	Sat, 03 Sep 2011 14:37:38 -0500
parents	addfc0ae69c0
children	e5aaba072d2b

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