#include <omp.h>
#include <stdio.h>
#include <stdlib.h>

#include <sys/time.h>

#define N 1024
#define P 2

double A[N*N],B[N*N];

int main(int argc,char *argv[])
{
	struct timeval s,e;
	int i,j,k,m;
	double sum=0.0;
	omp_set_num_threads(P);

	gettimeofday(&s,NULL);
	for(m=0;m<10;m++){
	#pragma omp parallel for private(j) /*A*/
	for(i=0;i<N;i++){
		for(j=0;j<N;j++){
			A[N*i+j]=0.0;
			B[N*i+j]=0.0;
			if(i==j) A[N*i+j]=1.0;
		}
	}
	
	#pragma omp parallel for private(j,k) /*B*/
	for(i=0;i<N;i++){
		for(j=0;j<N;j++){
			for(k=0;k<N;k++){
				B[N*i+j]=B[N*i+j]+A[N*i+k]*A[N*k+j];
			}
		}
	}
	
	#pragma omp parallel for private(j) reduction(+:sum) /*C*/
	for(i=0;i<N;i++){
		for(j=0;j<N;j++){
			sum+=B[N*i+j];
		}
	}
	}
	gettimeofday(&e,NULL);
	
	printf("sum = %lf\ntime = %lf[s]\n",sum,((e.tv_sec-s.tv_sec)+(e.tv_usec-s.tv_usec)*1.0E-6)/10.0);
}