CG.f90

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

C version above has been translated from fortran by Adrien Cassagne.

CG.c

#include<math.h>
#include<stdio.h>
#include<stdlib.h>
 
int main(int argc, char** argv)
{
// 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
 
	int    i, n, Nx, ndim = 1;
	double Lx, Dx, InvDxDx;
	double beta, rho, rho0, gam, maxResidut;
	double *T, *Residut, *p, *q; // pointers
	FILE   *f;
 
	Nx      = 100;
	Lx      = 0.1;
	Dx      = Lx / ((Nx -1) * 1.0);
	InvDxDx = 1.0 / (Dx * Dx);
 
	T       = (double*)malloc((Nx +1) * sizeof(double));
	Residut = (double*)malloc((Nx +1) * sizeof(double));
	p       = (double*)malloc((Nx +1) * sizeof(double));
	q       = (double*)malloc((Nx +1) * sizeof(double));
 
	// Initial Solution
	T[0] = 10;
	for(i = 1; i < Nx; ++i)
		T[i] = 0.0;
	T[Nx] = -1.0;
 
	// Initialisation
	for(i = 1; i < Nx; ++i)
		Residut[i] = -( T[i +1] - 2.0 * T[i] + T[i -1] ) * InvDxDx;
 
	for(i = 0; i < Nx +1; ++i)
		p[i] = 0.0;
 
	beta = 0.0;
	rho  = 0.0;
 
	for(i = 1; i < Nx; ++i)
		rho += Residut[i] * Residut[i];
 
	// Iteration
	n = 0;
	while(1)
	{
		n++;
		for(i = 1; i < Nx; ++i)
			p[i] = Residut[i] + beta * p[i];
 
		// Matrix Operation : q = A p
		for(i = 1; i < Nx; ++i)
			q[i] = (p[i +1] - 2.0 * p[i] + p[i -1]) * InvDxDx;
 
		gam = 0.0;
		for(i = 1; i < Nx; ++i)
			gam = gam + p[i] * q[i];
		gam = rho / gam;
 
		for(i = 1; i < Nx; ++i)
			T[i] = T[i] + gam * p[i];
 
		for(i = 1; i < Nx; ++i)
			Residut[i] = Residut[i] - gam * q[i];
 
		// find maxResidut
		maxResidut = 0;
		for(i = 1; i < Nx; ++i)
			if(fabs(Residut[i]) > maxResidut)
				maxResidut = fabs(Residut[i]);
 
		// Test if residut is < 1.e-7 so that solution is converged
		if(maxResidut < 1e-7)
			break;  // Converged Solution
 
		printf("CG step %d, maxResidut %lf.\n", n, maxResidut);
 
		rho0 = rho; // Save old rho value
		rho  = 0.0;
 
		for(i = 1; i < Nx; ++i)
			rho += Residut[i] * Residut[i];
		beta = rho / rho0;
	}
 
	printf("Solution converged in %d iterations, with a residut of %g\n", n, maxResidut);
 
	// Solution plot
	f = fopen("out-GC.dat", "w");
	for(i = 1; i < Nx; ++i)
		fprintf(f, "%lf %lf\n", Dx * (i -1), T[i]);
	fclose(f);
 
	free(T);
	free(Residut);
	free(p);
	free(q);
 
	return 0;
}