function [ x, err, total_iters ] = bicgstab ( x0, b, atv, params )

%*****************************************************************************80
%
%% bicgstab() is a bi-conjugate gradient stabilized solver for linear systems.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 April 2021
%
%  Author:
%
%    Tim Kelley
%
%  Reference:
%
%    Tim Kelley,
%    Iterative Methods for Linear and Nonlinear Equations,
%    SIAM, 2004,
%    ISBN: 0898713528,
%    LC: QA297.8.K45.
%    
%  Input:
%
%    real x0(n): the initial solution estimate.
%
%    real b(n): the right hand side.
%
%    function ax = atv(x):, a matrix-vector product routine computing A*x.
%
%    real params(2): iteration parameters.
%    params(1) = relative residual reduction factor.
%    params(2) = kmax = maximum number of iterations.
%
%  Output:
%
%    real x(n): the estimated solution.
%
%    real err(*): vector of iteration residual norms.
%
%    integer total_iters: number of iterations.
%
  x = x0;
  err = [];
  total_iters = 0;
%
%  Initialization.
%
  n = length ( b );
  errtol = params(1) * norm ( b );
  kmax = params(2);
  rho = zeros(kmax+1,1);

  if ( norm ( x ) ~= 0.0 )
    r = b - feval ( atv, x );
  else
    r = b;
  end

  hatr0 = r;
  k = 0;
  rho(1) = 1.0;
  alpha = 1.0;
  omega = 1.0;
  v = zeros ( n, 1 );
  p = zeros ( n, 1 );
  rho(2) = hatr0' * r;
  zeta = norm ( r );
  err = [ err, zeta ];
%
%  Bi-CGSTAB iteration
%
  while ( ( errtol < zeta ) && ( k < kmax ) )

    k = k + 1;

    if ( omega == 0.0 )
      error ( 'bicgstab: breakdown, omega=0' );
    end

    beta = ( rho(k+1) / rho(k) ) * ( alpha / omega );
    p = r + beta * ( p - omega * v );
    v = feval ( atv, p );
    tau = hatr0' * v;

    if ( tau == 0.0 )
      error ( 'bicgstab: breakdown, tau=0' );
    end 

    alpha = rho(k+1) / tau;
    s = r - alpha * v; 
    t = feval ( atv, s );
    tau = t' * t;

    if ( tau == 0.0 )
      error ( 'bicgstab: breakdown, |t|=0' );
    end

    omega = t' * s / tau; 
    rho(k+2) = - omega * ( hatr0' * t);
    x = x + alpha * p + omega * s;
    r = s - omega * t;
    zeta = norm ( r );
    total_iters = k;
    err = [ err, zeta ];

  end

  return
end

