Optimizing Vector Dot Product in CUDA: Exploring Shared Memory and Reduction Techniques
Introduction
In the realm of parallel computing, CUDA (Compute Unified Device Architecture) provides a powerful platform to harness the full potential of GPUs. This blog post dives into the intricacies of implementing a vector dot product using CUDA, focusing on the efficient use of shared memory, handling race conditions, and applying reduction techniques to optimize performance.
Understanding the Vector Dot Product
The dot product of two vectors is a fundamental operation in many scientific and engineering applications. For vectors a\mathbf{a}a and b\mathbf{b}b of size nnn, the dot product is calculated as:
dot_product= ∑i=(0 to n−1)a[i]×b[i]
In this post, we’ll explore how to perform this operation on a GPU, leveraging CUDA for parallel computation.
1. Code Overview
The primary goal of this CUDA program is to compute the dot product of two vectors a\mathbf{a}a and b\mathbf{b}b using GPU parallelism. The result of the dot product is stored in an array c\mathbf{c}c, where each element corresponds to a partial dot product computed by a block of threads.
2. Code Breakdown
Here’s a step-by-step explanation of the code:
#include<iostream>
#define n 10
using namespace std;
const int threadperblock = 256;
- Include and Definitions: The program includes the
iostream
header for input and output operations. It defines a constantn
for the size of the vectors andthreadperblock
to set the number of threads per block.
__global__ void dot_product(int *a, int *b, int *c) {
__shared__ float cache[threadperblock];
- Kernel Function:
__global__
indicates that this function will run on the GPU. Thedot_product
kernel takes three pointers:a
,b
(input vectors), andc
(output array for results). The__shared__
keyword declares an arraycache
in shared memory, accessible by all threads within a block.
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int cacheIndex = threadIdx.x;
float temp = 0; while (tid < n) {
temp += a[tid] * b[tid];
tid += blockDim.x * gridDim.x;
}
- Thread Indexing and Computation: Each thread calculates its unique index
tid
based on itsthreadIdx
andblockIdx
. Thetemp
variable accumulates the dot product result for each thread. Thewhile
loop ensures that threads handle elements beyond their block's initial range by incrementingtid
in steps of the total number of threads across all blocks.
cache[cacheIndex] = temp;
__syncthreads();
- Storing and Synchronizing: Each thread stores its computed result in the shared
cache
array. The__syncthreads()
function is a barrier that ensures all threads have finished writing tocache
before any thread proceeds.
int i = blockDim.x / 2;
while (i != 0) {
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
i /= 2;
}
- Reduction Technique: The reduction loop combines partial results from different threads within a block. Initially, each thread has its own
temp
value incache
. Threads iteratively add values fromcache
to perform a reduction, effectively summing up contributions from all threads to produce a final result per block.
if (cacheIndex == 0) {
c[blockIdx.x] = cache[0];
}
}
- Storing Final Result: The first thread in each block (
cacheIndex == 0
) writes the final reduced result of the block to the global output arrayc
.
3. Host Code
int main() {
int a[n], b[n], c[n], cpu_stored, stored_dot_val;
int *dev_a, *dev_b, *dev_c;
for (int i = 0; i < n; i++) {
a[i] = i + 2;
b[i] = i + 1;
} cpu_stored = 0;
for (int i = 0; i < n; i++) {
cpu_stored += a[i] * b[i];
} cout << "CPU Stored Value: " << cpu_stored << endl;
- Initialization: The host code initializes vectors
a
andb
with sample values and calculates the dot product on the CPU for comparison.
cudaMalloc((void **) &dev_a, n * sizeof(int));
cudaMalloc((void **) &dev_b, n * sizeof(int));
cudaMalloc((void **) &dev_c, n * sizeof(int));
cudaMemcpy(dev_a, a, n * sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(dev_b, b, n * sizeof(int), cudaMemcpyHostToDevice);
- Memory Allocation and Data Transfer: CUDA functions allocate memory on the GPU for vectors and copy the data from host to device.
int blockpergrid = (n + threadperblock - 1) / threadperblock;
dot_product<<<blockpergrid, threadperblock>>>(dev_a, dev_b, dev_c); cudaMemcpy(&c, dev_c, blockpergrid * sizeof(int), cudaMemcpyDeviceToHost);
- Kernel Launch and Result Transfer: The
dot_product
kernel is launched with a calculated number of blocks per grid. After execution, results are copied back from the GPU to the host.
stored_dot_val = 0;
for (int i = 0; i < blockpergrid; i++) {
stored_dot_val += c[i];
}
cout << "Stored Value is: " << stored_dot_val << endl; cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_c); return 0;
}
- Result Aggregation and Cleanup: The final dot product result is obtained by summing the partial results from each block. Memory allocated on the GPU is freed, and the program ends.
Conclusion
In this blog post, we’ve explored an efficient CUDA implementation for computing the dot product of vectors, emphasizing the use of shared memory and the reduction technique to enhance performance. By leveraging shared memory, we minimized global memory accesses, which significantly improves computational speed. The reduction technique, applied within each block, ensures that partial results are combined effectively, making the algorithm more efficient.
Handling race conditions and using synchronization barriers like __syncthreads()
is crucial for ensuring the correctness of parallel computations. This approach allows multiple threads to safely collaborate and achieve accurate results.
For a complete view of the CUDA code and further exploration, visit my GitHub repository. There, you can find the full implementation and experiment with the code to better understand the concepts discussed.
Happy coding, and feel free to reach out with any questions or feedback!