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

	Simple three level multi-grid program to
	solve Poisson equation. u1,r1 value/residual on fine grid
	u2,r2 correction/residual on coarse mesh.
	Most routines name explains their function.
	Program would be more elegant using some form
	of recursion to deal with V cycle. One
	way would be to define u[levels][MAX_X][MAX_Y]
	and then use the extra index to control
	V cycle level.

--------------------------------------------------*/
#include <stdio.h>
#include <math.h>
#define MAX_X		201
#define MAX_Y		201
double u1[MAX_X][MAX_Y];
double u2[MAX_X/2+1][MAX_Y/2+1];
double r1[MAX_X][MAX_Y];
double r2[MAX_X/2+1][MAX_Y/2+1];
double u3[MAX_X/4+1][MAX_Y/4+1];
double r3[MAX_X/4+1][MAX_Y/4+1];
double b[MAX_X][MAX_Y];
void initialise_ubr();
void calculate_r1();
void set_r2();
void gauss_seidel_1();
void gauss_seidel_2();
void interpolate_2_onto_1();
void interpolate_3_onto_2();
void gauss_seidel_3();
void set_r3();


int m,n,m1,n1,m2,n2,m3,n3,l,p1,p2,p3;
double hx1,hx2,hy1,hy2,hx3,hy3,stop_criterion;
double dx,dy;
double norm_r1_squared;

main()
{
int i,j,iterations;
	printf("Input number of x points:");
	scanf("%d",&m);
	m1=m-1;
	m2=m1/2;
	m3=m2/2;
	printf("Input number of y points:");
	scanf("%d",&n);
	n1=n-1;
	n2=n1/2;
	n3=n2/2;
	dx=1./(double)(m1);
	dy=1./(double)(n1);
	hx1=1./(dx*dx);
	hy1=1./(dy*dy);
	hx2=hx1/4;
	hy2=hy1/4;
	hx3=hx2/4;
	hy3=hy2/4;
	stop_criterion=0.000001*m*n;
	iterations=0;
	p1=2;
	p2=8;
	p3=16;

	initialise_ubr();
	calculate_r1();
	while(norm_r1_squared>stop_criterion)
		{
		iterations++;
		for(l=0;l<p1;l++)
			{
			gauss_seidel_1();
			}
		calculate_r1();
		set_r2();
		for(l=0;l<p2;l++)
			{
			gauss_seidel_2();
			}
		set_r3();
		for(l=0;l<p3;l++)
			{
			gauss_seidel_3();
			}
		interpolate_3_onto_2();
		for(l=0;l<p2;l++)
			{
			gauss_seidel_2();
			}

		interpolate_2_onto_1();
		for(l=0;l<p1;l++)
			{
			gauss_seidel_1();
			}

		calculate_r1();
	printf("%8d  %15.6le\n",iterations,norm_r1_squared/(m*n));
		}
	printf("Solution obtained in %8d iterations\n",iterations);
	printf("centre value %15.8le\n",u1[m1/2][n1/2]);
}
void initialise_ubr()
{
int i,j;
	for(i=0;i<=m1;i++)
		{for(j=0;j<=n1;j++)
			{u1[i][j]=0.;
			b[i][j]=-1.;
			r1[i][j]=0.;}
		}
	for(i=0;i<=m2;i++)
		{for(j=0;j<=n2;j++)
			{u2[i][j]=0.;
			r2[i][j]=0.;}
		}
	for(i=0;i<=m3;i++)
		{for(j=0;j<=n3;j++)
			{u3[i][j]=0.;
			r3[i][j]=0.;}
		}
}
void calculate_r1()
{int i,j;
	norm_r1_squared=0.;
	for(i=1;i<m1;i++)
		{
		for(j=1;j<n1;j++)
		  {
		  r1[i][j]=b[i][j]-hx1*(u1[i+1][j]-2*u1[i][j]+u1[i-1][j])
			       -hy1*(u1[i][j+1]-2*u1[i][j]+u1[i][j-1]);
		norm_r1_squared=norm_r1_squared+r1[i][j]*r1[i][j];
			}
		}

}
void set_r3()
{int i,j,i1,j1;
	
	for(i=1;i<m3;i++)
		{i1=2*i;
		for(j=1;j<n3;j++)
		  {j1=2*j;
	    u3[i][j]=0.;
	    r3[i][j]=r2[i1][j1]-hx2*(u2[i1+1][j1]-2*u2[i1][j1]+u2[i1-1][j1])
			       -hy2*(u2[i1][j1+1]-2*u2[i1][j1]+u2[i1][j1-1]);
			}
		}

}
void set_r2()
{int i,j;
	for(i=1;i<m2;i++)
		{
		for(j=1;j<n2;j++)
			{
			r2[i][j]=r1[2*i][2*j];
			u2[i][j]=0.;
			}
		}
}

void gauss_seidel_1()
{int i,j;
	for(i=1;i<m1;i++)
		{
		for(j=1;j<n1;j++)
		  {
		  u1[i][j]=(hx1*(u1[i+1][j]+u1[i-1][j])
			  +hy1*(u1[i][j+1]+u1[i][j-1])-b[i][j])/(2*hx1+2*hy1);
		  }
		}
}
void gauss_seidel_2()
{int i,j;
	for(i=1;i<m2;i++)
		{
		for(j=1;j<n2;j++)
		  {
		  u2[i][j]=(hx2*(u2[i+1][j]+u2[i-1][j])
			  +hy2*(u2[i][j+1]+u2[i][j-1])-r2[i][j])/(2*hx2+2*hy2);
		  }
		}
}
void gauss_seidel_3()
{int i,j;
	for(i=1;i<m3;i++)
		{
		for(j=1;j<n3;j++)
		  {
		  u3[i][j]=(hx3*(u3[i+1][j]+u3[i-1][j])
			  +hy3*(u3[i][j+1]+u3[i][j-1])-r3[i][j])/(2*hx3+2*hy3);
		  }
		}
}
void interpolate_2_onto_1()
{int i1,j1,i,j;
	for(i=1;i<=m2;i++)
		{
		for(j=1;j<=n2;j++)
		  {i1=2*i;
		   j1=2*j;
		  u1[i1][j1]=u1[i1][j1]+u2[i][j];
		  u1[i1-1][j1]=u1[i1-1][j1]+0.5*(u2[i-1][j]+u2[i][j]);
		  u1[i1][j1-1]=u1[i1][j1-1]+0.5*(u2[i][j-1]+u2[i][j]);
		  u1[i1-1][j1-1]=u1[i1-1][j1-1]+
			0.25*(u2[i][j]+u2[i-1][j]+u2[i][j-1]+u2[i-1][j-1]);
		  }
		}
}void interpolate_3_onto_2()
{int i1,j1,i,j;
	for(i=1;i<=m3;i++)
		{
		for(j=1;j<=n3;j++)
		  {i1=2*i;
		   j1=2*j;
		  u2[i1][j1]    =u2[i1][j1]+u3[i][j];
		  u2[i1-1][j1]  =u2[i1-1][j1]+0.5*(u3[i-1][j]+u3[i][j]);
		  u2[i1][j1-1]  =u2[i1][j1-1]+0.5*(u3[i][j-1]+u3[i][j]);
		  u2[i1-1][j1-1]=u2[i1-1][j1-1]+
			0.25*(u3[i][j]+u3[i-1][j]+u3[i][j-1]+u3[i-1][j-1]);
		  }
		}
}