function [x, error, total_iters] = ...
                     bicgstab(x0, b, atv, params)
% Bi-CGSTAB solver for linear systems
%
% C. T. Kelley, December 16, 1994
%
% This code comes with no guarantee or warranty of any kind.
%
% function [x, error, total_iters]
%                    = bicgstab(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;
rho=zeros(kmax+1,1);
%
if norm(x)~=0
    r=b-feval(atv,x);
else
    r=b;
end
error=[];
%
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=feval(atv,p);
    tau=hatr0'*v;
    if tau==0
        error('Bi-CGSTAB breakdown, tau=0');
    end 
    alpha=rho(k+1)/tau;
    s=r-alpha*v; 
    t=feval(atv,s);
    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