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

#define N 128

double get_dtime(void){
  struct timeval tv;
  gettimeofday(&tv, NULL);
  return ((double)(tv.tv_sec) + (double)(tv.tv_usec) * 0.001 * 0.001);
}

int main(int argc,char **argv){
  double A[N][N],B[N][N],C[N][N],C1[N/4][N],C2[N/4][N],C3[N/4][N];
  double d1,d2;
  int i,j,k;

  int procs,rank;
  MPI_Status stat;
  
  /*MPI初期化*/
  MPI_Init(&argc,&argv);
  /*全体のプロセス数所得*/
  MPI_Comm_size(MPI_COMM_WORLD,&procs);
  /*自分のランク番号取得*/
  MPI_Comm_rank(MPI_COMM_WORLD,&rank);

  if(rank==0){
    
    /*行列の初期化*/
    for(i = 0; i < N; i++) {
      for(j = 0; j < N; j++){
	if(i != j){
	  A[i][j] = 0.0;
	  B[i][j] = 0.0;
	}else{
	  A[i][j] = i+1;
	  B[i][j] = 1;
	} 
	C[i][j]=0.0;
	if(i<N/4){
	  C1[i][j] = 0.0;
	  C2[i][j] = 0.0;
	  C3[i][j] = 0.0;
	}
      }
    }
      
    d1 = get_dtime();
    /*データ送信 to P1*/
    MPI_Send(&A[N/4][0],N*N/4,MPI_DOUBLE,1,0,MPI_COMM_WORLD);
    MPI_Send(&B,N*N,MPI_DOUBLE,1,1,MPI_COMM_WORLD);
    /*データ送信 to P2*/
    MPI_Send(&A[2*(N/4)][0],N*N/4,MPI_DOUBLE,2,0,MPI_COMM_WORLD);
    MPI_Send(&B,N*N,MPI_DOUBLE,2,1,MPI_COMM_WORLD);
    /*データ送信 to P3*/
    MPI_Send(&A[3*(N/4)][0],N*N/4,MPI_DOUBLE,3,0,MPI_COMM_WORLD);
    MPI_Send(&B,N*N,MPI_DOUBLE,3,1,MPI_COMM_WORLD);

    /*行列積の計算(1/4)*/
    for(i=0;i<N/4;i++){
      for(j=0;j<N;j++){
	for(k=0;k<N;k++){
	  C[i][j] += A[i][k] * B[k][j];
	}
      }
    }
    
    /*他部分の計算結果受信*/
    MPI_Recv(&C[N/4][0],N*N/4,MPI_DOUBLE,1,2,MPI_COMM_WORLD,&stat);
    MPI_Recv(&C[2*(N/4)][0],N*N/4,MPI_DOUBLE,2,2,MPI_COMM_WORLD,&stat);
    MPI_Recv(&C[3*(N/4)][0],N*N/4,MPI_DOUBLE,3,2,MPI_COMM_WORLD,&stat);

    }
    else if(rank ==1){
      /*行列データ受信*/
      MPI_Recv(&A[N/4][0],N*N/4,MPI_DOUBLE,0,0,MPI_COMM_WORLD,&stat);
      MPI_Recv(&B,N*N,MPI_DOUBLE,0,1,MPI_COMM_WORLD,&stat);

      /*行列積の計算(2/4)*/
      for(i=N/4;i<2*(N/4);i++){
	for(j=0;j<N;j++){
	  for(k=0;k<N;k++){
	    C1[i][j] += A[i][k] * B[k][j];
	  }
	}
      }
      
      /*計算結果の送信*/
      MPI_Send(C1,N*N/4,MPI_DOUBLE,0,2,MPI_COMM_WORLD);
    }
    else if(rank ==2){
      /*行列データ受信*/
      MPI_Recv(&A[2*(N/4)][0],N*N/4,MPI_DOUBLE,0,0,MPI_COMM_WORLD,&stat);
      MPI_Recv(&B,N*N,MPI_DOUBLE,0,1,MPI_COMM_WORLD,&stat);

      /*行列積の計算(3/4)*/
      for(i=N/2;i<3*(N/4);i++){
	for(j=0;j<N;j++){
          for(k=0;k<N;k++){
            C2[i][j] += A[i][k] * B[k][j];
          }
	}
      }

      /*計算結果の送信*/
      MPI_Send(C2,N*N/4,MPI_DOUBLE,0,2,MPI_COMM_WORLD);
    }
    else if(rank ==3){
      /*行列データ受信*/
      MPI_Recv(&A[3*(N/4)][0],N*N/4,MPI_DOUBLE,0,0,MPI_COMM_WORLD,&stat);
      MPI_Recv(&B,N*N,MPI_DOUBLE,0,1,MPI_COMM_WORLD,&stat);

      /*行列積の計算(4/4)*/
      for(i=3*(N/4);i<N;i++){
	for(j=0;j<N;j++){
          for(k=0;k<N;k++){
            C3[i][j] += A[i][k] * B[k][j];
          }
	}
      }

      /*計算結果の送信*/
      MPI_Send(C3,N*N/4,MPI_DOUBLE,0,2,MPI_COMM_WORLD);
    }
  
  MPI_Finalize();

  d2 = get_dtime();
  printf("time = %f[sec]\n",d2-d1);

  return 0;
}