/* ----------------------------------------------

	straight forward conjugate gradient solution
        of Poisson equation on a unit square with 
	zero boundary conditions
		solving del^2 x =b

------------------------------------------------*/

#include <stdio.h>
#include <math.h>
#define MAX_X		201
#define MAX_Y		201
double x[MAX_X][MAX_Y];
double x_h[MAX_X/2][MAX_Y/2];
double r[MAX_X][MAX_Y];
double r_h[MAX_X][MAX_Y];
double b[MAX_X][MAX_Y];
double d[MAX_X][MAX_Y];
double q[MAX_X][MAX_Y];
int m,n;
double dx,dy;
double dx_h,dy_h;
double alpha,dn,d0,beta,hx,hy,eps,dtq,xm;
main()
{
int i,j,iterations;
	printf("Input number of x points:");
	scanf("%d",&m);
	printf("Input number of y points:");
	scanf("%d",&n);
	dx=1./(m-1.);
	dy=1./(n-1.);
	hx=1./(dx*dx);
	hy=1./(dy*dy);
	eps=0.000001*m*n;
	iterations=0;
/* initiailise x,b,r */
	for(i=0;i<m;i++)
		{
		for(j=0;j<n;j++)
			{
			x[i][j]=0.;
			b[i][j]=-1.;
			r[i][j]=0.;
			}
		}
/* set residuals in interior */
	dn=0.;
	for(i=1;i<m-1;i++)
		{
		for(j=1;j<n-1;j++)
			{
			r[i][j]=b[i][j]-hx*(x[i+1][j]-2*x[i][j]+x[i-1][j])
			               -hy*(x[i][j+1]-2*x[i][j]+x[i][j-1]);
			d[i][j]=r[i][j];
			dn=dn+r[i][j]*r[i][j];
			}
		}
	d0=dn;
/* loop until convergence criterion reached */
	while(dn>eps)
		{
		iterations++;
		dtq=0;
/* calculate q and d.q simultaneously */
		for(i=1;i<m-1;i++)
			{
			for(j=1;j<n-1;j++)
				{
				q[i][j]=hx*(d[i+1][j]-2*d[i][j]+d[i-1][j])
			               +hy*(d[i][j+1]-2*d[i][j]+d[i][j-1]);
				dtq=dtq+d[i][j]*q[i][j];
				}
			}
/* calculate search distance alpha in direction d */
		alpha=dn/dtq;
		d0=dn;
/* update x and residual */
		dn=0.;
		for(i=1;i<m-1;i++)
			{
			for(j=1;j<n-1;j++)
				{
				x[i][j]=x[i][j]+alpha*d[i][j];
				r[i][j]=r[i][j]-alpha*q[i][j];
				dn=dn+r[i][j]*r[i][j];
				}
			}
/* calculate next search direction */
		beta=dn/d0;
		for(i=1;i<m-1;i++)
			{
			for(j=1;j<n-1;j++)
				{
				d[i][j]=r[i][j]+beta*d[i][j];
				}
			}
		printf("%6d %15.8le\n",iterations,dn/(m*n));
		}
	printf("Solution obtained in %8d iterations\n",iterations);
	printf("centre value %15.8le\n",x[(m-1)/2][(n-1)/2]);
}