Skip to content

Instantly share code, notes, and snippets.

@jiahao
Last active October 12, 2017 19:51
Show Gist options
  • Save jiahao/7349818 to your computer and use it in GitHub Desktop.
Save jiahao/7349818 to your computer and use it in GitHub Desktop.
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)
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment