#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 %lf\n", n, maxResidut);

	// Solution plot
	f = fopen("out-GC.dat", "w");
	for(i = 1; i < Nx; ++i)
		fprintf(f, "%d %lf", Dx * (i -1), T[i]); 
	fclose(f);

	free(T);
	free(Residut);
	free(p);
	free(q);

	return 0;
}