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


double random_number(double *rmin, double *rmax){
  return (*rmin) + (float)rand()/(float)RAND_MAX*((*rmax)-(*rmin));
}


int main(int argc, char **argv){
  double **A;
  double *Af;
  int i,j,k;
  double r_min,r_max;
  int n_ges;

  int myrank,nproz;
  MPI_Status status;

  FILE *fp;

  n_ges = 500;
  r_min = 0.0;
  r_max = 1.0;

  srand(time(0));

  A = (double**)malloc(n_ges*sizeof(double));

  Af = (double*)malloc(n_ges*n_ges*sizeof(double));

  for (i=0; i<n_ges; i++){
    A[i] = &Af[i*n_ges];
  }

  MPI_Init(&argc,&argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
  MPI_Comm_size(MPI_COMM_WORLD, &nproz);

  if ( myrank==0 ){
    fp = fopen("original_matrix.dat","w");
    for(i=0; i<n_ges; i++){
      for(j=0; j<n_ges; j++){
        A[i][j] = random_number(&r_min,&r_max);
        while ( fabs(A[i][j])<1.0e-9 && i==j ) /* main diagonal elements */
          A[i][j] = random_number(&r_min,&r_max); /* must be other than zero */
        fprintf(fp,"%6.3f",A[i][j]);
      }
      fprintf(fp,"\n");
    }
    fclose(fp);
  }


  for(k=0; k<n_ges; k++){
    if ( myrank==0 ){
      MPI_Send(A[k],n_ges,MPI_DOUBLE,k%nproz,k,MPI_COMM_WORLD);
    }
    if ( myrank==k%nproz ){
      MPI_Recv(A[k],n_ges,MPI_DOUBLE,0,k,MPI_COMM_WORLD,&status);
    }
  }

  MPI_Barrier(MPI_COMM_WORLD);

  for(k=0; k<n_ges-1; k++){
    MPI_Bcast(&A[k][k],n_ges-k,MPI_DOUBLE,k%nproz,MPI_COMM_WORLD);
    for(i=k+1; i<n_ges; i++){
      if (myrank==i%nproz){
        A[i][k] = A[i][k]/A[k][k];
        for(j=k+1; j<n_ges; j++){
          A[i][j] = A[i][j] - A[i][k]*A[k][j];
        }
      }
    }
  }

  for(k=0; k<n_ges; k++){
    if ( myrank==k%nproz ){
      MPI_Send(A[k],n_ges,MPI_DOUBLE,0,k,MPI_COMM_WORLD);
    }
    if ( myrank==0 ){
      MPI_Recv(A[k],n_ges,MPI_DOUBLE,k%nproz,k,MPI_COMM_WORLD,&status);
    }
  }

  if ( myrank==0 ){
    fp = fopen("calculated_matrix.dat","w");
    for(i=0; i<n_ges; i++){
      for(j=0; j<n_ges; j++){
        tmp = 0.;
        k = 0;
        while( k<=i && k<=j ){
          if ( k==i ) tmp += A[i][j];
          else tmp += A[i][k]*A[k][j];
          k++;
        }
        fprintf(fp,"%6.3f",tmp);
      }
      fprintf(fp,"\n");
    }
    fclose(fp);
  }

  MPI_Finalize();
  return 0;
}