function [sol, it_hist, ierr] = brsola(x,f,tol, parms)% Broyden's Method solver, globally convergent% solver for f(x) = 0, Armijo rule, one vector storage%% C. T. Kelley, June 29, 1994%% This code comes with no guarantee or warranty of any kind.%% function [sol, it_hist, ierr] = brsola(x,f,tol,parms)%% inputs:%        initial iterate = x%        function = f%        tol = [atol, rtol] relative/absolute%            error tolerances for the nonlinear iteration%        parms = [maxit, maxdim]%            maxit = maxmium number of nonlinear iterations%                default = 40%            maxdim = maximum number of Broyden iterations%                before restart, so maxdim-1 vectors are %                stored%                default = 40%% output:%        sol = solution%        it_hist(maxit,3) = scaled l2 norms of nonlinear residuals%            for the iteration, number function evaluations,%            and number of steplength reductions%        ierr = 0 upon successful termination%        ierr = 1 if after maxit iterations%             the termination criterion is not satsified.%        ierr = 2 failure in the line search. The iteration%             is terminated if too many steplength reductions%             are taken.%%% internal parameter:%       debug = turns on/off iteration statistics display as%               the iteration progresses%%       alpha = 1.d-4, parameter to measure sufficient decrease%%       maxarm = 10, maximum number of steplength reductions before%                    failure is reported%% 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;  it_histx=zeros(maxit,3);maxarm=10;%if nargin == 4    maxit=parms(1); maxdim=parms(2)-1; endrtol=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;it_histx(itc+1,1)=fnrm; it_histx(itc+1,2)=0; it_histx(itc+1,3)=0;fnrmo=1;stop_tol=atol + rtol*fnrm;outstat(itc+1, :) = [itc fnrm 0 0];%% terminate on entry?%if fnrm < stop_tol    sol=x;    returnend%% initialize the iteration history storage matrices%stp=zeros(n,maxdim);stp_nrm=zeros(maxdim,1);lam_rec=ones(maxdim,1);%% Set the initial step to -F, compute the step norm%lambda=1;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%    xold=x; lambda=1; iarm=0; lrat=.5; alpha=1.d-4;    x = x + stp(:,nbroy);    fc=feval(f,x);    fnrm=norm(fc)/sqrt(n);    ff0=fnrmo*fnrmo; ffc=fnrm*fnrm; lamc=lambda;%%%   Line search, we assume that the Broyden direction is an%   ineact Newton direction. If the line search fails to%   find sufficient decrease after maxarm steplength reductions%   brsola returns with failure. %%   Three-point parabolic line search%    while fnrm >= (1 - lambda*alpha)*fnrmo & iarm < maxarm%       lambda=lambda*lrat;        if iarm==0            lambda=lambda*lrat;        else            lambda=parab3p(lamc, lamm, ff0, ffc, ffm);        end        lamm=lamc; ffm=ffc; lamc=lambda;        x = xold + lambda*stp(:,nbroy);        fc=feval(f,x);        fnrm=norm(fc)/sqrt(n);        ffc=fnrm*fnrm;        iarm=iarm+1;    end%%   set error flag and return on failure of the line search%    if iarm == maxarm        disp('Line search failure in brsola ')        ierr=2;        it_hist=it_histx(1:itc+1,:);        sol=xold;        return;    end%%   How many function evaluations did this iteration require?%    it_histx(itc+1,1)=fnrm;    it_histx(itc+1,2)=it_histx(itc,2)+iarm+1;    if(itc == 1) it_histx(itc+1,2) = it_histx(itc+1,2)+1; end;    it_histx(itc+1,3)=iarm;%%   terminate?%    if fnrm < stop_tol        sol=x;        rat=fnrm/fnrmo;        outstat(itc+1, :) = [itc fnrm iarm rat];        it_hist=it_histx(1:itc+1,:);%        it_hist(itc+1)=fnrm;        if debug==1            disp(outstat(itc+1,:))        end        return    end%%%   modify the step and step norm if needed to reflect the line %   search%    lam_rec(nbroy)=lambda;    if lambda ~= 1         stp(:,nbroy)=lambda*stp(:,nbroy);         stp_nrm(nbroy)=lambda*lambda*stp_nrm(nbroy);    end%%%    it_hist(itc+1)=fnrm;     rat=fnrm/fnrmo;    outstat(itc+1, :) = [itc fnrm iarm rat];        if debug==1            disp(outstat(itc+1,:))        end%%%   if there's room, compute the next search direction and step norm and%   add to the iteration history %    if nbroy < maxdim+1        z=-fc;        if nbroy > 1            for kbr = 1:nbroy-1                 ztmp=stp(:,kbr+1)/lam_rec(kbr+1);                 ztmp=ztmp+(1 - 1/lam_rec(kbr))*stp(:,kbr);                 ztmp=ztmp*lam_rec(kbr);                 z=z+ztmp*((stp(:,kbr)'*z)/stp_nrm(kbr));            end        end%%       store the new search direction and its norm%        a2=-lam_rec(nbroy)/stp_nrm(nbroy);        a1=1 - lam_rec(nbroy);        zz=stp(:,nbroy)'*z;        a3=a1*zz/stp_nrm(nbroy);        a4=1+a2*zz;        stp(:,nbroy+1)=(z-a3*stp(:,nbroy))/a4;        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 whileend%% We're not supposed to be here, we've taken the maximum% number of iterations and not terminated.%sol=x;it_hist=it_histx(1:itc+1,:);ierr=1;if debug==1    disp(outstat)end