A transcription of the BiCGSTAB(ell) algorithm (Sleijpen and Fokkema, 1993, http://igitur-archive.library.uu.nl/math/2006-1214-210728/sleijpen_93_bicgstab.pdf)

bicgstabell.jl

function bicgstabell(K::KrylovSubspace; tol::Real, l::Int=1)
    k=-l
    initrand!(K)
    #XXX Choose r̃[0]
    r[0]=b-nextvec(K)
    u[-1]=0
    x[0]=K.v0
    ρ₀=1
    α=0
    ω=1
    
    while norm(r[k+l]) >= tol
        k+=l
        û[0]=u[k-1]; r̂[0]=r[k]; x̂[0]=x[k]
        ρ₀*=-ω
        for j=0:l-1 #BiCG part
            ρ₁=r̂[j]⋅r̃[0]; β=β[k+j]=α*ρ1/ρ₀; ρ₀=ρ₁
            for i=0:j
                û[i]=r̂[i]-β*û[i]
            end
            û[j+1]=K.A*û[j]
            γ=û[j+1]⋅r̃[0]; α=α[k+j]=ρ0/γ
            for i=0:j
                r̂[i]-=α*û[i+1]
            end
            r̂[j+1]=K.A*r̂[j]; x̂0+=α*û[0]
        end
        for j=1:l #Mod G-S; MR part
            for i=1:j-1
                τ[i,j]=r̂[j]⋅r̂[i]/σ[i]
                r̂[j]-=τ[i,j]*r̂[i]
            end
            σ[j]=r̂[j]⋅r̂[j];γ′[j]=r̂[0]⋅r̂[j]/σ[j]
        end
        ω=γ[l]=γ′[l]
        for j=l-1:-1:1 #γ=T\γ′
            γ[j]=γ′[j]-sum([τ[j,i]*γ[i] for i=j+1:l])
        end
        for j=1:l-1 #γ″=TSγ
            γ″[j]=γ[j+1]-sum([τ[j,i]*γ[i+1] for i=j+1:l-1])
        end
        x̂[0]+=γ[1]*r̂[0]; r̂[0]-=γ′[l]*r̂[l]; û[0]-=γ[l]*û[l] #update
        for j=1:l-1
            û[0]-=γ[j]*û[j]
            x̂[0]+=γ″[j]*r̂[j]
            r̂[0]-=γ′[j]*r̂[j]
        end
        u[k+l-1]=û[0]; r[k+l]=r̂[0]; x[k+l]=x̂[0]
    end
    #TODO What to return?
end