--- 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)
+