function [sol, it_hist, ierr] = nsola(x,f,tol, parms) % Newton-Krylov solver, globally convergent % solver for f(x) = 0 % % Inexact-Newton-Armijo iteration % % Eisenstat-Walker forcing term % % Parabolic line search via three point interpolation. % % C. T. Kelley, June 29, 1994 % parms made truly optional, Feb 12, 1996 % % This code comes with no guarantee or warranty of any kind. % % function [sol, it_hist, ierr] = nsola(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, lmeth, restart_limit] % maxit = maxmium number of nonlinear iterations % default = 40 % maxitl = maximum number of inner iterations before restart % in GMRES(m), m=maxitl % default = 40 % % For iterative methods other than GMRES(m) maxitl % is the upper bound on linear iterations. % % |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 % % lmeth = choice of linear iterative method % 1 (GMRES), 2 GMRES(m), % 3 (BICGSTAB), 4 (TFQMR) % default = 1 (GMRES, no restarts) % % restart_limit = max number of restarts for GMRES if % lmeth = 2 % default=20 % % output: % sol = solution % it_hist(maxit,3) = scaled l2 norms of nonlinear residuals % for the iteration, number of 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. % ierr = 2 failure in the line search. The iteration % is terminated if too many steplength reductions % are taken. % % % internal parameters: % debug = turns on/off iteration statistics display as % the iteration progresses % % alpha = 1.d-4, parameter to measure sufficient decrease % % sigma0 = .1, sigma1=.5, safeguarding bounds for the linesearch % % maxarm = 50, maximum number of steplength reductions before % failure is reported % % % Requires dirder.m, fdkrylov.m, parab3p.m, fdgmres.m, % fdcgstab.m, fdtfqmr.m, givapp.m % % % % % set the debug parameter, 1 turns display on, otherwise off % debug=0; % % set alpha, sigma0, sigma1, maxarm, and restart_limit % alpha = 1.d-4; sigma0=.1; sigma1=.5; maxarm = 50; % % initialize it_hist, ierr, and set the iteration parameters % gamma=.9; ierr = 0; maxit=40; lmaxit=40; etamax=.9; it_histx=zeros(maxit,3); lmeth=1; restart_limit=20; % % initialize parameters for the iterative methods % gmparms=[abs(etamax), lmaxit]; if nargin == 4 maxit=parms(1); lmaxit=parms(2); etamax=parms(3); it_histx=zeros(maxit,3); gmparms=[abs(etamax), lmaxit]; if length(parms)>=4 lmeth=parms(4); end if length(parms)==5 gmparms=[abs(etamax), lmaxit, parms(5), 1]; end end % 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; it_histx(itc+1,3)=0; fnrmo=1; stop_tol=atol + rtol*fnrm; outstat(itc+1, :) = [itc fnrm 0 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; [step, errstep, inner_it_count,inner_f_evals]=... fdkrylov(f0, f, x, gmparms, lmeth); % % The line search starts here. % xold=x; lambda=1; lamm=1; lamc=lambda; iarm=0; xt = x + lambda*step; ft=feval(f,xt); nft=norm(ft); nf0=norm(f0); ff0=nf0*nf0; ffc=nft*nft; ffm=nft*nft; while nft >= (1 - alpha*lambda) * nf0; % % apply the three point parabolic model % if iarm == 0 lambda=sigma1*lambda; else lambda=parab3p(lamc, lamm, ff0, ffc, ffm); end % % update x; keep the books on lambda % xt=x+lambda*step; lamm=lamc; lamc=lambda; % % keep the books on the function norms % ft=feval(f,xt); nft=norm(ft); ffm=ffc; ffc=nft*nft; iarm=iarm+1; if iarm > maxarm disp(' Armijo failure, too many reductions '); ierr=2; disp(outstat) it_hist=it_histx(1:itc+1,:); sol=xold; return; end end x=xt; f0=ft; % % end of line search % fnrm=norm(f0)/sqrt(n); it_histx(itc+1,1)=fnrm; % % How many function evaluations did this iteration require? % it_histx(itc+1,2)=it_histx(itc,2)+inner_f_evals+iarm+1; if(itc == 1) it_histx(itc+1,2) = it_histx(itc+1,2)+1; end; it_histx(itc+1,3)=iarm; % rat=fnrm/fnrmo; % % adjust eta as per Eisenstat-Walker % 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 % outstat(itc+1, :) = [itc fnrm inner_it_count rat iarm]; % % if debug==1 % disp(outstat(itc+1,:)) % 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