// Solves Jacobi Computation
// 
// Given an array A, iteratively replace its
// elements with the average of their four nearest
// neighbors, until the largest change between two
// consecutive iterations is less than delta. 
//
// Fixed size version 
//
// Copyright by mpC team, 1999
//
//    http://www.ispras.ru/~cbr
//    mpc@ispras.ru
//
										       


#include<stdio.h>

#define TYPE double

#define N 10
#define DELTA 0.0001
#define PREC 0.01

typedef TYPE (*NxN) [N][N];
TYPE AA[N+2][N+2],T[N][N];
NxN An,As,Ae,Aw,A;
int reps;


void write_region( NxN r ) {
  int i,j;
  for(i=0;i<N;i++) {
    for(j=0;j<N;j++) {
      printf("%lf ",(*r)[i][j]);
    }
    printf("\n");
  }
  printf("\n\n");
}

void write_matrix(  ) {
  int i,j;
  for(i=0;i<N+2;i++) {
    for(j=0;j<N+2;j++) {
      printf("%lf ",AA[i][j]);
    }
    printf("\n");
  }
  printf("\n\n");
}




void jacobi() {
  TYPE err;
  reps=0;
  
  AA[]=0;
  AA[N+1][1:N]=1.0;

  do {
    reps++;
    
    T[]=0.25*( AA[0:N-1][1:N]+
               AA[2:N+1][1:N]+
               AA[1:N][0:N-1]+
               AA[1:N][2:N+1]
               );
    
    err=[+][+]((AA[1:N][1:N]-T[])*(AA[1:N][1:N]-T[]));
    AA[1:N][1:N]=T[];
  }
  while(err>DELTA);
}


void report() {
  printf("\nJacobi scheme used %d iterations.\n", reps);
}

main() {
  jacobi();
  report();
}