CUDA Kernel for parallel reduction

Summing a large array is another common problem that benefits from parallelization on a GPU. The

basic idea behind parallel reduction is to split the array into smaller chunks and sum them in parallel.

For example, given an array of N elements, you can have N/2 threads each sum two elements,

reducing the problem size by half. This process repeats until there’s only one element left—the sum of

the entire array.

Summing Arrays: Parallel Reduction Example Consider an array of 16 random numbers. The array

is reduced by summing pairs of elements in parallel until a single sum remains. The process is as

follows:

`` Array = [8, 3, 5, 7, 2, 9, 1, 6, 4, 10, 12, 15, 11, 14, 13, 16] ``

In this example, the array is reduced from 16 elements down to 8, then to 4, then to 2, and finally to

1, which is the sum of the entire array. At each step, pairs of elements are summed in parallel.

The final sum of the array is 136.

In actual GPU threads, each column corresponds to a thread, and the table shows how the values

change at each step of the reduction. Empty cells indicate that the thread has completed its task in

that step and is no longer active.

>## Cuda Parallel Reuction Source Code

Code Breakdown with Detailed Explanation

Kernel Definition

`` __global__ void reduceSum(float *input, float *output, int N) {

* __global__: This function runs on the GPU but is called from CPU
* input: Pointer to the large array we want to sum (in GPU memory)
* output: Pointer where each block will store its partial sum
* N: Total number of elements in the input array

### 1. Shared Memory Declaration

* What is Shared Memory?
    * Ultra-fast memory shared by all threads in the same block
    * ~100x faster than global GPU memory
    * Limited size (typically 16-48KB per block)
    * Used for thread communication within a block
* Why use shared memory here?
    * Each thread needs to access other threads' data during reduction
    * Global memory accesses would be too slow
### 2. Thread Identification

* Example: If we have 1024 elements and 256 threads per block:
    * Block 0: i = 0 to 255
    * Block 1: i = 256 to 511
    * Block 2: i = 512 to 767
    * Block 3: i = 768 to 1023
### 3. Loading Data into Shared Memory

Thread:   0  1  2  3  4  5  6  7  8  9  10 11 12 13 14 15
sdata:   [8, 3, 5, 7, 2, 9, 1, 6, 4, 10,12,15,11,14,13,16]

for (int stride = 1; stride < blockDim.x; stride *= 2) {
   int index = 2 * stride * tid;
   if (index < blockDim.x) {
       if (index + stride < blockDim.x) {
           sdata[index] += sdata[index + stride];
       }
   }
   __syncthreads();
}

const int N = 16;
float h_input[] = {8, 3, 5, 7, 2, 9, 1, 6, 4, 10, 12, 15, 11, 14, 13, 16};
float h_output[1];  // To store the final sum

float *d_input, *d_output;
cudaMalloc(&d_input, N * sizeof(float));   // 16 floats on GPU
cudaMalloc(&d_output, sizeof(float));      // 1 float on GPU

int threads = 16;
int sharedMem = threads * sizeof(float);  // 16 × 4 bytes = 64 bytes

reduceSum<<<1, threads, sharedMem>>>(d_input, d_output, N);

GPU VRAM (Device Memory)     CPU RAM (Host Memory)
---------------------        ---------------------
d_output[0] = 136.0    →    h_output[0] = 136.0