Translation from C.T. Kelley's Matlab code into Python (original http://www.mathworks.com/matlabcentral/fileexchange/2198-iterative-methods-for-linear-and-nonlinear-equations/content/kelley/bicgstab.m)

bicgstab.py

def bicgstab(x0, b, atv, eps, kmax):
    n = b.shape[0]
    errtol = eps*np.linalg.norm(b)
    error = []
    x = x0
    
    rho = np.zeros(kmax+1)
    if np.isclose(np.linalg.norm(x), 0.):
        r = b
    else:
        r = b - atv(x)
        
    hatr0=r
    k=-1.
    rho[0]=1.
    alpha=1.
    omega=1.
    
    v=np.zeros(n)
    p=np.zeros(n)
    rho[1]=hatr0.T.dot(r)
    
    zeta=np.linalg.norm(r)
    error.append(zeta)
    
    # Bi-CGSTAB iteration
    while((zeta > errtol) and (k < kmax)):
        k=k+1
        if np.isclose(omega, 0.):
           print 'Bi-CGSTAB breakdown, omega=0'
        
        beta=(rho[k+1]/rho[k])*(alpha/omega)
        
        p=r+beta*(p - omega*v)
        
        v = atv(p)       
        tau=hatr0.T.dot(v)
        if np.isclose(tau, 0.):
            print 'Bi-CGSTAB breakdown, tau=0'
         
        alpha=rho[k+1]/tau
        
        s=r-alpha*v 
        t=atv(s)
        
        tau=t.T.dot(t)
        if np.isclose(tau, 0):
            print 'Bi-CGSTAB breakdown, t=0'
        
        omega=t.T.dot(s)/tau
        
        rho[k+2]=-omega*(hatr0.T.dot(t))
                         
        x=x+alpha*p+omega*s
        r=s-omega*t
        zeta=np.linalg.norm(r)
        
        total_iters=k
        error.append(zeta)
        
    return x, error, total_iters