function [ x, ierr ] = broyden ( x, f, atol, rtol, maxit, maxdim )

%*****************************************************************************80
%
%% broyden() uses Broyden's method for a nonlinear system of equations.
%
%  Discussion:
%
%    Given an initial estimate X0, a point X* is sought for which 
%      ||F(X*)|| <= ATOL + RTOL * ||F(X0)||
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    16 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 X(N): the initial solution estimate.
%
%    function fx = F ( x ): MATLAB code that evaluates the function.
%
%    real ATOL, RTOL: error tolerances.
%
%    integer MAXIT: maximum number of nonlinear iterations
%
%    integer MAXDIM: maximum number of Broyden iterations before restart.,
%
%  Output:
%
%    real X(N): the estimated solution.
%
%    integer IERR: termination flag.
%    0: success
%    1: the termination criterion was not satisfied after MAXIT iterations,
%       or the ratio of successive nonlinear residuals exceeded 1. 
%
  x = x(:);
  n = size ( x, 1 ); 
%
%  Evaluate f at the initial iterate.
%
  fc = feval ( f, x );
  fnrm = norm ( fc ) / sqrt ( n );
%
%  Compute the stopping tolerance.
%
  stop_tol = atol + rtol * fnrm;
%
%  Initialize the step array.
%
  stp = zeros ( n, maxdim );
  stp(:,1) = - fc;
%
%  Initialize the step norm squared array.
%
  stp_nrm = zeros(maxdim,1);
  stp_nrm(1) = stp(:,1)' * stp(:,1);
%
%  Iteration:
%
  nbroy = 0;
  itc = 0; 

  while ( itc < maxit )

    nbroy = nbroy + 1;
    fnrmo = fnrm;
    itc = itc + 1;
%
%  Next iterate.
%
    x = x + stp(:,nbroy);
    fc = feval ( f, x );
    fnrm = norm ( fc ) / sqrt ( n );
%
%  Accept?
%
    if ( fnrm <= stop_tol )
      ierr = 0;
      return;
    end
%
%  Reject?
%
    if ( fnrmo <= fnrm )
      ierr = 1;
      fprintf ( 1, 'broyden: Residual ratio increase = %g\n', fnrm / fnrmo )
      return;
    end
%
%  Continue. 
%
    if ( nbroy < maxdim + 1 )

      z = - fc;
      if ( 1 < nbroy )
        for kbr = 1 : nbroy - 1
          z = z + stp(:,kbr+1) * ( ( stp(:,kbr)' * z ) / stp_nrm(kbr) );
        end
      end

      zz = stp(:,nbroy)' * z;
      zz = zz / stp_nrm(nbroy);
      stp(:,nbroy+1) = z / ( 1.0 - zz );
      stp_nrm(nbroy+1) = stp(:,nbroy+1)' * stp(:,nbroy+1);

    else
%
%  If out of room, restart.
%
      stp(:,1) = - fc;
      stp_nrm(1) = stp(:,1)' * stp(:,1);
      nbroy = 0;

    end

  end
%
%  Reached MAXIT iterations.  Did we make it?
%
  if ( stop_tol < fnrm )
    ierr = 1;
    disp ( 'broyden: Requested accuracy was not achieved.' )
  else
    ierr = 0;
  end

  return
end

