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

annotate scripts/sparse/bicgstab.m @ 13061:addfc0ae69c0

Make bicgstab interface more compatible * bicgstab.m: Add the possibility to pass a function handle for the coefficient matrix. Also add more tests.

author	Carlo de Falco <kingcrimson@tiscali.it>
date	Sat, 03 Sep 2011 20:06:30 +0200
parents	c792872f8942
children	a0d854f079d2

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
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	29 ## @item @var{rtol} is the relative tolerance, if not given or set to [] the default value 1e-6 is used.
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	30 ## @item @var{maxit} the maximum number of outer iterations, if not given or set to [] the default value @code{min (20, numel (b))} is used.
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	31 ## @item @var{x0} the initial guess, if not given or set to [] the default value @code{zeros (size (b))} is used.
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	32 ## @end itemize
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	33 ##
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	34 ## @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	35 ## 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	36 ##
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	37 ## 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	38 ## 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	39 ## 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	40 ##
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	41 ## If called with more than one output parameter
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	42 ##
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	43 ## @itemize @minus
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	44 ## @item @var{flag} indicates the exit status:
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	45 ## @itemize @minus
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	46 ## @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	47 ## @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	48 ## @item 3: the algorithm reached stagnation
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	49 ## @end itemize
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	50 ## (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	51 ## @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	52 ## @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	53 ## @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	54 ## @end itemize
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	55 ##
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	56 ## @seealso{pcg,cgs,bigc,gmres}
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	57 ##
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	58 ## @end deftypefn
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	59
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	60 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	61 M1, M2, x0)
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	62
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	63 if ((nargin >= 2) && (nargin <= 7) && isvector (full (b)))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	64
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	65 if (ischar (A))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	66 A = str2func (A);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	67 elseif (ismatrix (A))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	68 Ax = @(x) A * x;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	69 elseif (isa (A, "function_handle"))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	70 Ax = @(x) feval (A, x);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	71 else
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	72 error (["bicgstab: first argument is expected " ...
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	73 "to be a function or a square matrix"]);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	74 endif
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	75
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	76 if ((nargin < 3) \|\| (isempty (tol)))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	77 tol = 1e-6;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	78 endif
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	79
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	80 if ((nargin < 4) \|\| (isempty (maxit)))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	81 maxit = min (rows (b), 20);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	82 endif
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	83
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	84 if ((nargin < 5) \|\| isempty (M1))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	85 M1m1x = @(x) x;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	86 elseif (ischar (M1))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	87 M1m1x = str2func (M1);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	88 elseif (ismatrix (M1))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	89 M1m1x = @(x) M1 \ x;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	90 elseif (isa (M1, "function_handle"))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	91 M1m1x = @(x) feval (M1, x);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	92 else
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	93 error (["bicgstab: preconditioner is " ...
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	94 "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	95 endif
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	96
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	97 if ((nargin < 6) \|\| isempty (M2))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	98 M2m1x = @(x) x;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	99 elseif (ischar (M2))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	100 M2m1x = str2func (M2);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	101 elseif (ismatrix (M2))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	102 M2m1x = @(x) M2 \ x;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	103 elseif (isa (M2, "function_handle"))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	104 M2m1x = @(x) feval (M2, x);
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	105 else
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	106 error (["bicgstab: preconditioner is "...
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	107 "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	108 endif
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	109
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	110 precon = @(x) M2m1x (M1m1x (x));
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	111
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	112 if ((nargin < 7) \|\| (isempty (x0)))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	113 x0 = zeros (size (b));
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	114 endif
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	115
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	116
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	117 ## specifies initial estimate x0
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	118 if (nargin < 7)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	119 x = zeros (rows (b), 1);
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	120 else
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	121 x = x0;
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	122 endif
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	123
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	124 norm_b = norm (b);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	125
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	126 res = b - Ax (x);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	127 rr = res;
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	128
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	129 ## Vector of the residual norms for each iteration.
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	130 resvec = norm(res) / norm_b;
8416 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 ## 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	133 flag = 1;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	134
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	135 for iter = 1:maxit
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	136 rho_1 = res' * rr;
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	137
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	138 if (iter == 1)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	139 p = res;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	140 else
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	141 beta = (rho_1 / rho_2) * (alpha / omega);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	142 p = res + beta * (p - omega * v);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	143 endif
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	144
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	145 phat = precon (p);
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	146
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	147 v = Ax (phat);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	148 alpha = rho_1 / (rr' * v);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	149 s = res - alpha * v;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	150
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	151 shat = precon (s);
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 t = Ax (shat);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	154 omega = (t' * s) / (t' * t);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	155 x = x + alpha * phat + omega * shat;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	156 res = s - omega * t;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	157 rho_2 = rho_1;
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	158
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	159 relres = norm (res) / norm_b;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	160 resvec = [resvec; relres];
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	161
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	162 if (relres <= tol)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	163 ## We reach tolerance tol within maxit iterations.
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	164 flag = 0;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	165 break;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	166 elseif (resvec(end) == resvec(end - 1))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	167 ## The method stagnates.
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	168 flag = 3;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	169 break;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	170 endif
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	171 endfor
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	172
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	173 if (nargout < 2)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	174 if (flag == 0)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	175 printf ("bicgstab converged at iteration %i ", iter);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	176 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	177 elseif (flag == 3)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	178 printf ("bicgstab stopped at iteration %i ", iter);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	179 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	180 printf ("because the method stagnated.\n");
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	181 printf ("The iterate returned (number %i) ", iter);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	182 printf ("has relative residual %e\n", relres);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	183 else
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	184 printf ("bicgstab stopped at iteration %i ", iter);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	185 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	186 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	187 printf ("The iterate returned (number %i) ", iter);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	188 printf ("has relative residual %e\n", relres);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	189 endif
9278 2b35cb600d50 improve bicgstab and cgs Radek Salac <salac.r@gmail.com> parents: 8920 diff changeset	190 endif
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	191
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	192 else
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	193 print_usage ();
9278 2b35cb600d50 improve bicgstab and cgs Radek Salac <salac.r@gmail.com> parents: 8920 diff changeset	194 endif
2b35cb600d50 improve bicgstab and cgs Radek Salac <salac.r@gmail.com> parents: 8920 diff changeset	195
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	196 endfunction
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	197
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	198 %!demo
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	199 %! % Solve system of A*x=b
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	200 %! A = [5 -1 3;-1 2 -2;3 -2 3]
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	201 %! b = [7;-1;4]
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	202 %! [x, flag, relres, iter, resvec] = bicgstab(A, b)
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	203
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	204 %!shared A, b, n, M1, M2
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	205 %!
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	206 %!test
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	207 %! n = 100;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	208 %! 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	209 %! b = sum (A, 2);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	210 %! tol = 1e-8;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	211 %! maxit = 15;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	212 %! 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	213 %! 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	214 %! [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	215 %! assert (x, ones (size (b)), 1e-7);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	216 %!
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	217 %!test
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 %!
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	221 %! function y = afun (x, a)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	222 %! y = a * x;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	223 %! endfunction
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 %! [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	226 %! tol, maxit, M1, M2);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	227 %! assert (x, ones (size (b)), 1e-7);
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 %!test
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	230 %! n = 100;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	231 %! tol = 1e-8;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	232 %! a = sprand (n, n, .1);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	233 %! A = a'a + 100 eye (n);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	234 %! b = sum (A, 2);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	235 %! [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	236 %! assert (x, ones (size (b)), 1e-7);