 	
 	TODO check: m.getWidth() = v.getSize()
 */
-Vector productMatrixVectorGPU_shared(const Matrix &h_m, const Vector &h_v) {
+Vector productMatrixVectorGPU_shared(const Matrix &h_m, const Vector &h_v) throw (cudaError_t) {
 	Matrix d_m(h_m, GPU);
 	Vector d_v(h_v, GPU);
 	

	
	const dim3 num_threads(CHECKBOARD_BLOCK_MAX_HEIGHT, CHECKBOARD_BLOCK_MAX_WIDTH, 1);
	const dim3 num_blocks((d_m.getHeight() + CHECKBOARD_BLOCK_MAX_HEIGHT -1)/CHECKBOARD_BLOCK_MAX_HEIGHT, (d_m.getWidth() + CHECKBOARD_BLOCK_MAX_WIDTH -1)/CHECKBOARD_BLOCK_MAX_WIDTH, 1);
	productMatrixVectorGPU_shared_kernel<<<num_blocks, num_threads, CHECKBOARD_BLOCK_MAX_HEIGHT + CHECKBOARD_BLOCK_MAX_WIDTH>>>(d_m, d_v, d_r);
	
	Vector h_r(d_r, CPU);
	return h_r;
}

int main(int argc, char **argv) {
	srand(time(NULL));
	
	Matrix h_m(WIDTH, HEIGHT, CPU);
	for (unsigned int i(0) ; i != HEIGHT ; i++) {
		for (unsigned int j(0) ; j != WIDTH ; j++) {
			h_m[i][j] = rand() % 10;
		}
	}
	//h_m.print();
	
	Vector h_v(WIDTH, CPU);
	for (unsigned int j(0) ; j != WIDTH ; j++) {
		h_v[j] = rand() % 10;
	}
	//h_v.print();
	
	clock_t t_chrono;
	t_chrono = clock();
	for (int i(0) ; i != NUM_TESTS ; i++) {
		productMatrixVectorCPU(h_m, h_v);
	}
	t_chrono = clock() - t_chrono;
	std::cout << "Measured time for <CPU>: " << ((float) t_chrono)/CLOCKS_PER_SEC << "s" << std::endl;
	
	t_chrono = clock();
	for (int i(0) ; i != NUM_TESTS ; i++) {
		productMatrixVectorGPU_naive(h_m, h_v);
	}
	t_chrono = clock() - t_chrono;
	std::cout << "Measured time for <GPU naive>: " << ((float) t_chrono)/CLOCKS_PER_SEC << "s" << std::endl;
	
	t_chrono = clock();
	for (int i(0) ; i != NUM_TESTS ; i++) {
		productMatrixVectorGPU_shared(h_m, h_v);
	}
	t_chrono = clock() - t_chrono;
	std::cout << "Measured time for <GPU shared>: " << ((float) t_chrono)/CLOCKS_PER_SEC << "s" << std::endl;
	
	return 0;
}

	const dim3 num_threads(CHECKBOARD_BLOCK_MAX_HEIGHT, CHECKBOARD_BLOCK_MAX_WIDTH, 1);
	const dim3 num_blocks((d_m.getHeight() + CHECKBOARD_BLOCK_MAX_HEIGHT -1)/CHECKBOARD_BLOCK_MAX_HEIGHT, (d_m.getWidth() + CHECKBOARD_BLOCK_MAX_WIDTH -1)/CHECKBOARD_BLOCK_MAX_WIDTH, 1);
	productMatrixVectorGPU_naive_kernel<<<num_blocks, num_threads>>>(d_m, d_v, d_r);
	cudaThreadSynchronize(); // block until the device is finished
	
	// check for error
	cudaError_t error = cudaGetLastError();
  	if(error != cudaSuccess) {
    		std::cout << "CUDA error: " << cudaGetErrorString(error) << std::endl;
		throw error;
	}
	
	Vector h_r(d_r, CPU);
	return h_r;
}

__global__ void productMatrixVectorGPU_shared_kernel(const Matrix d_m, const Vector d_v, Vector d_r) {
	extern __shared__ float block_result[]; // need: blockDim.x + blockDim.y
	
	int i = blockIdx.x * blockDim.x + threadIdx.x;
	int j = blockIdx.y * blockDim.y + threadIdx.y;
	
	if (i >= d_m.getHeight() || j >= d_m.getWidth())
		return;
	
	__syncthreads();
		
	block_result[threadIdx.x] = 0;
	block_result[blockDim.x + threadIdx.y] = d_v[j];
	__syncthreads();
	
	atomicAdd(&block_result[threadIdx.x], d_m.get(i, j) * block_result[blockDim.x + threadIdx.y]);
	__syncthreads();
	
