/*
  solving discrete poisson equation with Jacobi iteration
*/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define nx 32
#define ny 32

int main()
{
    int i, j, iter, itmax=1000;
    double a=1.0, b=1.0, tol=1.0e-6;
    double hx, hy, d, dx, dy, xi, yj, res, rtmp, err, tc;
    double (*u)[ny+1], (*ut)[ny+1], (*unew)[ny+1], (*f)[ny+1];

    u    = malloc((nx+1)*(ny+1)*sizeof(double));
    ut   = malloc((nx+1)*(ny+1)*sizeof(double));
    unew = malloc((nx+1)*(ny+1)*sizeof(double));
    f    = malloc((nx+1)*(ny+1)*sizeof(double));
    if (!u | !ut | !unew | !f) 
	{
        printf("Out of memory!\n"); exit(1);
    }

//  initial data
    hx = a / nx;
    hy = b / ny;
    d  = hx*hx * hy*hy / ( 2*(hx*hx + hy*hy) );
    dx = d / (hx*hx);
    dy = d / (hy*hy);
  
//  true solution : u = - (x^2 + y^2)/4
    for (i=0; i<=nx; i++) 
	{
	    for (j=0; j<=ny; j++) 
		{
		    xi = hx * i;
		    yj = hy * j;
		    ut[i][j] = -0.25*(xi*xi + yj*yj);
	    }
    }
  
//  set initial guess and right hand side
    for (i=1; i<nx; i++) 
	    for (j=1; j<ny; j++) 
		{
		    u[i][j] = 0.0;
		    f[i][j] = d;
	    }

//  boundary condition
    for (i=0; i<=nx; i++)  
    {
	  	u[i][0]  = ut[i][0];
	  	u[i][ny] = ut[i][ny];
    }
    for (j=0; j<=ny; j++)  
    {
	  	u[0][j]  = ut[0][j];
	  	u[nx][j] = ut[nx][j];
    }

//  ---------------------------------
//  Jacobi iteration
    tc = clock()/CLOCKS_PER_SEC;
    for (iter=1; iter<=itmax; iter++)
    {
        for (i=1; i<nx; i++)
        for (j=1; j<ny; j++)
	        unew[i][j] = f[i][j] + dx * (u[i+1][j] + u[i-1][j])
                                 + dy * (u[i][j+1] + u[i][j-1]);
  
//  update approximation
        for (i=1; i<nx; i++)
        for (j=1; j<ny; j++)
            u[i][j] = unew[i][j];

//  check residual
        res = 0.0;
        for (i=1; i<nx; i++)
        for (j=1; j<ny; j++) 
		{
	        rtmp = f[i][j] - u[i][j] 
		           + dx * (u[i+1][j] + u[i-1][j])
                   + dy * (u[i][j+1] + u[i][j-1]);
            if (res < fabs(rtmp))  res = fabs(rtmp);
        }
        res = res / d;
        printf("iter=%3d, residual=%.4e\n", iter, res);
	    if (res < tol)  break;
    }
    tc = clock()/CLOCKS_PER_SEC - tc;
    printf("Total iteration=%3d, residual=%.4e, time=%.4e\n", iter-1, res, tc);

//  check solution
    err = 0.0;
    for (i=1; i<nx; i++)
    for (j=1; j<ny; j++) 
	{
       rtmp = fabs(ut[i][j] - u[i][j]);
       if (err < rtmp)  err = rtmp;
    }
    printf("||u-ut|| = %.4e\n", err);

	free(u); free(ut); free(unew); free(f);
	
	return 0;
}