On this post I’m going to solve a common problem on CUDA, and that is a Reduction problem.
But using the algorithms posted by Mark Harris “Optimizing Parallel Reduction in CUDA“
We know on GPU is not very effective perform a reduction operation, and trying to make the most optimal code can be very difficult and not very intuitive.
The post by Mark Harris only gives the algorithm but is not quite clear on how to use it, I’m going to do that.
Let us solve the following integral
this is not a calculus class, but in the old days, this was solved liked this
the sum of all the areas of the rectangles gives the approximate solution to the integral and is not other than PI . And yes you are indeed right, this is a reduction problem.
So… how are we going to tackle this one?
first we need to calculate all the rectangle areas. Since one rectangle does not depend on any other rectangle, is perfect for our GPU to solve this one.
__global__ void calcRectangles(float * rectangles,int n){
int i = blockIdx.x * blockDim.x + threadIdx.x;
float x;
if (i >= 0 && i < n){
x = (i + 0.5) / n;
rectangles[i] = 4.0 / (1.0 + x*x);
}
}
where in the input
float * rectangles is the array to store all the areas.
int n the precision, bigger the n the more accurate the solution will be.
this kernel will give us an array of float values, this result will be out input to the reduce kernel.
Now let us prepare for the real deal
for each reduction I can only use 32 threads, this is because of my graphics card limitation, son in order to make large reduction I needed to use recursive kernel calls.
For example, to reduce an array of 2^21 elements
in the image we can see how an array of 2^21 elements is reduced to an array of 2^16 elements.
Now the same process with this new array, resulting in a 2^11 array, until it has 31 elements or less.
float calcReductionGeneral(int n, float * rectangles,int total , int type){
clock_t start; //vars to measure time
clock_t end;
float secondsTotal = 0;
if (n<32){ //If the lenght of the array is less than 32, we make the reduction in the host
start = clock();
float * resultadoGlobaL = (float*)malloc(sizeof(float)*n);
float pi;
for (int i=0;i<n;i++)
pi = pi + rectangles[i];
pi = pi / total;
end = clock();
secondsTotal = ((float)(end - start) / CLOCKS_PER_SEC);
printf ("pi = %f " ,pi);
free(resultadoGlobaL);
return secondsTotal;
}
else{ //reduction on GPU
float * resParcial = (float*)malloc(sizeof(float)*n/32); //partial solution, array of 32 floats
for (int i=0;i<n/32;i++){
/*i.e. if we have 64 rectangles we are going to make 2 kernel calls
first call with 0..31 elemets
second call with 32..63 elements*/
float * inputParcial = (float*)malloc(sizeof(float)*32); //partial input to the kernel allocation
float * input_d;
gpuErrchk(cudaMalloc((void**)&input_d, 32*sizeof(float))); //GPU allocation
copia(inputParcial,rectangles,i);//copy the relative 32 elements of the rectangles to the partial input
/* i.e if we have 64 elements in the first loop iteration we copy the 0..31 elements and in the second iteration from 32..63 elements
*/
cudaMemcpy(input_d, inputParcial, sizeof(float)*32, cudaMemcpyHostToDevice);
//We move the parcial input to the GPU
float * output;
gpuErrchk(cudaMalloc((void**)&output, 32*sizeof(float))); //allocating memory on GPU
float * resParcial132 = (float*)malloc(sizeof(float)*32); //allocating memory for parcial result
start = clock();//we measure time
switch (type){ //calling the kernel, using the type parameter
case 0: reduce0 <<<1,32,32 >>>(input_d, output); break;
case 1: reduce1 <<<1,32,32 >>>(input_d, output); break;
case 2: reduce2 <<<1,32,32 >>>(input_d, output); break;
default: break;
}
end = clock();
secondsTotal = secondsTotal + ((float)(end - start) / CLOCKS_PER_SEC);
gpuErrchk (cudaThreadSynchronize());
gpuErrchk (cudaMemcpy(resParcial132, output, sizeof(float)*32, cudaMemcpyDeviceToHost)); // we move the result to the host
//printf ("resParcial132[0] = %f" ,resParcial132[0]);
gpuErrchk (cudaThreadSynchronize());
resParcial[i]= resParcial132[0];
//clean up
free(inputParcial);
free(resParcial132);
cudaFree(input_d);
cudaFree(output);
}
return calcReductionGeneral((int)n/32,resParcial,total,type) + secondsTotal; //recursive call.
}
}
__global__ void reduce0(float *g_idata, float *g_odata) {
extern __shared__ float sdata[];
// each thread loads one element from global to shared mem
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
sdata[tid] = g_idata[i];
__syncthreads();
// do reduction in shared mem
for
(unsigned int s = 1; s < blockDim.x; s *= 2) {
if (tid % (2 * s) == 0) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
// write result for this block to global mem
if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}
__global__ void reduce1(float *g_idata, float *g_odata) {
extern __shared__ float sdata[];
// each thread loads one element from global to shared mem
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
sdata[tid] = g_idata[i];
__syncthreads();
for (unsigned int s=1; s < blockDim.x; s *= 2) {
int index = 2 * s * tid;
if (index < blockDim.x) { sdata[index] += sdata[index + s]; } __syncthreads(); } // write result for this block to global mem if (tid == 0) g_odata[blockIdx.x] = sdata[0]; }
__global__ void reduce2(float *g_idata, float *g_odata) { extern __shared__ float sdata[]; // each thread loads one element from global to shared mem unsigned int tid = threadIdx.x; unsigned int i = blockIdx.x*blockDim.x + threadIdx.x; sdata[tid] = g_idata[i]; __syncthreads(); for (unsigned int s=blockDim.x/2; s>0; s>>=1) {
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
// write result for this block to global mem
if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}
where reduce0, reduce1, and reduce2 are exactly the same algorithms used by Mark Harris.
Now, something to point out
pi = 3.141593 Time reduce0: 0.456775 seconds
pi = 3.141593 Time reduce1: 0.459540 seconds
pi = 3.141593 Time reduce2: 0.458745 seconds
the times of each algorithm for 4M elements is pretty much the same, where is this article “Optimizing Parallel Reduction in CUDA” the speedup from reduce0 to reduc1 is 2.33.
So… yeah.. I think they are selling a little bit of hot air in here, BUT, my implementation is very different from Nvidia’s, I try to use all the threads all the time and in their case they don’t, sticking to 128 threads because that is when the achieve peak performance and stuff like that.
The main structure of the program is the same, a tree based calling kernel function in a recursive manner and the results are very different still.
You can clone my entire code
git clone https://adri1988@bitbucket.org/adri1988/picudareduction.git
