function [x, error, total_iters] = ... fdtfqmr(f0, f, xc, params, xinit) % Forward difference TFQMR solver for use in nsola % % C. T. Kelley, December 30, 1994 % % This code comes with no guarantee or warranty of any kind. % % function [x, error, total_iters] % = fdtfqmr(f0, f, xc, params, xinit) % % % % Input: f0 = function at current point % f = nonlinear function % the format for f is function fx = f(x) % Note that for Newton-GMRES we incorporate any % preconditioning into the function routine. % xc = current point % params = two dimensional vector to control iteration % params(1) = relative residual reduction factor % params(2) = max number of iterations % % xinit = initial iterate. xinit=0 is the default. This % is a reasonable choice unless restarts are needed. % % % Output: x=solution % error = vector of residual norms for the history of % the iteration % total_iters = number of iterations % % Requires: dirder.m % % % initialization % b=-f0; n=length(b); errtol = params(1)*norm(b); kmax = params(2); error=[]; x=zeros(n,1); r=b; if nargin == 5 x=xinit; r=-dirder(xc, x, f, f0)-f0; end % u=zeros(n,2); y=zeros(n,2); w = r; y(:,1) = r; k=0; d=zeros(n,1); v=dirder(xc, y(:,1),f,f0); u(:,1)=v; theta=0; eta=0; tau=norm(r); error=[error,tau]; rho=tau*tau; % % TFQMR iteration % while( k < kmax) k=k+1; sigma=r'*v; % if sigma==0 error('TFQMR breakdown, sigma=0') end % alpha=rho/sigma; % % % for j=1:2 % % Compute y2 and u2 only if you have to % if j==2 y(:,2)=y(:,1)-alpha*v; u(:,2)=dirder(xc, y(:,2),f,f0); end m=2*k-2+j; w=w-alpha*u(:,j); d=y(:,j)+(theta*theta*eta/alpha)*d; theta=norm(w)/tau; c=1/sqrt(1+theta*theta); tau=tau*theta*c; eta=c*c*alpha; x=x+eta*d; % % Try to terminate the iteration at each pass through the loop % if tau*sqrt(m+1) <= errtol error=[error, tau]; total_iters=k; return end end % % % if rho==0 error('TFQMR breakdown, rho=0') end % rhon=r'*w; beta=rhon/rho; rho=rhon; y(:,1)=w + beta*y(:,2); u(:,1)=dirder(xc, y(:,1),f,f0); v=u(:,1)+beta*(u(:,2)+beta*v); error=[error, tau]; total_iters=k; end