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


#define DIMENSION 1500

void print_matrix(double **A, int d) {
    int i, j;

    for(i=0; i<d ; i++) {
       for(j=0; j<d; j++) {
          printf("% 6.2e  ",A[i][j]);
       }
       printf("\n");
    }
}


int main(int argc, char *argv[]) {
   double **a , **l;
   int   i, j, k, d;
   int my_rank, p;
   int dest;
   MPI_Status status;


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

   d = DIMENSION;

   a = (double **)malloc( d * sizeof(double*) );
   l = (double **)malloc( d * sizeof(double*) );   

   for(i=0; i<d; i++) {
      a[i] = (double *)malloc( d * sizeof(double) );
      l[i] = (double *)malloc( d * sizeof(double) );
   }

   for(i=0; i<d ; i++) {
      for(j=0; j<d ; j++) { 
         a[i][j] = 0.0;
         l[i][j] = 0.0;
      }
   }
   
   if(my_rank == 0) {
      srand( (unsigned)time( (long *)0 ));
      for(i=0; i<d ; i++) {
         for(j=0; j<d ; j++) { 
            a[i][j] = rand() /2147483647.0;   /* RAND_MAX = 2^31-1 = 2147483647 */
         }
      }
      for(i=0; i<d ; i++) {
         if(a[i][i]==0) { 
            printf("One diagonal element is zero! Exit!\n");
            MPI_Finalize();
            exit(1);
         }
      }
      for(i=0; i<d ; i++) {
         dest = i % p;
         if( dest != 0 ) {
            MPI_Send(a[i], d, MPI_DOUBLE, dest, i, MPI_COMM_WORLD);  
            for(j=0; j<d ; j++)
               a[i][j]= 0.0; 
         }
      }
   } else {
     for(i=0; i<d ; i++) {
        if(my_rank == (i % p) )
	   MPI_Recv(a[i], d, MPI_DOUBLE, 0, i, MPI_COMM_WORLD, &status);
      }
   }
/*
   printf("Process %d, Matrix a: \n",my_rank);
   print_matrix(a,d);
   printf("\n"); 
   fflush(stdout);  
*/
   for(k=0; k < d-1 ; k++) {
      MPI_Bcast(a[k], d, MPI_DOUBLE, k % p, MPI_COMM_WORLD);
      for(i=k+1; i < d ; i++) {
         if( my_rank == (i % p) ) {
            l[i][k] = a[i][k]/a[k][k];
            for(j=k; j < d; j++) {
               a[i][j] = a[i][j] - l[i][k]*a[k][j];
	    } 
         }
      }
   }

   /* in the end process ( d%p-1 + p ) % p   has the whole matrix U */
/*
   if( my_rank == ( d%p-1 + p ) % p ) {
      printf("Process %d, Matrix a: \n",my_rank);
      print_matrix(a,d);
      printf("\n");
      fflush(stdout);
   }
*/
   MPI_Barrier(MPI_COMM_WORLD);
   MPI_Finalize();

}