/*
nvcc MM-GPU-GM.cu -o MM-GPU-GM.out
./MM-GPU-GM.out
*/

#include "cuda.h"
#include <stdio.h>
#include <sys/time.h>
#include <stdlib.h>
void MatrixMultiplication(float *, float *, float *, int);
void verify(float *P, int size){
int i,j;
float sum;
for(i=0;i<size;i++)
	for(j=0;j<size;j++){
		sum += P[i*size+j];
		}
printf("Verification Sum: %f\n",sum);
}

main(void){

    
    
    const int Width = 512;

    int size = Width * Width * sizeof(float);
    float *M, *N, *P;

    // allocate memory on the CPU
    M = (float*)malloc(size);
    N = (float*)malloc(size);
    P = (float*)malloc(size);

    // initialize the matrices
    for (int y=0; y<Width; y++) {
		for (int x=0; x<Width; x++){
			M[y*Width + x] = x + y*Width;
			N[y*Width + x] = x + y*Width; 
		}
    }
	
    MatrixMultiplication(M, N, P, Width);
	 verify(P, Width);
    // free the memory allocated on the CPU
    free( M );
    free( N );
    free( P );

    return 0;
}

__global__ void Kernel(float *Md, float *Nd, float *Pd, int Width) {

  // Calculate the column index of the Pd element, denote by x
  int x = threadIdx.x + blockIdx.x * blockDim.x; 
  // Calculate the row index of the Pd element, denote by y
  int y = threadIdx.y + blockIdx.y * blockDim.y; 

  float Pvalue = 0;
  // each thread computes one element of the output matrix Pd.      
  for (int k = 0; k < Width; ++k) {
    Pvalue += Md[y*Width + k] * Nd[k*Width + x];
  }
  // write back to the global memory
  Pd[y*Width + x] = Pvalue;
 
}

void MatrixMultiplication(float *M, float *N, float *P, int Width) {

    int size = Width * Width * sizeof(float);
    float *Md, *Nd, *Pd;
    

    // allocate memory on the GPU
     cudaMalloc((void**)&Md, size);
     cudaMalloc((void**)&Nd, size);
     cudaMalloc((void**)&Pd, size);
	  
	  // capture the start time
	struct timeval timer_start, timer_end;
	gettimeofday(&timer_start, NULL);
    // transfer M and N to device memory
     cudaMemcpy(Md, M, size, cudaMemcpyHostToDevice);
     cudaMemcpy(Nd, N, size, cudaMemcpyHostToDevice);

    // kernel invocation code
    dim3 dimBlock(16, 16);
    dim3 dimGrid(Width/16, Width/16);
	 	
    Kernel<<<dimGrid, dimBlock>>>( Md, Nd, Pd, Width);
	
     cudaThreadSynchronize();
	 // transfer P from device  
	 cudaMemcpy(P,Pd,size,cudaMemcpyDeviceToHost);
	cudaThreadSynchronize();		
	 gettimeofday(&timer_end, NULL);   
	 double timer_spent = timer_end.tv_sec - timer_start.tv_sec +
         (timer_end.tv_usec - timer_start.tv_usec) / 1000000.0;
	 printf("Time spent: %.6f\n", timer_spent);
    // free the memory allocated on the GPU
     cudaFree(Md);
     cudaFree(Nd);
     cudaFree(Pd);

}