program CG

	! Conjugate Gradient
	! Solve 1D heat equation
	! Solve A T = 0, with A the matrix corresponding to the laplacian
	!
	! Heat equation : dT/dt = alpha d^2 T/dx^2
	! When converged, dT/dt = 0, the equation becomes : 0 = alpha d^2 T/dx^2 <--> 0 = d^2 T/dx^2
	! Numericaly, this means : (T(i+1)-2.0*T(i)+T(i-1)) * 1/(Dx*Dx) = 0 <--> A T = 0 with the matrix A, a diagonal symetric positiv definit matrix
	! The matrix is like :
	! -2/(dx*dx)   1/(dx*dx)     0           0           0          0          ....
	!  1/(dx*dx)   -2/(dx*dx)    1/(dx*dx)   0           0          0          ....
	!  0            1/(dx*dx)	-2/(dx*dx)   1/(dx*dx)   0          0          ....
	!  0		    0            1/(dx*dx)  -2/(dx*dx)   1/(dx*dx)  0          ....
	!  0		    0            0           1/(dx*dx)  -2/(dx*dx)  1/(dx*dx)  0
	! etc

	implicit none

	integer, parameter :: ndim = 1
	real(8), dimension(1:ndim) :: Lx,Dx,InvDxDx
	integer, dimension(1:ndim) :: Nx
	real(8) :: beta, rho, rho0, gam
	integer :: n,i
	real(8), allocatable, dimension(:) :: T,Residut,p,q

	Nx(1) = 100

	Lx(:) = 0.1d0
	Dx(:) = Lx(:)/((Nx(:)-1)*1.0d0)
	InvDxDx(:) = 1.0d0/(Dx(:)*Dx(:))

	allocate(T(0:Nx(1)+1),Residut(0:Nx(1)+1),p(0:Nx(1)+1),q(0:Nx(1)+1))

	! Initial Solution
	T(1:Nx(1)) = 0.0d0
	T(0) = 10.0d0
	T(Nx(1)+1) = -1.0d0

	! Initialisation

	do i=1,Nx(1)
		Residut(i) = - (T(i+1)-2.0d0*T(i)+T(i-1)) * InvDxDx(1)
	end do

	p(:) = 0.0d0
	beta = 0.0d0

	rho = 0.0d0
	do i=1,Nx(1)
		rho = rho + Residut(i)*Residut(i)
	end do

	! Iteration

	n = 0
	do
		n = n + 1

		do i=1,Nx(1)
			p(i) = Residut(i) + beta * p(i)
		end do

		! Matrix Operation : q = A p
		do i=1,Nx(1)
			q(i) = (p(i+1)-2.0d0*p(i)+p(i-1)) * InvDxDx(1)
		end do

		gam = 0.0d0
		do i=1,Nx(1)
			gam = gam + p(i)*q(i)
		end do
		gam = rho / gam

		do i=1,Nx(1)
			T(i) = T(i) + gam * p(i)
		end do

		do i=1,Nx(1)
			Residut(i) = Residut(i) - gam * q(i)
		end do

		! Test if residut is < 1.e-7 so that solution is converged
		if(maxval(dabs(Residut(1:Nx(1)))) < 1d-7) exit ! Converged Solution

		write(*,*) "CG step ",n,"Residut",maxval(dabs(Residut(1:Nx(1))))

		rho0 = rho ! Save old rho value
		rho = 0.0d0
		do i=1,Nx(1)
			rho = rho + Residut(i)*Residut(i)
		end do

		beta = rho / rho0

	end do

	write(*,*) "Solution converged in ", n ,"Iterations, with a res of ", maxval(abs(Residut(1:Nx(1))))

	! Solution plot

	open(10,file="out-GC.dat")
	do i=1,Nx(1)
		write(10,*) Dx(1)*(i-1),T(i)
	end do
	close(10)

	deallocate(T,Residut,p,q)

end program CG