function [x, error, total_iters] = ... fdcgstab(f0, f, xc, params, xinit) % Forward difference Bi-CGSTAB solver for use in nsola % % C. T. Kelley, December 29, 1994 % % This code comes with no guarantee or warranty of any kind. % % function [x, error, total_iters] % = fdcgstab(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=[]; rho=zeros(kmax+1,1); % % Use zero vector as initial iterate for Newton step unless % the calling routine has a better idea (useful for GMRES(m)). % x=zeros(n,1); r=b; if nargin == 5 x=xinit; r=-dirder(xc, x, f, f0)-f0; end % hatr0=r; k=0; rho(1)=1; alpha=1; omega=1; v=zeros(n,1); p=zeros(n,1); rho(2)=hatr0'*r; zeta=norm(r); error=[error,zeta]; % % Bi-CGSTAB iteration % while((zeta > errtol) & (k < kmax)) k=k+1; if omega==0 error('Bi-CGSTAB breakdown, omega=0'); end beta=(rho(k+1)/rho(k))*(alpha/omega); p=r+beta*(p - omega*v); v=dirder(xc,p,f,f0); tau=hatr0'*v; if tau==0 error('Bi-CGSTAB breakdown, tau=0'); end alpha=rho(k+1)/tau; s=r-alpha*v; t=dirder(xc,s,f,f0); tau=t'*t; if tau==0 error('Bi-CGSTAB breakdown, t=0'); end omega=t'*s/tau; rho(k+2)=-omega*(hatr0'*t); x=x+alpha*p+omega*s; r=s-omega*t; zeta=norm(r); total_iters=k; error=[error, zeta]; end