#include<stdio.h>
#include<stdlib.h>
#include<math.h>

/* This function prints the matrix A to the screen. */
void print_matrix(double **A, int imax, int jmax) {
  int i, j;
  for(i=0; i<imax; i++) {
    for(j=0; j<jmax; j++){
      printf("% 6.2e  ",A[i][j]);
    }
    printf("\n");
  }
  printf("\n");
}


/* This function prints the matrix A to the harddrive. */
void print_matrix_disc(double **A, int imax, int jmax, 
                       double dx, double dy, int k) {
  int i,j;
  char file[500];
  FILE *fp;

  sprintf(file,"p.%.4d.dat",k);

  fp = fopen(file,"w");

  for(i=1; i<=imax; i++) {
    for(j=1; j<=jmax; j++){
      fprintf(fp,"%6.2e %6.2e %6.2e\n",(i-0.5)*dx,(j-0.5)*dy,A[i][j]);
    }
  }
  fclose(fp);
}


int main(int argc, char **argv){

  double **p,**d;
  double dx,dy;
  double res,eps,tmp;
  int i,j,imax,jmax;
  int k,k_stop;

  imax = 24; jmax = 16;
  eps = 1.0e-10;
  dx = 1./(float)imax; dy = 1./(float)jmax;
  k_stop = 50;

  p = (double**)malloc((imax+2)*sizeof(double*));
  d = (double**)malloc((imax+2)*sizeof(double*));

  for (i=0; i<imax+2; i++){
    p[i] = (double*)malloc((jmax+2)*sizeof(double));
    d[i] = (double*)malloc((jmax+2)*sizeof(double));
  }

  /* Initialisation: */

  for(i=0; i<imax+2; i++){
    for(j=0; j<jmax+2; j++){
      d[i][j] = 0.;
      if( i==0 || i==imax+1 || j==0 || j==jmax+1 )
        p[i][j] = 0.;
      else
        p[i][j] = 1.;
    }
  }

  res = 10.*eps;
  k = 0;

  /* the iteration loop */

  while( k<=k_stop && res>eps ){

    print_matrix_disc(p,imax,jmax,dx,dy,k);

    res = 0.;

    /* Set boundary conditions. */
    for(i=1; i<=imax; i++){
      p[i][0]      = -p[i][1];
      p[i][jmax+1] = -p[i][jmax];
    }
    for(j=1; j<=jmax; j++){
      p[0][j]      = -p[1][j];
      p[imax+1][j] = -p[imax][j];
    }

    for(i=1; i<=imax; i++){
      for(j=1; j<=jmax; j++){
        p[i][j] = -1./ (2./(dx*dx)+2./(dy*dy)) * ( d[i][j] -
            (p[i+1][j]+p[i-1][j])/(dx*dx) + (p[i][j+1]+p[i][j-1])/(dy*dy) );
        tmp = (p[i+1][j]-2.*p[i][j]+p[i+1][j])/(dx*dx)
            + (p[i][j+1]-2.*p[i][j]+p[i][j-1])/(dy*dy) - d[i][j];
        res += tmp*tmp;
      }
    }
    res = sqrt(res);
    printf("%d Iterationen, Residuum: %f\n",k,res);
    k++;
  }

  return 0;
}