CUDA Inclusive SCAN Hillis-Steele

remember starting my first class in CUDA programming. I was really excited, so many new information about hardware, threads and some grids … I really wanted to understand whats happening out there in the GPU core, but I had lots of other work unfortunately. My first homework was to implement a SCAN operation on GPU using CUDA-C. I am laughing all the time, when I remember what my code looked like, it was a nightmare … And was everything but a good, efficient and usefull code. In fact the GPU CUDA implementations was about 10 times slower than the CPU version …

Recently, I have been studying at udacity and you know, I dont like homeworks … But I decided to even the score by implemnting a good SCAN algorithm :) For those who have never met SCAN, its an operation on one array. The trick is that each output element is a sum of the previous elements in the input array. (input = {0,1,2,3,4,5,6 …} –> output ={0,1,3,6,10,15,21 …}). There are two major versions of this algorithm : INCLUSIVE and EXCLUSIVE SCAN. Today, we will talk about the first one and how to implement it on the GPU using the Hillis-Steele algorithm. Let us see how it works:

At first, take a look at the number of steps. Its a log[2](N) dependence, therefore we need 3 steps to calculate SCAN of 8 elements. Now how it works, its easy. We begin with a variable ,,space” = 1. During each iteration,we multiply this variable by 2. Now this is a very usefull variable, beause if we have ,,N” threads,each thread can do the correct operation: if ( ,,threadIdx” < space) –> Just rewrite last result for next step. If ( ,,threadIdx” > space) –> add an element positioned at [threadIdx – space] to the current one. Continue to next step.

Now this is really very easy, but note, that we have to synchronize writing and reading from memory. This is a problem in many algorithms (FFT, Bitonic Sort, SCAN …) … As you know we can only synchronize withing a block, so we can calculate a SCAN of up to 1024 elements without any problems:

__global__ void Kernel(int *input, int *output)
{
int tid = threadIdx.x;
//USE SHARED MEMORY - COMON WE ARE EXPERIENCED PROGRAMMERS
__shared__ int Cache[1024];
Cache[tid] = input[tid];
__syncthreads();
int space = 1;
//BEGIN 
for(int i = 0;i<10;i++){
int temp = Cache[tid];
int neighbor = 0;
if((tid-space)>= 0){
neighbor = Cache[tid-space];
}
__syncthreads(); //AFTER LOADING
if(tid < space){
//DO NOTHING
}else{
Cache[tid] = temp + neighbor; 
}
space = space*2;
__syncthreads();
}
//REWRITE RESULTS TO MAIN MEMORY
output[tid] = Cache[tid];
}

__global__ void Kernel(int *input, int *output)

{

int tid = threadIdx.x;

//USE SHARED MEMORY - COMON WE ARE EXPERIENCED PROGRAMMERS

__shared__ int Cache[1024];

Cache[tid] = input[tid];

__syncthreads();

int space = 1;

//BEGIN

for(int i = 0;i<10;i++){

int temp = Cache[tid];

int neighbor = 0;

if((tid-space)>= 0){

neighbor = Cache[tid-space];

}

__syncthreads(); //AFTER LOADING

if(tid < space){

//DO NOTHING

}else{

Cache[tid] = temp + neighbor;

}

space = space*2;

__syncthreads();

}

//REWRITE RESULTS TO MAIN MEMORY

output[tid] = Cache[tid];

}

And of course,because accessing main memory takes too much time, we use shared memory instead to implement the problem more effectively :) …. There is still however one big problem: GPU Algorithms calculate effectively, when they have huge amounts of data to process … And for sure, 1024 elements is simply not enough and therefore is better to leave the calculation to the CPU. On the other hand i like to use the GPU :) Furthermore, I decided to implement this algorithm to be able to handle these huge amounts of data unlike in my previous post, where we have been limited by 65536 elements :) The trick is in using 2D grid, I mean, when we will have X and Y dimensions of the running blocks :

dim3 THREADS_PER_BLOCK(1024,1,1);
dim3 BLOCKS_PER_GRID;
if(N>65536){
BLOCKS_PER_GRID = dim3(64,N/65536,1);
}else{
BLOCKS_PER_GRID = dim3(N/1024,1,1);
}

dim3 THREADS_PER_BLOCK(1024,1,1);

dim3 BLOCKS_PER_GRID;

if(N>65536){

BLOCKS_PER_GRID = dim3(64,N/65536,1);

}else{

BLOCKS_PER_GRID = dim3(N/1024,1,1);

}

This is just a simple modification added to our host code, however we actually need to know how to synchronize the memory reading/writing. I am curious about this for a long time,but yet I havent found anything better than simply using 2(thats enough) or 3 memory arrays and calling multiple times a kernel with different arguments. This means something like ,,Put the synchronization part on the Host“ … It is a good method I believe and I have updated it from last time (FFT long arrays) because i do not need to copy the whole array back to input after the work is finished. Instead during each step, we simply replace input and output (Step 0: read from input,write to output. Step 1: read from output,write to input …):

__global__ void Kernel(int *input, int *output,int* results, int space, int step, int steps, bool Direction)
{
int tix = threadIdx.x + blockDim.x * blockIdx.x; 
int tiy = threadIdx.y + blockDim.y * blockIdx.y;
int tid = tix + tiy*gridDim.x*blockDim.x;
// space =      GAP BETWEEN NEIGHBORS
// direction =  COPY FROM INPUT TO OUTPUT OR FROM OUTPUT TO INPUT
// step =       CURRENT PROCESSING STEP
// steps =      TOTAL AMOUNT OF STEPS NEEDED TO CALCULATE SCAN
int res = 0;
if(Direction){
//OUTPUT -> INPUT
if(tid<space){
res  = output[tid]; //ONLY REWRITE TO CORRECT MEMORY ADDRESS
input[tid] = res;
}else{
res = output[tid] + output[tid-space];
input[tid] = res;
}
}else{
//INPUT -> OUTPUT
if(tid<space){
res = input[tid]; //ONLY REWRITE TO CORRECT MEMORY ADRESS
output[tid] = res;
}else{
res = input[tid] + input[tid-space];
output[tid] = res;
}
}
//THE FINAL STEP: WRITE RESULTS INTO CORRECT LOCATION ON GPU
if(step == (steps-1)){
results[tid] = res;
}
}

__global__ void Kernel(int *input, int *output,int* results, int space, int step, int steps, bool Direction)

{

int tix = threadIdx.x + blockDim.x * blockIdx.x;

int tiy = threadIdx.y + blockDim.y * blockIdx.y;

int tid = tix + tiy*gridDim.x*blockDim.x;

// space = GAP BETWEEN NEIGHBORS

// direction = COPY FROM INPUT TO OUTPUT OR FROM OUTPUT TO INPUT

// step = CURRENT PROCESSING STEP

// steps = TOTAL AMOUNT OF STEPS NEEDED TO CALCULATE SCAN

int res = 0;

if(Direction){

//OUTPUT -> INPUT

if(tid<space){

res = output[tid]; //ONLY REWRITE TO CORRECT MEMORY ADDRESS

input[tid] = res;

}else{

res = output[tid] + output[tid-space];

input[tid] = res;

}

}else{

//INPUT -> OUTPUT

if(tid<space){

res = input[tid]; //ONLY REWRITE TO CORRECT MEMORY ADRESS

output[tid] = res;

}else{

res = input[tid] + input[tid-space];

output[tid] = res;

}

//THE FINAL STEP: WRITE RESULTS INTO CORRECT LOCATION ON GPU

if(step == (steps-1)){

results[tid] = res;

}

We also need a minor modification to the host code as I have said already:

//LAUNCH KERNELS IN LOOP
int space = 1;
int steps = static_cast<int>(log2(static_cast<float>(N)));
for (size_t step = 0; step<steps; step++){
bool direction = (step % 2) !=0;  
Kernel << <BLOCKS_PER_GRID, THREADS_PER_BLOCK >> >(dev_input, dev_output, dev_results, space, direction, step, steps);
space = space * 2;
}

//LAUNCH KERNELS IN LOOP

int space = 1;

int steps = static_cast<int>(log2(static_cast<float>(N)));

for (size_t step = 0; step<steps; step++){

bool direction = (step % 2) !=0;

Kernel << <BLOCKS_PER_GRID, THREADS_PER_BLOCK >> >(dev_input, dev_output, dev_results, space, direction, step, steps);

space = space * 2;

}

The most important thing in implementing an algorithm is to always make a picture how it works and think about it for a while :) Thats really all. Have a good day and feel free to email me back at : [email protected] … you can alternatively join the Intro to parallel Computing at udacity or post some comments :)

EDIT: 18.4.2016

The old code is a mess and can be easily replaced with a modern version (See last link that implements the SCAN on large arrays).The Major Modifications includes:

No additional variables during BLOCK/THREAD Initialization
Automatic calculation of the number of steps
Kernel Calls with parameters rather than memcpy(4*int)
CPU Version of the same algorithm added for comparison
More user-Friendly For-loops, results displaying and program end
Malloc has been replaced with “new” operator