// Compute scalar product using GPU
# include "/home/Tit5/lecture27/gpu-lecture/book.h"
# define imin(a,b) (a<b?a:b)

const int N = 400000 * 1024; // size of input arrays
const int threadsPerBlock = 256;
const int blocksPerGrid = imin( 32, (N+threadsPerBlock-1) / threadsPerBlock ); // keep #blocks as small as possible



__global__ void dot( float *a, float *b, float *c ) {
	// __shared__ to make this variable resident in shared memory
        // So that a copy of cache will be created in each block
	// each copy of the cake can be modified only from the threads in the block
	__shared__ float cache[threadsPerBlock];
	int tid = threadIdx.x + blockIdx.x * blockDim.x;
	int cacheIndex = threadIdx.x;
	float temp = 0;

	while (tid < N) {
		temp += a[tid] * b[tid];
		tid += blockDim.x * gridDim.x;
	}

	// set the cache values
	cache[cacheIndex] = temp;
	// synchronize threads in this block
	// Wait that each thread finished its job

        //__syncthreads();

	// for reductions, threadsPerBlock must be a power of 2
	// because of the following code
	int i = blockDim.x/2;
	while (i != 0) {
		if (cacheIndex < i)
			cache[cacheIndex] += cache[cacheIndex + i];
	    //__syncthreads();
	    i /= 2;
	}
	
	// Copy the sum of the cache values in the output c
	// if statament is used to use only one thread per block to avoid
	// unnecessarily large amount of memory traffic 
	if (cacheIndex == 0)
		c[blockIdx.x] = cache[0];
	

	// Summary
	// 1) For each block you created a cache
	// 2) You summ all the element of the chache and stored in cache[0]
	// 3) you copy all those values in the output array c	
}


int main( void ) {
	float *a, *b, c, *partial_c;
	float *dev_a, *dev_b, *dev_partial_c;

	// allocate memory on the CPU side
	a = (float*)malloc( N*sizeof(float) );
	b = (float*)malloc( N*sizeof(float) );
	partial_c = (float*)malloc( blocksPerGrid*sizeof(float) );

	// allocate the memory on the GPU
	HANDLE_ERROR( cudaMalloc( (void**)&dev_a, N*sizeof(float) ) );
	HANDLE_ERROR( cudaMalloc( (void**)&dev_b, N*sizeof(float) ) );
	HANDLE_ERROR( cudaMalloc( (void**)&dev_partial_c, blocksPerGrid*sizeof(float) ) );

	// fill in the host memory with data
	for (int i=0; i<N; i++) {
		a[i] = i;
		b[i] = i*2;
	}


	// copy the arrays ‘a’ and ‘b’ to the GPU
	HANDLE_ERROR( cudaMemcpy( dev_a, a, N*sizeof(float), cudaMemcpyHostToDevice ) );
	HANDLE_ERROR( cudaMemcpy( dev_b, b, N*sizeof(float), cudaMemcpyHostToDevice ) );

	dot<<<blocksPerGrid,threadsPerBlock>>>( dev_a, dev_b, dev_partial_c );

	// copy the array 'c' back from the GPU to the CPU
	HANDLE_ERROR( cudaMemcpy( partial_c, dev_partial_c,
	blocksPerGrid*sizeof(float), cudaMemcpyDeviceToHost ) );
	// finish up on the CPU side
	c = 0;
	for (int i=0; i<blocksPerGrid; i++) {
		c += partial_c[i];
	}

	# define sum_squares(x) (x*(x+1)*(2*x+1)/6)
	printf( "GPU value %.12g Correct value %.12g?\n", c, 2 * sum_squares( (float)(N - 1) ) );
	
	// free memory on the GPU side
	cudaFree( dev_a );
	cudaFree( dev_b );
	cudaFree( dev_partial_c );

	// free memory on the CPU side
	free( a );
	free( b );
	free( partial_c );
}