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

annotate scripts/sparse/bicgstab.m @ 14868:5d3a684236b0

maint: Use Octave coding conventions for cuddling parentheses in scripts directory * lin2mu.m, loadaudio.m, wavread.m, accumarray.m, bicubic.m, celldisp.m, colon.m, cplxpair.m, dblquad.m, divergence.m, genvarname.m, gradient.m, int2str.m, interp1.m, interp1q.m, interp2.m, interpn.m, loadobj.m, nthargout.m, __isequal__.m, __splinen__.m, quadgk.m, quadl.m, quadv.m, rat.m, rot90.m, rotdim.m, saveobj.m, subsindex.m, triplequad.m, delaunay3.m, griddata.m, inpolygon.m, tsearchn.m, voronoi.m, get_first_help_sentence.m, which.m, gray2ind.m, pink.m, dlmwrite.m, strread.m, textread.m, textscan.m, housh.m, ishermitian.m, issymmetric.m, krylov.m, logm.m, null.m, rref.m, compare_versions.m, copyfile.m, dump_prefs.m, edit.m, fileparts.m, getappdata.m, isappdata.m, movefile.m, orderfields.m, parseparams.m, __xzip__.m, rmappdata.m, setappdata.m, swapbytes.m, unpack.m, ver.m, fminbnd.m, fminunc.m, fsolve.m, glpk.m, lsqnonneg.m, qp.m, sqp.m, configure_make.m, copy_files.m, describe.m, get_description.m, get_forge_pkg.m, install.m, installed_packages.m, is_architecture_dependent.m, load_package_dirs.m, print_package_description.m, rebuild.m, repackage.m, save_order.m, shell.m, allchild.m, ancestor.m, area.m, axes.m, axis.m, clabel.m, close.m, colorbar.m, comet.m, comet3.m, contour.m, cylinder.m, ezmesh.m, ezsurf.m, findobj.m, fplot.m, hist.m, isocolors.m, isonormals.m, isosurface.m, isprop.m, legend.m, mesh.m, meshz.m, pareto.m, pcolor.m, peaks.m, plot3.m, plotmatrix.m, plotyy.m, polar.m, print.m, __add_datasource__.m, __add_default_menu__.m, __axes_limits__.m, __bar__.m, __clabel__.m, __contour__.m, __errcomm__.m, __errplot__.m, __ezplot__.m, __file_filter__.m, __fltk_print__.m, __ghostscript__.m, __gnuplot_print__.m, __go_draw_axes__.m, __go_draw_figure__.m, __interp_cube__.m, __marching_cube__.m, __patch__.m, __pie__.m, __plt__.m, __print_parse_opts__.m, __quiver__.m, __scatter__.m, __stem__.m, __tight_eps_bbox__.m, __uigetdir_fltk__.m, __uigetfile_fltk__.m, __uiputfile_fltk__.m, quiver.m, quiver3.m, rectangle.m, refreshdata.m, ribbon.m, scatter.m, semilogy.m, shading.m, slice.m, subplot.m, surface.m, surfl.m, surfnorm.m, text.m, uigetfile.m, uiputfile.m, whitebg.m, deconv.m, mkpp.m, pchip.m, polyaffine.m, polyder.m, polygcd.m, polyout.m, polyval.m, ppint.m, ppjumps.m, ppval.m, residue.m, roots.m, spline.m, splinefit.m, addpref.m, getpref.m, setpref.m, ismember.m, setxor.m, arch_fit.m, arch_rnd.m, arch_test.m, autoreg_matrix.m, diffpara.m, fftconv.m, filter2.m, hanning.m, hurst.m, periodogram.m, triangle_sw.m, sinc.m, spectral_xdf.m, spencer.m, stft.m, synthesis.m, unwrap.m, yulewalker.m, bicgstab.m, gmres.m, pcg.m, pcr.m, __sprand_impl__.m, speye.m, spfun.m, sprandn.m, spstats.m, svds.m, treelayout.m, treeplot.m, bessel.m, factor.m, legendre.m, perms.m, primes.m, magic.m, toeplitz.m, corr.m, cov.m, mean.m, median.m, mode.m, qqplot.m, quantile.m, ranks.m, zscore.m, logistic_regression_likelihood.m, bartlett_test.m, chisquare_test_homogeneity.m, chisquare_test_independence.m, kolmogorov_smirnov_test.m, run_test.m, u_test.m, wilcoxon_test.m, z_test.m, z_test_2.m, bin2dec.m, dec2base.m, mat2str.m, strcat.m, strchr.m, strjust.m, strtok.m, substr.m, untabify.m, assert.m, demo.m, example.m, fail.m, speed.m, test.m, now.m: Use Octave coding conventions for cuddling parentheses in scripts directory.

author	Rik <octave@nomad.inbox5.com>
date	Tue, 17 Jul 2012 07:08:39 -0700
parents	72b8b39e12be
children	12005245b645

rev	line source
14138 72c96de7a403 maint: update copyright notices for 2012 John W. Eaton <jwe@octave.org> parents: 13929 diff changeset	1 ## Copyright (C) 2008-2012 Radek Salac
72c96de7a403 maint: update copyright notices for 2012 John W. Eaton <jwe@octave.org> parents: 13929 diff changeset	2 ## Copyright (C) 2012 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 --
13929 9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	21 ## @deftypefn {Function File} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	22 ## @deftypefnx {Function File} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
13929 9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	23 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicgstab (@var{A}, @var{b}, @dots{})
9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	24 ## Solve @code{A x = b} using the stabilizied Bi-conjugate gradient iterative
9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	25 ## method.
13070 a0d854f079d2 codesprint: Fix building of docs for new bicg functions Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 13061 diff changeset	26 ##
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	27 ## @itemize @minus
e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	28 ## @item @var{rtol} is the relative tolerance, if not given or set to
e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	29 ## [] the default value 1e-6 is used.
13929 9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	30 ##
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	31 ## @item @var{maxit} the maximum number of outer iterations, if not
e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	32 ## given or set to [] the default value @code{min (20, numel (b))} is
e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	33 ## used.
13929 9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	34 ##
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	35 ## @item @var{x0} the initial guess, if not given or set to [] the
e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	36 ## default value @code{zeros (size (b))} is used.
e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	37 ## @end itemize
13070 a0d854f079d2 codesprint: Fix building of docs for new bicg functions Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 13061 diff changeset	38 ##
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	39 ## @var{A} can be passed as a matrix or as a function handle or
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	40 ## 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	41 ##
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	42 ## The preconditioner @var{P} is given as @code{P = M1 * M2}.
13929 9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	43 ## Both @var{M1} and @var{M2} can be passed as a matrix or as a function
9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	44 ## handle or inline function @code{g} such that @code{g(x) = M1 \ x} or
9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	45 ## @code{g(x) = M2 \ x}.
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	46 ##
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	47 ## If called with more than one output parameter
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	48 ##
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	49 ## @itemize @minus
e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	50 ## @item @var{flag} indicates the exit status:
14853 72b8b39e12be doc: Periodic grammarcheck of documentation. Rik <octave@nomad.inbox5.com> parents: 14609 diff changeset	51 ##
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	52 ## @itemize @minus
e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	53 ## @item 0: iteration converged to the within the chosen tolerance
13929 9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	54 ##
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	55 ## @item 1: the maximum number of iterations was reached before convergence
13929 9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	56 ##
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	57 ## @item 3: the algorithm reached stagnation
e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	58 ## @end itemize
e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	59 ## (the value 2 is unused but skipped for compatibility).
13929 9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	60 ##
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	61 ## @item @var{relres} is the final value of the relative residual.
13929 9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	62 ##
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	63 ## @item @var{iter} is the number of iterations performed.
13929 9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	64 ##
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	65 ## @item @var{resvec} is a vector containing the relative residual at each iteration.
e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	66 ## @end itemize
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	67 ##
13929 9cae456085c2 Grammarcheck of documentation before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 13837 diff changeset	68 ## @seealso{bicg, cgs, gmres, pcg}
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	69 ##
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	70 ## @end deftypefn
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	71
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	72 function [x, flag, relres, iter, resvec] = bicgstab (A, b, tol, maxit,
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	73 M1, M2, x0)
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	74
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	75 if (nargin >= 2 && nargin <= 7 && isvector (full (b)))
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	76
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	77 if (ischar (A))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	78 A = str2func (A);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	79 elseif (ismatrix (A))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	80 Ax = @(x) A * x;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	81 elseif (isa (A, "function_handle"))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	82 Ax = @(x) feval (A, x);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	83 else
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	84 error (["bicgstab: first argument is expected " ...
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	85 "to be a function or a square matrix"]);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	86 endif
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	87
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	88 if (nargin < 3 \|\| isempty (tol))
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	89 tol = 1e-6;
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
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	92 if (nargin < 4 \|\| isempty (maxit))
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	93 maxit = min (rows (b), 20);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	94 endif
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	95
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	96 if (nargin < 5 \|\| isempty (M1))
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	97 M1m1x = @(x) x;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	98 elseif (ischar (M1))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	99 M1m1x = str2func (M1);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	100 elseif (ismatrix (M1))
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	101 M1m1x = @(x) M1 \ x;
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	102 elseif (isa (M1, "function_handle"))
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	103 M1m1x = @(x) feval (M1, x);
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	104 else
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	105 error (["bicgstab: preconditioner is " ...
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	106 "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	107 endif
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	108
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	109 if (nargin < 6 \|\| isempty (M2))
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	110 M2m1x = @(x) x;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	111 elseif (ischar (M2))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	112 M2m1x = str2func (M2);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	113 elseif (ismatrix (M2))
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	114 M2m1x = @(x) M2 \ x;
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	115 elseif (isa (M2, "function_handle"))
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	116 M2m1x = @(x) feval (M2, x);
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	117 else
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	118 error (["bicgstab: preconditioner is "...
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	119 "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	120 endif
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	121
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	122 precon = @(x) M2m1x (M1m1x (x));
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	123
13091 e5aaba072d2b maint: style fixes in iterative linear solvers Carlo de Falco <kingcrimson@tiscali.it> parents: 13070 diff changeset	124 if (nargin < 7 \|\| isempty (x0))
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	125 x0 = zeros (size (b));
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	126 endif
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	127
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	128 ## specifies initial estimate x0
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	129 if (nargin < 7)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	130 x = zeros (rows (b), 1);
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	131 else
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	132 x = x0;
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	133 endif
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	134
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	135 norm_b = norm (b);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	136
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	137 res = b - Ax (x);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	138 rr = res;
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 ## Vector of the residual norms for each iteration.
14868 5d3a684236b0 maint: Use Octave coding conventions for cuddling parentheses in scripts directory Rik <octave@nomad.inbox5.com> parents: 14853 diff changeset	141 resvec = norm (res) / norm_b;
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	142
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	143 ## 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	144 flag = 1;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	145
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	146 for iter = 1:maxit
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	147 rho_1 = res' * rr;
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	148
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	149 if (iter == 1)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	150 p = res;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	151 else
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	152 beta = (rho_1 / rho_2) * (alpha / omega);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	153 p = res + beta * (p - omega * v);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	154 endif
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	155
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	156 phat = precon (p);
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	157
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	158 v = Ax (phat);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	159 alpha = rho_1 / (rr' * v);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	160 s = res - alpha * v;
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 shat = precon (s);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	163
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	164 t = Ax (shat);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	165 omega = (t' * s) / (t' * t);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	166 x = x + alpha * phat + omega * shat;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	167 res = s - omega * t;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	168 rho_2 = rho_1;
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	169
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	170 relres = norm (res) / norm_b;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	171 resvec = [resvec; relres];
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	172
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	173 if (relres <= tol)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	174 ## We reach tolerance tol within maxit iterations.
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	175 flag = 0;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	176 break;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	177 elseif (resvec(end) == resvec(end - 1))
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	178 ## The method stagnates.
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	179 flag = 3;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	180 break;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	181 endif
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	182 endfor
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	183
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	184 if (nargout < 2)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	185 if (flag == 0)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	186 printf ("bicgstab converged at iteration %i ", iter);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	187 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	188 elseif (flag == 3)
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	189 printf ("bicgstab stopped at iteration %i ", iter);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	190 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	191 printf ("because the method stagnated.\n");
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	192 printf ("The iterate returned (number %i) ", iter);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	193 printf ("has relative residual %e\n", relres);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	194 else
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	195 printf ("bicgstab stopped at iteration %i ", iter);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	196 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	197 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	198 printf ("The iterate returned (number %i) ", iter);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	199 printf ("has relative residual %e\n", relres);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	200 endif
9278 2b35cb600d50 improve bicgstab and cgs Radek Salac <salac.r@gmail.com> parents: 8920 diff changeset	201 endif
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	202
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	203 else
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	204 print_usage ();
9278 2b35cb600d50 improve bicgstab and cgs Radek Salac <salac.r@gmail.com> parents: 8920 diff changeset	205 endif
2b35cb600d50 improve bicgstab and cgs Radek Salac <salac.r@gmail.com> parents: 8920 diff changeset	206
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	207 endfunction
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	208
14237 11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	209
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	210 %!demo
0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	211 %! % Solve system of A*x=b
14237 11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	212 %! A = [5 -1 3;-1 2 -2;3 -2 3];
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	213 %! b = [7;-1;4];
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	214 %! [x, flag, relres, iter, resvec] = bicgstab (A, b)
8416 0a56e5c21c29 Add the bicgstab function Radek Salac <salac.r@gmail.com> parents: diff changeset	215
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	216 %!shared A, b, n, M1, M2
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	217 %!
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	218 %!test
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	219 %! n = 100;
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	220 %! A = spdiags ([-2ones(n,1) 4ones(n,1) -ones(n,1)], -1:1, n, n);
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	221 %! b = sum (A, 2);
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	222 %! tol = 1e-8;
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	223 %! maxit = 15;
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	224 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n);
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	225 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	226 %! [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	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
13837 2c80bbd87f5d don't define functions in test and demo blocks John W. Eaton <jwe@octave.org> parents: 13141 diff changeset	230 %!function y = afun (x, a)
2c80bbd87f5d don't define functions in test and demo blocks John W. Eaton <jwe@octave.org> parents: 13141 diff changeset	231 %! y = a * x;
2c80bbd87f5d don't define functions in test and demo blocks John W. Eaton <jwe@octave.org> parents: 13141 diff changeset	232 %!endfunction
2c80bbd87f5d don't define functions in test and demo blocks John W. Eaton <jwe@octave.org> parents: 13141 diff changeset	233 %!
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	234 %! tol = 1e-8;
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	235 %! maxit = 15;
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	236 %!
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	237 %! [x, flag, relres, iter, resvec] = bicgstab (@(x) afun (x, A), b,
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	238 %! tol, maxit, M1, M2);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	239 %! assert (x, ones (size (b)), 1e-7);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	240
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	241 %!test
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	242 %! n = 100;
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	243 %! tol = 1e-8;
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	244 %! a = sprand (n, n, .1);
addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	245 %! A = a'a + 100 eye (n);
13141 e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 13091 diff changeset	246 %! b = sum (A, 2);
13061 addfc0ae69c0 Make bicgstab interface more compatible Carlo de Falco <kingcrimson@tiscali.it> parents: 11587 diff changeset	247 %! [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	248 %! assert (x, ones (size (b)), 1e-7);
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14237 diff changeset	249