using IterativeSolvers

export bicgstabl, bicgstabl!

bicgstabl(A, b, Pl=1; kwargs...) =  bicgstabl!(IterativeSolvers.zerox(A,b), A, b, Pl; kwargs...)

function bicgstabl!(x, A, b, Pl=1; tol::Real=size(A,2)*eps(), maxiter::Int=size(A,2))
    K = KrylovSubspace(A, length(b), 1, Vector{IterativeSolvers.Adivtype(A,b)}[])
    IterativeSolvers.init!(K, x)
    bicgstabell(K; tol=tol, l=1)
end

function bicgstabell(K::KrylovSubspace; tol::Real=eps(), l::Int=1)
    k=-l
    r = zeros(3)
    IterativeSolvers.initrand!(K)
    #XXX Choose r̃[0]
    r[0]=b-IterativeSolvers.nextvec(K)
    u[-1]=0
    x[0]=K.v0
    ρ₀=1
    α=0
    ω=1
    
    while norm(r[k+l]) >= tol
        k+=l
        û[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
                û[i]=r̂[i]-β*û[i]
            end
            û[j+1]=K.A*û[j]
            γ=û[j+1]⋅r̃[0]; α=α[k+j]=ρ0/γ
            for i=0:j
                r̂[i]-=α*û[i+1]
            end
            r̂[j+1]=K.A*r̂[j]; x̂0+=α*û[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]; û[0]-=γ[l]*û[l] #update
        for j=1:l-1
            û[0]-=γ[j]*û[j]
            x̂[0]+=γ″[j]*r̂[j]
            r̂[0]-=γ′[j]*r̂[j]
        end
        u[k+l-1]=û[0]; r[k+l]=r̂[0]; x[k+l]=x̂[0]
    end
    #TODO What to return?
end