//This is a code for solving the Poisson equation with
//homogeneous boundary conditions 
// A subroutine called proc_cg.c needs to be used (a CG code)
// A description of this code can be found in a PDF file (in Chinese)
// titled multigrid.pdf

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "time.h"

#define cengshu 7   //最细的网格的是第cengshu 层
#define N 128        //最细的网格剖分是N*N
#define minN 16     //最粗的网格剖分是minN*minN

double newz[(N+1)*(N+1)],newg[(N+1)*(N+1)];
double error[(N+1)*(N+1)],exact[(N+1)*(N+1)];

#include "proc_cg.c"

void main()
{ int i,j;
  double h,maxerror;
double sum;
  int k,cishu;
  clock_t st,ct;

  extern void cg();
void gauss();
  void limit();
  void extend();
  void MG();

h=1.0/N;

  j=0; 
  for(i=0;i<(N+1)*(N+1);i++)
  {newz[i]=0;
   if((i%(N+1))==0&&i!=0) j++;
   newg[i]=-2*(j*h)*(j*h-1)-2*(i%(N+1)*h)*(i%(N+1)*h-1);         //给出右端项
}
 
  cishu=0;
  st=clock();

  do
  {sum=0;
   MG(cengshu,newz,newg);           //一次v-cyle

  for(i=1;i<N;i++)
  for(j=1;j<N;j++)
  {k=(N+1)*i+j;
    sum+=(newz[k+1]+newz[k-1]+newz[k-N-1]+newz[k+N+1]+h*h*newg[k]-4.0*newz[k])*
		    (newz[k+1]+newz[k-1]+newz[k-N-1]+newz[k+N+1]+h*h*newg[k]-4.0*newz[k]);
  } 
  cishu++;
 } while(sqrt(sum)>1.0e-10);
 
 ct=clock();
 printf("%e\n",(double) (ct-st) / CLOCKS_PER_SEC);

 j=0;
 for(i=0;i<(N+1)*(N+1);i++)
 {if((i%(N+1))==0&&i!=0) j++;
  exact[i]=(i%(N+1)*h)*(j*h)*(j*h-1)*(i%(N+1)*h-1);
  error[i]=fabs(exact[i]-newz[i]);
}

maxerror=0;
 for(i=0;i<(N+1)*(N+1);i++)
 {if(maxerror<error[i]) maxerror=error[i];}

 printf("The maximum error is : %e",maxerror);
}

void MG(int k,double *z,double *g)
{int i,j,kk;
 double h;
 double *r,*rr,*q,*qq,zz[minN+1][minN+1],gg[minN+1][minN+1];
 int n;                       //剖分次数
 int it;

 n=(int) pow(2,k);
 h=1.0/n;
 
if(k==4)                   //最粗的网格是16*16剖分的网格，在上面用共轭梯度法求解
 {for(i=0;i<(minN+1)*(minN+1);i++)
     z[i]=0;	 
	 
  it=0;
  for(i=0;i<(minN+1);i++)
  {for(j=0;j<(minN+1);j++)
  {zz[i][j]=z[it*(minN+1)+j];
   gg[i][j]=g[it*(minN+1)+j];
  }
   it++;
  }
 
  cg(zz,gg);

  j=0;
  for(i=0;i<(minN+1)*(minN+1);i++)
  {if((i%(minN+1))==0&&i!=0) j++;
   z[i]=zz[j][i%(minN+1)];
  }
 }

 else 
 {if(k!=cengshu){
	 for(i=0;i<(n+1)*(n+1);i++)
		 z[i]=0;}

  gauss(z,g,n);                     //Gauss 迭代１次求解
  
  r=malloc((n+1)*(n+1)*sizeof(double));
  for(i=0;i<=n;i++)
	  for(j=0;j<=n;j++)
	  {kk=(n+1)*i+j;
	   r[kk]=0;
	  }

  for(i=1;i<n;i++)
	  for(j=1;j<n;j++)
	  {kk=(n+1)*i+j;
       r[kk]=g[kk]-(4.0*z[kk]-z[kk+1]-z[kk-1]-z[kk+n+1]-z[kk-n-1])/(h*h);
	  }

  rr=malloc(((n/2)+1)*((n/2)+1)*sizeof(double));
  limit(rr,r,n);                             //限制余量到粗网格上，用到限制算子
  q=malloc(((n/2)+1)*((n/2)+1)*sizeof(double));

  MG(k-1,q,rr);                           //在k-1层网格上调用MG求解余量方程

  qq=malloc((n+1)*(n+1)*sizeof(double));
  extend(qq,q,n/2);                        //把余量方程求得的解延拓到细网格上

  for(i=0;i<(n+1)*(n+1);i++)
	  z[i]+=qq[i];
   
  gauss(z,g,n);                           //以新的初值在细网格上进行迭代求解
 }
}

void limit(double *rr,double *r,int n)
{int i,j,k1,k2;
 
 for(i=0;i<=n/2;i++)
	 for(j=0;j<=n/2;j++)
	 {k1=(n/2+1)*i+j;
	  rr[k1]=0;
	 }

 for(i=1;i<n/2;i++)
	 for(j=1;j<n/2;j++)
	 {k1=(n/2+1)*i+j;
	  k2=2*(n+1)*i+2*j;
	  rr[k1]=(1.0/4)*r[k2]+(1.0/8)*r[k2-1]+(1.0/8)*r[k2+1]+
		     (1.0/8)*r[k2+n+1]+(1.0/8)*r[k2-n-1]+
			 (1.0/16)*r[k2+n]+(1.0/16)*r[k2+n+2]+(1.0/16)*r[k2-n]+(1.0/16)*r[k2-n-2];
	 }
}


void extend(double *rr,double *r,int n)
{int i,j,k1,k2,k,m;
 
 for(i=0;i<=2*n;i++)
	 for(j=0;j<=2*n;j++)
	 {k1=(2*n+1)*i+j;
	  rr[k1]=0;
	 }

 for(i=1;i<n;i++)
	 for(j=1;j<n;j++)
	 {k1=(n+1)*i+j;
	  k2=2*(2*n+1)*i+2*j;
	  rr[k2]=r[k1];
	 }

 m=1;
 for(k=(2*n+1)+2;k<=(2*n+1)*(2*n-1)+(2*n-1);k=k+2)
 {m++;
  rr[k]=(1.0/2)*(rr[k+2*n+1]+rr[k-2*n-1]);
  if(m%n==0) {k=k+2*(2*n);m=1;}
 }

 m=0;
 for(k=(2*n+1)+1;k<=(2*n+1)*(2*n-1)+(2*n-1);k=k+2)
 {m++;
  rr[k]=(1.0/2)*(rr[k+1]+rr[k-1]);
  if(m%n==0) {k++;m=0;}
 }
}

void gauss(double *z,double *g,int n)
{double h;
 int i,j,k,m;
 h=1.0/n;

 for(m=0;m<1;m++)
{for(i=1;i<n;i++)
  for(j=1;j<n;j++)
	  {k=(n+1)*i+j;
       z[k]=(z[k+1]+z[k-1]+z[k-n-1]+z[k+n+1]+h*h*g[k])/4.0;
	  }
 }
}