function [sol, it_hist, ierr] = brsol(x,f,tol, parms)
% Broyden's Method solver, locally convergent
% solver for f(x) = 0
%
% C. T. Kelley, June 29, 1994
%
% This code comes with no guarantee or warranty of any kind.
%
% function [sol, it_hist, ierr] = brsol(x,f,tol,parms)
%
% inputs:
% initial iterate = x
% function = f
% tol = [atol, rtol] relative/absolute
% error tolerances for the nonlinear iteration
% parms = [maxit, maxdim, linprob]
% maxit = maxmium number of nonlinear iterations
% default = 40
% maxdim = maximum number of Broyden iterations
% before restart, so maxdim+3 vectors are
% stored (see text). By making the code a bit more
% subtle (putting z and F in the same place) one can
% reduce this overhead to maxdim+2 vectors.
% default = 40
% linprob = 0 for nonlinear problem
% = 1 for linear problem
% if linprob = 1 an increase in the residual does
% not result in termination
% default = 0
%
% output:
% sol = solution
% it_hist(maxit) = scaled l2 norms of nonlinear residuals
% ierr = 0 upon successful termination
% ierr = 1 if either after maxit iterations
% the termination criterion is not satsified
% or the ratio of successive nonlinear residuals
% exceeds 1. In this latter case, the iteration
% is terminated.
%
%
% internal parameter:
% debug = turns on/off iteration statistics display as
% the iteration progresses
%
% set the debug parameter, 1 turns display on, otherwise off
%
debug=1;
%
% initialize it_hist, ierr, and set the iteration parameters
%
ierr = 0; maxit=40; maxdim=39; linprob = 0; it_histx=zeros(maxit);
%
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; itc=0; nbroy=0;
%
% evaluate f at the initial iterate
% compute the stop tolerance
%
f0=feval(f,x);
fc=f0;
fnrm=norm(f0)/sqrt(n);
it_hist(itc+1)=fnrm;
fnrmo=1;
stop_tol=atol + rtol*fnrm;
outstat(itc+1, :) = [itc fnrm 0];
%
% 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
%
stp(:,1) = -fc;
stp_nrm(1)=stp(:,1)'*stp(:,1);
%
% main iteration loop
%
while(itc < maxit)
%
nbroy=nbroy+1;
%
% keep track of successive residual norms and
% the iteration counter (itc)
%
fnrmo=fnrm; itc=itc+1;
%
% compute the new point, test for termination before
% adding to iteration history
%
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==1
disp(outstat(itc+1,:))
end
%
% test for termination before computing the next w
%
if fnrm <= stop_tol
sol=x;
return;
end
%
% if residual norms increase, terminate, set error flag
%
if rat >= 1 & linprob == 0
ierr=1;
disp('increase in residual')
disp(outstat)
return;
end
%
% if there's room, compute the step and step norm and
% add to the iteration history
%
if nbroy < maxdim+1
z=-fc;
if nbroy > 1
for kbr = 1:nbroy-1
z=z+stp(:,kbr+1)*((stp(:,kbr)'*z)/stp_nrm(kbr));
end
end
%
% store the new step and step norm
%
zz=stp(:,nbroy)'*z;
zz=zz/stp_nrm(nbroy);
stp(:,nbroy+1)=z/(1-zz);
stp_nrm(nbroy+1)=stp(:,nbroy+1)'*stp(:,nbroy+1);
%
%
%
else
%
% out of room, time to restart
%
stp(:,1)=-fc;
stp_nrm(1)=stp(:,1)'*stp(:,1);
nbroy=0;
%
%
%
end
%
% end while
end
sol=x;
if debug==1
disp(outstat)
end
%
% on failure, set the error flag
%
if fnrm > stop_tol
ierr = 1;
end
