Mercurial > hg > octave-nkf
changeset 8416:0a56e5c21c29
Add the bicgstab function
author | Radek Salac <salac.r@gmail.com> |
---|---|
date | Tue, 23 Dec 2008 01:08:59 +0100 |
parents | fa12c67683d3 |
children | 654bcfb937bf |
files | scripts/ChangeLog scripts/sparse/Makefile.in scripts/sparse/bicgstab.m |
diffstat | 3 files changed, 185 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2008-11-21 Radek Salac <salac.r@gmail.com> + + * sparse/bicgstab.m: New function. + * sparse/Makefile.in (SOURCES): Add it here. + 2008-12-18 Daniel J Sebald <daniel.sebald@ieee.org> * time/datevec.m (__date_vfmt2sfmt__): New helper function.
--- a/scripts/sparse/Makefile.in +++ b/scripts/sparse/Makefile.in @@ -32,7 +32,7 @@ INSTALL_PROGRAM = @INSTALL_PROGRAM@ INSTALL_DATA = @INSTALL_DATA@ -SOURCES = cgs.m colperm.m etreeplot.m gplot.m nonzeros.m normest.m \ +SOURCES = bicgstab.m cgs.m colperm.m etreeplot.m gplot.m nonzeros.m normest.m \ pcg.m pcr.m spalloc.m spaugment.m spconvert.m spdiags.m speye.m \ spfun.m sphcat.m spones.m sprand.m sprandn.m sprandsym.m spstats.m \ spvcat.m spy.m treelayout.m treeplot.m
new file mode 100644 --- /dev/null +++ b/scripts/sparse/bicgstab.m @@ -0,0 +1,179 @@ +## Copyright (C) 2008 Radek Salac +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {} bicgstab (@var{A}, @var{b}) +## @deftypefnx {Function File} {} bicgstab (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0}) +## This procedure attempts to solve a system of linear equations A*x = b for x. +## The @var{A} must be square, symmetric and positive definite real matrix N*N. +## The @var{b} must be a one column vector with a length of N. +## The @var{tol} specifies the tolerance of the method, the default value is 1e-6. +## The @var{maxit} specifies the maximum number of iterations, the default value is min(20,N). +## The @var{M1} specifies a preconditioner, can also be a function handler which returns M\X. +## The @var{M2} combined with @var{M1} defines preconditioner as preconditioner=M1*M2. +## The @var{x0} is the initial guess, the default value is zeros(N,1). +## +## The value @var{x} is a computed result of this procedure. +## The value @var{flag} can be 0 when we reach tolerance in @var{maxit} iterations, 1 when +## we don't reach tolerance in @var{maxit} iterations and 3 when the procedure stagnates. +## The value @var{relres} is a relative residual - norm(b-A*x)/norm(b). +## The value @var{iter} is an iteration number in which x was computed. +## The value @var{resvec} is a vector of @var{relres} for each iteration. +## +## @end deftypefn + +function [x, flag, relres, iter, resvec] = bicgstab (A, b, tol, maxit, M1, M2, x0) + + if (nargin < 2 || nargin > 7 || nargout > 5) + print_usage (); + elseif (!isnumeric (A) || rows (A) != columns (A)) + error ("bicgstab: the first argument must be a n-by-n matrix"); + elseif (!isvector (b)) + error ("bicgstab: b must be a vector"); + elseif (!any (b)) + error ("bicgstab: b shuldn't be a vector of zeros"); + elseif (rows (A) != rows (b)) + error ("bicgstab: the first and second argument must have the same number of rows"); + elseif (nargin > 2 && !isscalar (tol)) + error ("bicgstab: tol must be a scalar"); + elseif (nargin > 3 && !isscalar (maxit)) + error ("bicgstab: maxit must be a scalar"); + elseif (nargin > 4 && ismatrix (M1) && (rows (M1) != rows (A) || columns (M1) != columns (A))) + error ("bicgstab: M1 must have the same number of rows and columns as A"); + elseif (nargin > 5 && (!ismatrix (M2) || rows (M2) != rows (A) || columns (M2) != columns (A))) + error ("bicgstab: M2 must have the same number of rows and columns as A"); + elseif (nargin > 6 && !isvector (x0)) + error ("bicgstab: x0 must be a vector"); + elseif (nargin > 6 && rows (x0) != rows (b)) + error ("bicgstab: xO must have the same number of rows as b"); + endif + + ## default toleration + if (nargin < 3) + tol = 1e-6; + endif + + ## default maximum number of iteration + if (nargin < 4) + maxit = min (rows (b), 20); + endif + + ## left preconditioner + if (nargin == 5) + precon = M1; + elseif (nargin > 5) + if (isparse(M1) && issparse(M2)) + precon = @(x) M1 * (M2 * x); + else + precon = M1 * M2; + endif + endif + + if (nargin > 4 && isnumeric (precon)) + ## precon can by also function + if (det (precon) != 0) + ## we can compute inverse preconditioner and use quicker algorithm + precon=inv (precon); + else + error ("bicgstab: preconditioner is ill conditioned"); + endif + + if (isinf (cond (precon))); + ## we must make test if preconditioner isn't ill conditioned + error ("bicgstab: preconditioner is ill conditioned"); + endif + endif + + ## specifies initial estimate x0 + if (nargin < 7) + x = zeros (rows (b), 1); + else + x = x0; + endif + + norm_b = norm (b); + + res = b - A*x; + rr = res; + + ## vector of the residual norms for each iteration + resvec = [norm(res)]; + + ## default behaviour we don't reach tolerance tol within maxit iterations + flag = 1; + + for iter = 1:maxit + rho_1 = res' * rr; + + if (iter == 1) + p = res; + else + beta = (rho_1 / rho_2) * (alpha / omega); + p = res + beta * (p - omega * v); + endif + + if (nargin > 4 && isnumeric (precon)) + phat = precon * p; + elseif (nargin > 4) + ## our preconditioner is a function + phat = feval (precon, p); + else + phat = p; + endif + + v = A * phat; + alpha = rho_1 / (rr' * v); + s = res - alpha * v; + + if (nargin > 4 && isnumeric (precon)) + shat = precon * s; + elseif (nargin > 4) + ## our preconditioner is a function + shat = feval (precon, s); + else + shat = s; + endif + + t = A * shat; + omega = (t' * s) / (t' * t); + x = x + alpha * phat + omega * shat; + res = s - omega * t; + rho_2 = rho_1; + + relres = norm (res) / norm_b; + resvec = [resvec; relres]; + + if (relres <= tol) + ## we reach tolerance tol within maxit iterations + flag = 0; + break; + elseif ( resvec (end) == resvec (end - 1)) + ## the method stagnates + flag = 3; + break; + endif + endfor + +endfunction + +%!demo +%! % Solve system of A*x=b +%! A = [5 -1 3;-1 2 -2;3 -2 3] +%! b = [7;-1;4] +%! [x, flag, relres, iter, resvec] = bicgstab(A, b) +