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;enderror=[];%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