function conjugateGradient(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
real(wp)::alpha,beta,rTr,rTr_old
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)
p = r
rTr = sum(r*r)
rss0 = sqrt(sum(r**2))
call startReport('Conjugate Gradient Solver')
do k=1,lmaxIts
Ap = matmul(A,p)
alpha = rTr/sum(p*Ap)
x = x+alpha*p
r = r-alpha*Ap
rTr_old = rTr
rTr = sum(r*r)
beta = rTr/rTr_old
p = r+beta*p
rss = sqrt(sum(r**2))
call solverProgress(k,rss/rss0,ltol,'Conjugate Gradient Solver')
if(rss/rss0<ltol) exit
end do
end function conjugateGradient