//      program laplace

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

const int imax=500; //2001;
const int jmax=500; //2001;

double max (double x, double y)
{
   if (x >= y)
      return x;
   else return y;
}

int main()
{
    int im1,im2,jm1,jm2,it,itmax=100;
    double u[imax][jmax];
    double du[imax][jmax];
    double umax=10.0, dumax, tol=1.0e-2; //tol=1.0e-3;
    int i, j;
    struct timeval tv;
    struct timezone tz;
    long starttimeSec, starttimeMicro;
    long endtimeSec, endtimeMicro;
    long totalSec, totalMicro;   
     
//    parameter (imax=2001,jmax=2001)
//    parameter (im1=imax-1,im2=imax-2,jm1=jmax-1,jm2=jmax-2)
//    parameter (itmax=100)
//    real*8 u(imax,jmax),du(imax,jmax),umax,dumax,tol
//    parameter (umax=10.0,tol=1.0e-3)
      
      im1=imax-1;
      im2=imax-2;
      jm1=jmax-1;
      jm2=jmax-2;
      
//! Initialize
      for (j=0; j < jmax; j++)
      {
         for (i=0; i < imax-1; i++)
	 {
            u[i][j] = 0.0;
            du[i][j] = 0.0;
         }
         u[imax-1][j] = umax;
      }

      dumax = 1000.0;
      it=1;
      
     gettimeofday(&tv, &tz);
     starttimeSec = tv.tv_sec;
     starttimeMicro = tv.tv_usec;
     
//! Main computation loop
      while (dumax >= tol)
      {
         dumax = 0.0;
         for (j=1; j < jm1; j++)
	 {
            for (i=1; i < im1; i++)
	    {
               du[i][j]=0.25*(u[i-1][j] + u[i+1][j] +
	                u[i][j-1] + u[i][j+1]) - u[i][j];
	    }
	 }
         for (j=1; j < jm1; j++)
	 {
            for (i=1; i < im1; i++)
	    {
               dumax = max(dumax, fabs(du[i][j]));
	    }
	 }
         for (j=1; j < jm1; j++)
	 {
            for (i=1; i < im1; i++)
	    {
               u[i][j] = u[i][j] + du[i][j];
            }
         }
//         printf("* %d  %f\n", it,dumax);
         it=it+1;
      }
      
          gettimeofday(&tv, &tz);
     endtimeSec = tv.tv_sec;
     endtimeMicro = tv.tv_usec;
     if (endtimeSec > starttimeSec)
     {
        totalSec = endtimeSec-starttimeSec;
	totalMicro = endtimeMicro + (1000000 - starttimeMicro);
     }
     else {
        totalSec = 0;
	totalMicro = endtimeMicro - starttimeMicro;
     }	
 
  
       printf("* %d  %f\n", it,dumax);
   printf("Time = %d sec, %d microsec\n", totalSec, totalMicro);        
       
      return 0;
}