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

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

void main(int argc, char *argv[]) {
   double **a , *d, *x, *x_old, s;
   double om = 1.7;      /* omega */
   int i, j, k, n;
   int my_rank, p;
   double eps=1.0e-10;    /* epsilon */
   double n_rhs, n_lhs, n_rhs_complete, n_lhs_complete;   
	/* norm right- and left hand side of inequality */ 

   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) ); 
   x_old = (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 */

   do {
      memcpy( x_old, x, n * sizeof(double) );
      n_lhs = 0.0 ; 
      n_rhs = 0.0 ; 
      for(i = 0; i < n; i++) 
         if( my_rank == (i % p) ) {
            n_rhs += x_old[i] * x_old[i];
         }
      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]);
            n_lhs += (x[i]-x_old[i]) * (x[i]-x_old[i]);
	 }
      }
      MPI_Allreduce( &n_lhs, &n_lhs_complete, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
      MPI_Allreduce( &n_rhs, &n_rhs_complete, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );

      n_lhs_complete = sqrt( n_lhs_complete );
      n_rhs_complete = sqrt( n_rhs_complete );
      if( my_rank == 0 ) {
         printf("n_lhs_complete = % 6.2e   n_rhs_complete = % 6.2e\n", n_lhs_complete, n_rhs_complete);
      }
   } while ( n_lhs_complete > eps * n_rhs_complete );

   for(i = 0; i < n ; i++) {
      if( my_rank == (i % p) ) {
        printf("x[%d] = % 6.2e \n",i,x[i]); fflush(stdout);
      }
   }

   MPI_Finalize();
}

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