Homec4science
Diffusion Project_multiproc (master) Browse

multiproc/A4/implementation_trivial.cu
11cbb29bbfe2master

Last Change
Raw File

implementation_trivial.cu
View Options

/*
============================================================================
Filename : algorithm.c
Author : Your name goes here
SCIPER : Your SCIPER number
============================================================================
*/
#include <iostream>
#include <iomanip>
#include <sys/time.h>
#include <cuda_runtime.h>
using namespace std;
#define SHARED_MEM_SIZE 98304 // 96 KB
#define WARP_SIZE 32
#define COL_THREADS 16
#define ROW_THREADS 16
// CPU Baseline
void array_process(double *input, double *output, int length, int iterations)
{
double *temp;
for(int n=0; n<(int) iterations; n++)
{
for(int i=1; i<length-1; i++)
{
for(int j=1; j<length-1; j++)
{
output[(i)*(length)+(j)] = (input[(i-1)*(length)+(j-1)] +
input[(i-1)*(length)+(j)] +
input[(i-1)*(length)+(j+1)] +
input[(i)*(length)+(j-1)] +
input[(i)*(length)+(j)] +
input[(i)*(length)+(j+1)] +
input[(i+1)*(length)+(j-1)] +
input[(i+1)*(length)+(j)] +
input[(i+1)*(length)+(j+1)] ) / 9;
}
}
output[(length/2-1)*length+(length/2-1)] = 1000;
output[(length/2)*length+(length/2-1)] = 1000;
output[(length/2-1)*length+(length/2)] = 1000;
output[(length/2)*length+(length/2)] = 1000;
temp = input;
input = output;
output = temp;
}
}
/*
* General idea:
* - Allocate double precision shared memory for each block
* - 16 threads per block, because want full occupancy and each thread takes 2 banks
* - sizeof(shared) = 128b, then a block should not exceed sizeof(shared) / sizeof(double) = 16 values
*
* Idea before column reduction:
* - Start a new kernel
* - Avoid conflicts, sum over columns: reduct of 3 times the original matrix
* - Store transposed result into output
* - Swap input and output
*
* Idea after column reduction:
* - Start a new kernel
* - Avoid conflicts, sum over rows with padding of 1
* - Store result into output
*/
__global__ void Heat(double *input, double *output, int length)
{
unsigned int x = (blockIdx.x * blockDim.x) + threadIdx.x;
unsigned int y = (blockIdx.y * blockDim.y) + threadIdx.y;
unsigned int tid = threadIdx.x;
// grid boundary
if (((x > 0) && (y > 0) && (x < length-1) && (y < length-1))) {
if (((x == length/2-1) || (x == length/2)) && ((y == length/2-1) || (y == length/2))) {
output[x*length + y] = 1000;
} else {
output[y * length + x] =
input[(y-1) * length + (x-1)] +
input[(y-1) * length + (x)] +
input[(y-1) * length + (x+1)] +
input[y * length + (x-1)] +
input[y * length + x] +
input[y * length + (x+1)] +
input[(y+1) * length + (x-1)] +
input[(y+1) * length + (x)] +
input[(y+1) * length + (x+1)];
}
}
}
// GPU Optimized function
void GPU_array_process(double *input, double *output, int length, int iterations)
{
//Cuda events for calculating elapsed time
cudaEvent_t cpy_H2D_start, cpy_H2D_end, comp_start, comp_end, cpy_D2H_start, cpy_D2H_end;
cudaEventCreate(&cpy_H2D_start);
cudaEventCreate(&cpy_H2D_end);
cudaEventCreate(&cpy_D2H_start);
cudaEventCreate(&cpy_D2H_end);
cudaEventCreate(&comp_start);
cudaEventCreate(&comp_end);
/* Preprocessing goes here */
int SIZE = length*length*sizeof(double);
double* gpu_input;
double* gpu_output;
double* temp; //used to swith pointers between gpu_input and gpu_output
cudaMalloc((void**)&gpu_input, SIZE);
cudaMalloc((void**)&gpu_output, SIZE);
cudaEventRecord(cpy_H2D_start);
/* Copying array from host to device goes here */
cudaMemcpy((void*)gpu_input, (void*)input, SIZE, cudaMemcpyHostToDevice);
cudaMemcpy((void*)gpu_output, (void*)output, SIZE, cudaMemcpyHostToDevice);
dim3 thrsPerBlock(ROW_THREADS, COL_THREADS);
dim3 nBlks((length-1)/ROW_THREADS + 1,(length-1)/COL_THREADS + 1); //if length == 0,...,16 -> 1 block, if length == 17,...,32 -> 2 blocks etc...
cudaEventRecord(cpy_H2D_end);
cudaEventSynchronize(cpy_H2D_end);
//Copy array from host to device
cudaEventRecord(comp_start);
/* GPU calculation goes here */
for(size_t i = 0; i<iterations; i++)
{
// start kernel and compute evolution
Heat<<<nBlks, thrsPerBlock, SHARED_MEM_SIZE>>>(gpu_input, gpu_output, length);
// wait until GPU has finished
cudaDeviceSynchronize();
temp = gpu_input;
gpu_input = gpu_output;
gpu_output = temp;
}
gpu_output = gpu_input; //unswap from the last iteration
cudaEventRecord(comp_end);
cudaEventSynchronize(comp_end);
cudaEventRecord(cpy_D2H_start);
/* Copying array from device to host goes here */
cudaMemcpy((void*)output, (void*)gpu_output, SIZE, cudaMemcpyDeviceToHost);
cudaEventRecord(cpy_D2H_end);
cudaEventSynchronize(cpy_D2H_end);
/* Postprocessing goes here */
cudaFree(gpu_input);
cudaFree(gpu_output);
float time;
cudaEventElapsedTime(&time, cpy_H2D_start, cpy_H2D_end);
cout<<"Host to Device MemCpy takes "<<setprecision(4)<<time/1000<<"s"<<endl;
cudaEventElapsedTime(&time, comp_start, comp_end);
cout<<"Computation takes "<<setprecision(4)<<time/1000<<"s"<<endl;
cudaEventElapsedTime(&time, cpy_D2H_start, cpy_D2H_end);
cout<<"Device to Host MemCpy takes "<<setprecision(4)<<time/1000<<"s"<<endl;
}
c4science · Help