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

/*******************************************************************/

void main(int argc, char *argv[]) {
   double **a , *d, *x, s;
   double om = 1.7;      /* omega */
   int i, j, k, n;
   int k_stop = 100; 
   int my_rank, p;

   n = 8;   /* dimension of the problem */

   MPI_Init(&argc, &argv);
   MPI_Comm_size(MPI_COMM_WORLD, &p);
   MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

   a = (double **)malloc( n * sizeof(double*) );
   d = (double *) malloc( n * sizeof(double) );   
   x = (double *) malloc( n * sizeof(double) ); 

   for(i=0; i < n; i++)
      a[i] = (double *)malloc( n * sizeof(double) );

   if(my_rank == 0) {                           /* initialization */
      /* srand( (unsigned)time( (long *)0 )); */
      for(i = 0; i < n ; i++) {
         for(j = 0; j < n ; j++)
            a[i][j] = rand() /2147483647.0;
         a[i][i] += 11;  
         d[i] = rand() /2147483647.0;   
      }
   }

   for(i = 0; i < n ; i++)
      MPI_Bcast( a[i], n, MPI_DOUBLE, 0, MPI_COMM_WORLD ); 
   MPI_Bcast( d, n, MPI_DOUBLE, 0, MPI_COMM_WORLD );

   for(i = 0; i < n ; i++) 
      x[i] = 0.0;           /* start vector */

   for(k = 0; k <= k_stop ; k++) {
      for(i = 0; i < n; i++) {    
         s = 0.0;
         for(j = 0; j < n; j++) {
            if( my_rank == (j % p) ) {
               s += ( (i!=j) ? (-a[i][j]) : ((1-om)/om * a[i][j]) ) * x[j];
            } 
	 }
         MPI_Reduce( &s, &x[i], 1, MPI_DOUBLE, MPI_SUM, i % p, MPI_COMM_WORLD );
         if( my_rank == (i % p) ) {
            x[i] = om / a[i][i] * (x[i] + d[i]);
         }
      }
   }
   for(i = 0; i < n ; i++) {
      if( my_rank == (i % p) ) {
        printf("x[%d] = % 6.2e \n",i,x[i]); fflush(stdout);
      }
   }

   MPI_Finalize();
}

/*******************************************************************/