function [ sol, it_hist, ierr ] = broyden ( x, f, tol, parms )

%*****************************************************************************80
%
%% broyden() applies a locally convergent version of Broyden's Method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 June 1994
%
%  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 TOL(2): [ atol, rtol ] error tolerances.
%    atol: absolute error tolerance.
%    rtol: relative error tolerance.
%
%    integer PARMS(3): [ maxit, maxdim, linprob ]
%    integer MAXIT: maximum number of nonlinear iterations
%    integer MAXDIM: maximum number of Broyden iterations before restart.,
%    logical LINPROB: false for nonlinear problem, true for linear.
%    On a nonlinear problem, the iteration halts immediately if the
%    residual increases.
%
%  Output:
%
%    real SOL(N): the estimated solution.
%
%    integer IT_HIST(MAXIT): scaled l2 norms of nonlinear residuals.
%
%    integer IERR: termination flag.
%    0: success
%    1: the termination criterion was not satsified after MAXIT iterations,
%       or the ratio of successive nonlinear residuals exceeded 1. 
%
%  Local:
%
%    logical DEBUG: if TRUE, prints iteration statistics.
%
  debug = false;
%
%  Initialize.
%
  ierr = 0; 
  maxit = 40; 
  maxdim = 39; 
  linprob = false; 
  it_histx = zeros ( maxit, 1 );

  if ( nargin == 4 )
    maxit = parms(1); 
    maxdim = parms(2)-1; 
    linprob = parms(3);
  end

  rtol = tol(2); 
  atol = tol(1); 
  n = length ( x ); 
  fnrm = 1.0; 
  itc = 0; 
  nbroy = 0;
%
%  Evaluate f at the initial iterate.
%
  f0 = feval ( f, x );
  fc = f0;
  fnrm = norm ( f0 ) / sqrt ( n );
  it_hist(itc+1) = fnrm;
  fnrmo = 1.0;
  outstat(itc+1,:) = [ itc, fnrm, 0.0 ];
%
%  Compute the stop tolerance.
%
  stop_tol = atol + rtol * fnrm;
%
%  Initialize the iteration history storage matrices.
%
  stp = zeros(n,maxdim);
  stp_nrm = zeros(maxdim,1);
%
%  Set the initial step to -F, compute the step norm squared.
%
  stp(:,1) = - fc;
  stp_nrm(1) = stp(:,1)' * stp(:,1);
%
%  Iteration:
%
  while ( itc < maxit )

    nbroy = nbroy + 1;
    fnrmo = fnrm;
    itc = itc + 1;
%
%   Compute the new point.
%
    x = x + stp(:,nbroy);
    fc = feval ( f, x );
    fnrm = norm ( fc ) / sqrt ( n );
    it_hist(itc+1) = fnrm; 
    rat = fnrm / fnrmo;
    outstat(itc+1,:) = [ itc, fnrm, rat ];
    if ( debug )
      disp(outstat(itc+1,:))
    end
%
%  Stop if we satisfy the error tolerance.
%
    if ( fnrm <= stop_tol )
      sol = x;
      return;
    end
%
%  If the residual norm increased, and not a linear problem, terminate.
%
    if ( 1.0 <= rat && ~ linprob )
      ierr = 1;
      disp('broyden: increase in residual')
      disp(outstat)
      return;
    end
%
%  If room, compute step and step norm. 
%
    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

  sol = x;

  if ( debug )
    disp ( outstat )
  end
%
%  on failure, set the error flag
%
  if ( stop_tol < fnrm )
    ierr = 1;
    disp ( 'broyden: Requested accuracy was not achieved.' )
  end

  return
end

