function [ x, error, it, flag ] = bicg_mine ( A, x, b, it_max, tol )

%*****************************************************************************80
%
%% bicg_mine() implements the biconjugate gradient method.
%
%  Discussion:
%
%    In order to avoid conflicts with the built-in bicg() function,
%    this function has been renamed "bicg_mine".
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 December 2022
%
%  Reference:
%
%    Richard Barrett, Michael Berry, Tony Chan, James Demmel,
%    June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo,
%    Charles Romine, Henk van der Vorst,
%    Templates for the Solution of Linear Systems:
%    Building Blocks for Iterative Methods,
%    SIAM, 1994,
%    ISBN: 0898714710,
%    LC: QA297.8.T45.
%
%  Input:
%
%    real A(*,*): the matrix;
%
%    real x(:): an initial solution estimate.
%
%    real b(:): the right hand side;
%
%    integer it_max: the maximum number of iterations.
%
%    real tol: an eror tolerance.
%
%  Output:
%
%    real x(:): the updated solution estimate.
%
%    real error: the norm of (b-A*x)/|b|.
%
%    integer it: the number of iterations used.
%
%    integer flag:
%    0: a solution was found with error less than the tolerance.
%    1: the maximum number of steps were take without convergence.
%   -1: the algorithm broke down.
%
  it = 0;
  flag = 0;

  bnrm = norm ( b );
  if ( bnrm == 0.0 )
    bnrm = 1.0; 
  end

  r = b - A * x;
  error = norm ( r ) / bnrm;
  if ( error < tol ) 
    return;
  end

  r_tld  = r;

  for it = 1 : it_max

    z = r;
    z_tld = r_tld;
    rho = ( z' * r_tld );
    if ( rho == 0.0 ),
      break
    end

    if ( it == 1 )
      p = z;
      p_tld = z_tld;
    else
      beta = rho / rho_old;
      p     = z     + beta * p;
      p_tld = z_tld + beta * p_tld;
    end
%
%  Compute the residual pair.
%
    q = A * p;
    q_tld = A' * p_tld;
    alpha = rho / ( p_tld' * q );
%
%  Update the approximation.
%
    x = x + alpha * p;
    r = r - alpha * q;
    r_tld = r_tld - alpha * q_tld;
%
%  Check for convergence.
%
    error = norm ( r ) / bnrm;
    if ( error <= tol )
      break
    end

    rho_old = rho;

  end
%
%  Set the flag.
%
  if ( error <= tol )
    flag =  0;
  elseif ( rho == 0.0 )
    flag = -1;
  else
    flag = 1;
  end

  return
end

