...

function [x, error, total_iters] = ...
                     tfqmr(x0, b, atv, params)
% TFQMR solver for linear systems
%
% C. T. Kelley, December 28, 1994
%
% This code comes with no guarantee or warranty of any kind.
%
% function [x, error, total_iters]
%                    = tfqmr(x0, b, atv, params)
%
%
% Input:        x0=initial iterate
%               b=right hand side
%               atv, a matrix-vector product routine
%			atv must return Ax when x is input
%			the format for atv is
%                       function ax = atv(x)
%               params = two dimensional vector to control iteration
%                        params(1) = relative residual reduction factor
%                        params(2) = max number of iterations
%
% Output:       x=solution
%               error = vector of iteration residual norms
%               total_iters = number of iterations
%
% 
%

%
% initialization
%
n=length(b); errtol = params(1)*norm(b); kmax = params(2); error=[]; x=x0;
%
if norm(x)~=0
    r=b-feval(atv,x);
else
    r=b;
end
error=[];
%
u=zeros(n,2); y=zeros(n,2); w = r; y(:,1) = r; 
k=0; d=zeros(n,1); v=feval(atv,y(:,1)); 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)=feval(atv,y(:,2));
        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)=feval(atv,y(:,1));
    v=u(:,1)+beta*(u(:,2)+beta*v);
    error=[error, tau];
    total_iters=k;
end