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

#define max(a,b) ((a)>(b)?(a):(b))

#define SIZE 10			/* Size of matrices */
#define tolerance 0.001		/* defines when to stop */

int N;			   	/* number of threads */

double grid[SIZE+1][SIZE+1], temp[SIZE+1][SIZE+1];

typedef struct {		/* used to pass args to threads */
  int from, to;
  double diff;
} Args;


void fill_matrix()
{
  int i, j;

  for (i=0; i<SIZE+1; i++) {
    grid[i][0] = temp[i][0] = 0.0;
    grid[i][SIZE] = temp[i][SIZE] = 0.0;
    grid[0][i] = temp[0][i] = 1.0;
    grid[SIZE][i] = temp[SIZE][i] = 0.0;
  }

  for (i=1; i<SIZE; i++)
    for (j=1; j<SIZE; j++)
      grid[i][j] = 0.0;
}


void print_matrix()
{
  int i, j = 0;
  printf("\n");
  for (i=0; i<SIZE+1; i++) {
    for (j=0; j<SIZE+1; j++)
      printf("% 5.3f", grid[i][j]);
    printf("\n");
  }
  printf("\n");
}


void* sor_compute (void* arg)
{
  Args *s = (Args *)arg;
  int i,j;

  for (i=s->from; i < s->to; i++) 
    for (j=1; j<SIZE; j++)
      temp[i][j]=0.25* (grid[i-1][j] + grid[i+1][j] + grid[i][j-1] + grid[i][j+1]);
  
}


void* sor_update (void* arg)
{
  Args *a = (Args *)arg;
  int i,j;

  a->diff = 0;
  for (i=a->from; i < a->to; i++) 
    for (j=1; j<SIZE; j++) {
      a->diff = max (a->diff, fabs(grid[i][j]-temp[i][j]));
      grid[i][j]=temp[i][j];
    }
}


int main(int argc, char *argv[])
{
  pthread_t *thread;
  int i;
  double d;
  Args *args;

  if (argc!=2) {
    printf("Usage: %s number_of_threads\n",argv[0]);
    exit(-1);
  }  

  N=atoi(argv[1]);
  fill_matrix();
  args = malloc(N*sizeof(Args));
  thread = malloc(N*sizeof(pthread_t));

  for (i=0; i<N; i++) {
    args[i].from = i * (SIZE-1)/N + 1;
    args[i].to= ((i+1) * (SIZE-1))/N +1;
  }
    

  /* 
   * this example also shows how to pass multiple arguments 
   * (both in and out) to the threads 
   */

  do {
    for (i=1; i<N; i++) {
      if (pthread_create (&thread[i], NULL, sor_compute, &args[i]) != 0 ) {
	perror("Can't create thread");
	exit(-1);
      }
    }
    sor_compute(&args[0]);

    for (i=1; i<N; i++) pthread_join (thread[i], NULL);

    for (i=1; i<N; i++) {
      if (pthread_create (&thread[i], NULL, sor_update, &args[i]) != 0 ) {
	perror("Can't create thread");
	exit(-1);
      }
    }
    sor_update(&args[0]);

    d = 0;
    for (i=0; i<N; i++) {
      pthread_join (thread[i], NULL);
      d = max(d, args[i].diff);
    }

  } while (d>tolerance);
  

  print_matrix();

  return 0;

}