function [sol, it_hist, ierr] = nsolgm(x,f,tol, parms)
% Newton-GMRES locally convergent solver for f(x) = 0
%
% Uses Eisenstat-Walker forcing term
%
% C. T. Kelley, July 1, 1994
%
% This code comes with no guarantee or warranty of any kind.
%
% function [sol, it_hist, ierr] = nsolgm(x,f,tol,parms)
%
% inputs:
% initial iterate = x
% function = f
% tol = [atol, rtol] relative/absolute
% error tolerances for the nonlinear iteration
% parms = [maxit, maxitl, etamax]
% maxit = maxmium number of nonlinear iterations
% default = 40
% maxitl = maximum number of inner iterations
% default = 40
% |etamax| = Maximum error tolerance for residual in inner
% iteration. The inner iteration terminates
% when the relative linear residual is
% smaller than eta*| F(x_c) |. eta is determined
% by the modified Eisenstat-Walker formula if etamax > 0.
% If etamax < 0, then eta = |etamax| for the entire
% iteration.
% default: etamax=.9
%
% 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 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 terminted.
%
% internal parameters:
% debug = turns on/off iteration statistics display as
% the iteration progresses
%
% Requires fdgmres.m and givapp.m
%
%
% set the debug parameter, 1 turns display on, otherwise off
%
debug=1;
%
% initialize it_hist, ierr, and set the iteration parameters
%
gamma=.9;
ierr = 0; maxit=40; lmaxit=40; etamax=.9; it_histx=zeros(maxit,2);
%
if nargin == 4
maxit=parms(1); lmaxit=parms(2); etamax=parms(3);
end
gmparms=[abs(etamax), lmaxit];
rtol=tol(2); atol=tol(1); n = length(x); fnrm=1; itc=0;
%
% evaluate f at the initial iterate
% compute the stop tolerance
%
f0=feval(f,x);
fnrm=norm(f0)/sqrt(n);
it_histx(itc+1,1)=fnrm; it_histx(itc+1,2)=0;
fnrmo=1;
stop_tol=atol + rtol*fnrm;
outstat(itc+1, :) = [itc fnrm 0 0];
%
% main iteration loop
%
while(fnrm > stop_tol & itc < maxit)
%
% keep track of the ratio (rat = fnrm/frnmo)
% of successive residual norms and
% the iteration counter (itc)
%
rat=fnrm/fnrmo;
fnrmo=fnrm;
itc=itc+1;
%
% compute the step using a GMRES routine especially designed
% for this purpose
%
[step, errstep, inner_it_count]=fdgmres(f0, f, x, gmparms);
xold=x;
x = x + step;
f0=feval(f,x);
fnrm=norm(f0)/sqrt(n);
it_histx(itc+1,1)=fnrm;
rat=fnrm/fnrmo;
outstat(itc+1, :) = [itc fnrm inner_it_count rat];
%
if debug==1
disp(outstat(itc+1,:))
end
%
% How many function evaluations did this iteration require?
%
it_histx(itc+1,2)=it_histx(itc,2)+inner_it_count+1;
if(itc == 1) it_histx(itc+1,2) = it_histx(itc+1,2)+1; end;
%
% adjust eta
%
if etamax > 0
etaold=gmparms(1);
etanew=gamma*rat*rat;
if gamma*etaold*etaold > .1
etanew=max(etanew,gamma*etaold*etaold);
end
gmparms(1)=min([etanew,etamax]);
gmparms(1)=max(gmparms(1),.5*stop_tol/fnrm);
end
% end while
end
sol=x;
it_hist=it_histx(1:itc+1,:);
if debug==1
disp(outstat)
it_hist=it_histx(1:itc+1,:);
end
%
% on failure, set the error flag
%
if fnrm > stop_tol; ierr = 1; end
