function biConjugateGradientStabilized(A,b,x0,tol,maxIts) result(x)
class(sparse_t),intent(in)::A
real(wp),dimension(:),intent(in)::b
real(wp),dimension(:),intent(in),optional::x0
real(wp),intent(in),optional::tol
integer,intent(in),optional::maxIts
real(wp),dimension(:),allocatable::x
real(wp),dimension(:),allocatable::D,r,Ap,p,r0,s,As
real(wp)::orho,rho,alpha,beta,omega
integer::lmaxIts
real(wp)::ltol
real(wp)::rss0,rss
integer::k
lmaxIts = 1000000
ltol = 1.0E-6_wp
if(present(maxIts)) lmaxIts = maxIts
if(present(tol )) ltol = tol
D = A%getDiagonal()
x = b/D
if(present(x0)) x = x0
r = b-matmul(A,x)
r0 = r
rho = 1.0_wp
alpha = 1.0_wp
omega = 1.0_wp
Ap = [(0.0_wp,k=1,size(b))]
p = [(0.0_wp,k=1,size(b))]
rss0 = sqrt(sum(r**2))
call startReport('Bi-Conjugate Gradient (Stab) Solver')
do k=1,lmaxIts
orho = rho
rho = sum(r0*r)
beta = (rho/orho)*(alpha/omega)
p = r+beta*(p-omega*Ap)
Ap = matmul(A,p)
alpha = rho/sum(r0*Ap)
s = r-alpha*Ap
As = matmul(A,s)
omega = sum(As*s)/sum(As*As)
x = x+alpha*p+omega*s
r = s-omega*As
rss = sqrt(sum(r**2))
call solverProgress(k,rss/rss0,ltol,'Bi-Conjugate Gradient (Stab) Solver')
if(rss/rss0<ltol) exit
end do
end function biConjugateGradientStabilized