function [ x, error, total_iters ] = fdcgstab ( f0, f, xc, params, xinit )

%*****************************************************************************80
%
%% fdcgstab() implements a forward difference Bi-CGSTAB solver for nsola().
%
%  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 F0(N): the function value at the current point.
%
%    function fx = F ( x ): MATLAB code that evaluates the function.
%
%    real XC(N): the current point.
%
%    real PARAMS(2): iteration parameters.
%    params(1): relative residual reduction factor
%    params(2): max number of iterations
%
%    real XINIT(N): initial iterate.
%
%  Output:
%
%    real X(N): the estimated solution.
%
%    real ERROR(total_iters): residual norms for the  iteration.
%
%    integer TOTAL_ITERS: the number of iterations taken.
%
  b = - f0; 
  n = length ( b );
  errtol = params(1) * norm ( b );
  kmax = params(2);
  error = []; 
  rho = zeros ( kmax+1, 1 );
%
%  Use zero vector as initial iterate for Newton step unless
%  the calling routine has a better idea (useful for GMRES(m)).
%
  x = zeros ( n, 1 );
  r = b;
  if ( nargin == 5 )
    x = xinit;
    r = - dirder ( xc, x, f, f0 ) - f0;
  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 );
  error = [ error, zeta ];
%
%  Bi-CGSTAB iteration
%
  while ( ( errtol < zeta ) & ( k < kmax ) )
    k = k + 1;
    if ( omega == 0.0 )
      error('fdcgstab: breakdown, omega=0');
    end
    beta = ( rho(k+1) / rho(k) ) * ( alpha / omega );
    p = r + beta * ( p - omega * v );
    v = dirder ( xc, p, f, f0 );
    tau = hatr0' * v;
    if ( tau == 0.0 )
      error('fdcgstab: breakdown, tau=0');
    end 
    alpha = rho(k+1) / tau;
    s = r - alpha * v; 
    t = dirder ( xc, s, f, f0 );
    tau = t' * t;
    if ( tau == 0.0 )
      error('fdcgstab: 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;
    error = [ error, zeta ];
  end

  return
end

