#include<stdio.h>
#include<time.h>
#include<sys/time.h>
#include<sys/resource.h>

#define N 1000				//行列のサイズ


double A[N][N], B[N][N];
double sec(void);

int main(void){
  int i, j, k;
  double s, sum;
  double start1, end1, start2, end2;
  
  printf("N = %d\n", N);
  
  
  /*行列A・Bの初期化*/
  for(i = 0; i < N;i++){  //行
    for(j = 0; j < N;j++){//列
      B[i][j] = 0.0;
      if(i == j){
	A[i][j] = 1.0;//対角成分
      }else{
	A[i][j] = 0.00000001;//それ以外の成分
      }
    }
  }
  
  //行列計算（シングル）
  start1 = sec();			//開始時間を取得
  
  for(i = 0; i < N;i++){
    for(j = 0; j < N;j++){
      for(k = 0; k < N; k++){
	B[i][j] += A[i][k] * A[k][j];
      }
    }
  }
  
  end1 = sec();				//終了時間を取得
  
  for(sum = 0.0, i = 0; i < N; i++){
    for(j = 0; j < N; j++)
      sum += B[i][j];
  }

  printf("single\n");
  printf("time    : %e\n", (end1 - start1));
  printf("time/N^3: %e\n", (end1 - start1)/N/N/N);
  printf("sum     : %f\n\n", sum);
  
  
  /*行列Bの初期化*/
  for(i = 0; i < N;i++){  //行
    for(j = 0; j < N;j++){//列
      B[i][j] = 0.0;
    }
  }
  
  //行列計算（並列化）
  start2 = sec();
#pragma omp parallel for
  for(i = 0; i < N;i++){
#pragma omp parallel for
    for(j = 0; j < N;j++){
      s = 0.0;
      
#pragma omp parallel for reduction(+:s)
      for(k = 0; k < N; k++)s += A[i][k] * A[k][j];
      
      B[i][j] = s;
    }
  }
  
  end2 = sec();
  
  for(sum = 0.0, i = 0; i < N; i++){
    for(j = 0; j < N; j++)
      sum += B[i][j];
  }
  printf("parallel\n");
  printf("time    : %e\n", (end2 - start2));
  printf("time/N^3: %e\n", (end2 - start2)/N/N/N);
  printf("sum     : %f\n\n", sum);

  printf("parallel / single : %e\n", (end2 - start2)/(end1 - start1)); 

  return 0;
}

//現在時間の取得
double sec()
{
    struct timeval tv;
    gettimeofday(&tv, NULL);
    return tv.tv_sec + (double)tv.tv_usec*1e-6;
}