|| Jai Sri Gurudev||

Sri Adichunchanagiri Shikshana Trust(R)

SJB INSTITUTE OF TECHNOLOGY

BGS Health & Education City, Dr. Vishnuvardhan Road, Kengeri, Bengaluru–560060
An Autonomous Institute under Visvesvaraya Technological University, Belagavi
Affiliated to Visvesvaraya Technological University, Belagavi & Approved by AICTE, New Delhi, Certified by ISO 9001-2015
Accredited by NBA & NAAC, New Delhi with ‘A+’ Grade, Recognized by UGC, New Delhi with 2(f) and 12(B)

Notes On
Course Name: Parallel Computing

Course Code: B C S 7 0 2

Module – 5

By

Faculty Name: Dr Ranjith J , Prof. Kiran Kumar V , Prof. Veeresh K M

Semester: VIIth Semester

Department of Information Science & Engineering

Aca. Year ODD SEM /2025-26

Parallel Computing-BCS702 Module 5

GPU programming with CUDA

5.1 GPUs and GPGPU
The rise of realistic graphics in the late 1990s and 2000s led to the development of
powerful Graphics Processing Units (GPUs).
GPGPU Definition: General Purpose computing on Graphics Processing Units
(GPGPU) refers to using GPUs for tasks beyond graphics rendering.
GPU Programming: Originally, GPUs could only be programmed with graphics-
specific APIs like Direct3D and OpenGL. Early developers had to reformulate general
problems into graphics concepts like vertices and triangles.
CUDA and OpenCL: CUDA (NVIDIA) and OpenCL (cross-vendor) emerged to
simplify GPU programming for general-purpose tasks. CUDA Popularity: CUDA is
widely used due to its ease of use and support for NVIDIA GPUs. OpenCL
Compatibility: OpenCL works on multiple platforms but often requires more complex
setup.

5.2 GPU architectures

GPUs vs CPUs: GPUs consist of many simpler processors compared to a CPU’s
fewer but more powerful cores. SIMD Architecture: GPUs often use SIMD (Single
Instruction, Multiple Data) model for parallel processing. SIMD Execution: One
instruction is executed on multiple data points simultaneously.
Example Program: Consider a data path executing:
1. if (a[i] > 0)
2. x[i] = 1;
3. else
4. x[i] = 0;
Branching in SIMD: If the data paths differ in condition result, execution gets more
complex. When different SIMD paths need to execute different instructions,
performance can suffer.
SIMD vs SISD: SISD (Single Instruction, Single Data) contrasts SIMD by only
executing one instruction on one data point. GPUs use a variation called SIMT
(Single Instruction, Multiple Threads), combining aspects of SIMD and
multithreading. SIMD paths diverge (e.g., a[i] > 0 is true for some threads and false
for others), execution splits. Conditional branching reduces efficiency in SIMD due to
inactive threads during some instruction paths.
Control Unit Role: A control unit fetches a single instruction and broadcasts it to all
active data paths. GPUs are well-suited for tasks with massive data-level parallelism
but less so for heavily branched logic. Understanding GPU architecture is essential
for writing efficient CUDA programs.

1
Dept. of ISE
Parallel Computing-BCS702 Module 5

In the late 1990s, GPUs became powerful for rendering graphics in games and
animations. Developers began using GPUs for non-graphics applications, leading to
GPGPU (General Purpose computing on GPUs). GPUs were originally programmed
using APIs like Direct3D and OpenGL. These APIs were limited to graphics;
developers had to reformulate problems using graphics constructs (e.g., vertices,
textures).
Emergence of CUDA and OpenCL
CUDA (developed by NVIDIA) allows direct programming of GPUs using a C/C++-
like language. OpenCL is a cross-platform API supporting multiple devices including
GPUs and CPUs. CUDA is widely adopted due to its simplicity but only supports
NVIDIA GPUs. CPUs follow SISD (Single Instruction, Single Data), executing one
instruction on one data at a time. GPUs use SIMD (Single Instruction, Multiple Data),
applying the same instruction across many data elements. GPU hardware is
optimized for parallel data processing with many lightweight cores.
Example: Branching in SIMD
if (a[i] > 0)
x[i] = 1;
else
x[i] = 0;
 If a[i] > 0 varies across data, different execution paths must be taken.
 SIMD struggles with such divergence, causing inefficiencies.
 Shows how execution splits when different data paths require different instructions.
Paths with a[i] > 0 are executed separately from those with a[i] <= 0.
SIMT – Single Instruction, Multiple Threads
 Used in CUDA architecture.
 Threads in a group (warp) execute the same instruction concurrently unless control
divergence occurs.
 Threads within a warp follow different execution paths when conditions vary.
 Reduces efficiency as some threads must wait for others to complete.

GPU SM Architecture and Heterogeneous Computing

Streaming Multiprocessors (SMs)
 Modern NVIDIA GPUs contain multiple SMs.
 Each SM has many scalar processors (SPs) controlled by a control unit.
 SMs execute instructions in SIMT fashion.
SIMT vs SIMD
 SIMT is similar to SIMD but uses thread-based execution.

2
Dept. of ISE
Parallel Computing-BCS702 Module 5

 Threads are scheduled independently, but divergence still impacts performance.

Each SM has Shared Memory – fast memory accessible to threads in the same
block.
 GPUs also have Global Memory – accessible by all SMs but slower than shared
memory.
 Accessing shared memory is much faster than global memory.
 If some threads need to execute x[i] = 1 and others x[i] = 0, SMs handle these on
different schedules.
 Increases overhead when multiple divergent branches are handled across SMs.
CUDA Warps
 A warp is a group of typically 32 threads executed together on an SM.
 Threads within a warp execute instructions in lockstep unless control divergence
occurs.
GPU Memory Hierarchy
 Registers (fastest, private to threads).
 Shared Memory (fast, shared among threads in a block).
 Global Memory (large but slow, shared across all blocks).
 L2 Cache acts as a buffer between SMs and Global Memory.
 CPU = Host; GPU = Device.
 Programs are split between host (CPU code) and device (GPU kernels).
 CPU manages memory and launches kernels to run on the GPU.

5.3 Heterogeneous Computing

 Modern systems often include both CPU and GPU.
 Work is divided so that CPUs handle sequential logic, while GPUs handle parallel
computation.
 Developers must write code that orchestrates cooperation between both devices.
 Data is transferred from CPU main memory to GPU global memory before kernel
execution.
 After execution, results must be copied back to CPU memory.
Execution Diagram (Figure 5.1)
 Shows SM architecture: Each SM has multiple SPs and shared memory.
 Illustrates hierarchical structure of GPU computing units.
CPU-GPU Block Diagram (Figure 5.2)

3
Dept. of ISE
Parallel Computing-BCS702 Module 5

 Depicts connection between CPU (host) and GPU (device).

 Data flows through global memory and shared memory for computation.
 CUDA enables general-purpose programming on NVIDIA GPUs.
 Understanding memory hierarchy and SM execution model is crucial.
 Proper design avoids branch divergence and optimizes memory use for performance.

Heterogeneous Computing Overview

 Heterogeneous computing uses both CPUs and GPUs to execute a program.
 CPUs handle general-purpose sequential logic; GPUs execute data-parallel portions.

4
Dept. of ISE
Parallel Computing-BCS702 Module 5

 CUDA enables this model, allowing the GPU to handle parts of the program using its
many cores. GPUs use SIMT (Single Instruction, Multiple Threads), allowing
thousands of threads to execute the same instructions with separate data. This is
ideal for large data-parallel problems.
 CUDA allows users to write GPU-accelerated code in a C/C++-like syntax.
 Developers must manage device memory, define kernel functions (code that runs on
the GPU), and launch them from the CPU (host).
Introduction to CUDA API
 CUDA provides several libraries and APIs to facilitate GPU programming.
 Common CUDA components include:
o nvcc: NVIDIA’s CUDA compiler.
o CUDA runtime API: Provides functions for memory management and kernel
launch.
o CUDA driver API: Lower-level interface compared to the runtime API.
Example of CUDA Memory Allocation
c
CopyEdit
float *a = (float *)malloc(sizeof(float));
 This C-style memory allocation works on the host (CPU).
 On the device (GPU), you must use CUDA-specific memory allocation functions like
cudaMalloc().

5.4 CUDA HELLO

Writing a CUDA Program – “Hello, World”
The Source Code
 The first CUDA example prints a greeting from each GPU thread.
 It consists of:
A kernel function (marked with __global__) that runs on the GPU.
A main function that runs on the host, which launches the kernel.

Kernel Function
c
CopyEdit
__global__ void hello()

5
Dept. of ISE
Parallel Computing-BCS702 Module 5

{
printf("Hello from GPU thread!\n");
}
 __global__ indicates the function will execute on the GPU and can be called from the
host.
 printf() is supported within device code starting with CUDA 5.0+.
Launching a Kernel
c
CopyEdit
hello<<<1, 5>>>();
 This launches 1 block with 5 threads.
 The syntax <<<blocks, threads>>> controls the GPU execution configuration.
Compiling a CUDA Program
CUDA Compiler
 nvcc is the compiler used to compile CUDA programs.
 The command:
bash
CopyEdit
nvcc -o hello [Link]
 This compiles the CUDA source file [Link] and outputs an executable hello.
Running the Program
 Once compiled, run the binary like a standard program:
bash
CopyEdit
./hello
 Output should show 5 print statements, one from each thread.
 CUDA simplifies GPU programming by using familiar C/C++ syntax.
 Key components: kernel definition, memory management, and kernel launch
configuration.
 Proper compilation and execution setup is required to run CUDA code effectively.

6
Dept. of ISE
Parallel Computing-BCS702 Module 5

5.5 A Closer Look at CUDA Execution

1. Kernel Launch Review

7
Dept. of ISE
Parallel Computing-BCS702 Module 5

 When a CUDA kernel (e.g., hello<<<1,5>>>();) is launched:

o <<<1, 5>>> specifies launching 5 GPU threads in 1 block.
o Each thread runs the same kernel function, e.g., hello.
2. Thread Identification
 Each thread identifies itself using:
int tid = threadIdx.x;
 threadIdx.x is a built-in CUDA variable that gives the thread’s index within a block.
3. Single Program, Multiple Data (SPMD)
 CUDA follows the SPMD model where each thread runs the same code but operates
on different data.
4. Asynchronous Execution
 CUDA kernel launches are asynchronous:
o The host program continues immediately after the kernel launch.
o The kernel runs in the background unless synchronization is explicitly
required.
5. Host-Device Synchronization
 To ensure the host waits for the device to finish:
cudaDeviceSynchronize();
 This blocks the CPU until all GPU threads have completed.
6. Execution Flow Recap
 CUDA flow:
1. Kernel is launched.
2. GPU schedules and executes threads in parallel.
3. Threads identify themselves via threadIdx.
4. Optional: Synchronize host using cudaDeviceSynchronize().

5.6 Threads, Blocks, and Grids

Why <<<1,5>>>?
 The angle bracket syntax specifies:
o 1 block
o 5 threads per block
 Syntax: kernel<<<numBlocks, threadsPerBlock>>>();
8CUDA Hardware Layout

8
Dept. of ISE
Parallel Computing-BCS702 Module 5

 NVIDIA GPUs contain:

o Multiple SMs (Streaming Multiprocessors)
o Each SM executes many threads in groups called warps
CUDA Execution Hierarchy
 3 main components:
o Thread: Basic execution unit.
o Block: A group of threads (e.g., 256, 512, etc.)
o Grid: A collection of blocks.
Example
hello<<<2, 512>>>();
 Launches a grid with 2 blocks, each with 512 threads → 1024 total threads.
Thread Indexing
 Threads are identified by:
o threadIdx.x: Index within a block.
o blockIdx.x: Index of the block within the grid.
o Combine them for a global thread ID:
int globalID = blockIdx.x * blockDim.x + threadIdx.x;
Why Use More Than 1 Block?
 A single block has a limit on thread count (typically 1024 threads max).
 To scale beyond this, multiple blocks are needed.
Summary of CUDA Hierarchy

Level Description

Thread Executes kernel code

Block Group of threads, share local memory

Grid Group of blocks

Thread Cooperation
 Threads within a block can:
o Share shared memory
o Use __syncthreads() to synchronize
 Threads from different blocks cannot directly communicate.

9
Dept. of ISE
Parallel Computing-BCS702 Module 5

Note that all the blocks must have the same dimensions. More importantly, CUDA requires
that thread blocks be independent. So one thread block must be able to complete its
execution, regardless of the states of the other thread blocks: the thread blocks can be
executed sequentially in any order.

5.7 Nvidia compute capabilities and device architectures

1. What is Compute Capability?
 Compute capability is a version number (format: a.b) used by NVIDIA to describe
the hardware features supported by a GPU.
 Indicates compatibility with CUDA features and performance limits.

10
Dept. of ISE
Parallel Computing-BCS702 Module 5

2. Format and Range

 Valid values of a (major revision): 1, 2, 3, 5, 6, 7, 8.
 CUDA no longer supports GPUs with compute capability < 3.0.
 Minor revisions (the .b part) depend on GPU model and architecture.
3. Thread and Block Limits
 For compute capability ≥ 1.0:
o Max threads per block: 1024.
 For compute capability 2.b:
o Max threads per Streaming Multiprocessor (SM): 1536.
 For compute capability > 2.0:
o Max threads per SM: 2048.
 Max x, y, and z dimensions of thread blocks and grids are limited and vary by
compute capability.
4. Table 6.3 – GPU Architectures & Compute Capabilities

Architecture Compute Capability

Ampere 8.0

Turing 7.5

Volta 7.0

Pascal 6.b

Maxwell 5.b

Kepler 3.b

Fermi 2.b

Tesla 1.b

5. Key Notes
 Tesla is both a product name and an architecture family, which can be confusing.
 Some GPUs are part of NVIDIA SoCs (e.g., Jetson Nano), combining GPU + CPU on
a single chip.
 Compute capability affects kernel execution limits and memory usage.
 Knowing your GPU’s compute capability is essential for writing compatible and
efficient CUDA code.
 CUDA C++ programming guide provides detailed compute capability limits and
specs.

11
Dept. of ISE
Parallel Computing-BCS702 Module 5

5.8 Vector addition

GPUs and CUDA are highly efficient for data-parallel programs, where the same operation
is performed across many data elements. A classic example is vector addition, which is
embarrassingly parallel—each element can be computed independently.
The program involves:
 Creating three arrays: x, y, and z, each of size n.
 Initializing x and y on the host (CPU).
 Launching a CUDA kernel with at least n threads.
 Each thread computes:
 z[i] = x[i] + y[i];
 The arrays use float (32-bit) instead of double for better GPU efficiency.
 After kernel execution:
o The program checks the result,
o Frees allocated memory,
o And exits.
This example illustrates the simplicity and effectiveness of CUDA in handling parallel
numerical operations.

12
Dept. of ISE
Parallel Computing-BCS702 Module 5

5.8.1 The Kernel

1. Purpose
 Demonstrates an embarrassingly parallel CUDA program: adding two vectors.
 Each thread adds one pair of elements from arrays x and y to compute z[i] = x[i] +
y[i].
2. Thread Indexing in CUDA Kernels
 Each thread uses a global index to identify which data element it should process.
 Global index formula:

13
Dept. of ISE
Parallel Computing-BCS702 Module 5

 int rank = blockIdx.x * blockDim.x + threadIdx.x;

 This allows threads across blocks to process distinct vector elements.
3. Example Setup
 With 4 blocks and 5 threads per block:
o Global thread ranks range from 0 to 19 (see Table 6.4).

o
4. Handling Uneven Vector Sizes
 Always check bounds to avoid accessing beyond array size:
 if (rank < n)
 z[rank] = x[rank] + y[rank];
5. Serial vs. Parallel
 If CUDA is unavailable, a serial version loops through each index:
 for (int i = 0; i < n; i++)
 z[i] = x[i] + y[i];
 The CUDA version replaces this loop by assigning each iteration to a separate
thread.
6. Design Consideration
 This design simplifies the mapping of computations.
 Makes it easier to convert existing serial loops into GPU-parallel kernels.
7. Program 6.4 Summary
 Serial_vec_add() performs vector addition on the CPU.
 Highlights contrast between CPU loops and GPU thread-based execution.

14
Dept. of ISE
Parallel Computing-BCS702 Module 5

Here’s a concise summary of the provided section on Get_args, Allocate_vectors, and

unified memory in CUDA:

5.8.2 Get_args Function

 Purpose: Parses command-line arguments and returns:
o n: Number of elements in arrays
o blk_ct: Number of thread blocks
o th_per_blk: Threads per block
o i_g: Character flag to decide if user inputs arrays or uses random numbers
 Error handling:
o Prints usage and exits if arguments are incorrect
o Exits if n > total threads available
 Note: Written in standard C and runs on the host (CPU).

5.8.3 Allocate_vectors and Managed Memory

Allocates memory for four n-element float arrays: x, y, z, and cz.
o x, y, and z: Used by both host and device
o cz: Used only on the host (for verifying GPU results)
Allocation Methods
 cz is allocated using malloc() (host memory).
 x, y, and z are allocated using:
cudaMallocManaged(&ptr, size);
 This uses CUDA’s unified memory system, allowing both the CPU and GPU to
access the same memory.

Key Points on Unified Memory

 Simplifies CUDA development—no need to manually transfer memory between host
and device.
 Requirements:
1. Device must support compute capability ≥ 3.0
2. Requires a 64-bit host OS
 Limitations:
o Devices with compute capability < 6.0 can’t access unified memory from
both host and device at the same time.

15
Dept. of ISE
Parallel Computing-BCS702 Module 5

o May be slower than manually managed memory due to hidden data

transfers.
o Performance depends on when and how the system moves data between
host and device.

4. Performance Tip
 For better performance, developers may still prefer explicit memory transfers (e.g.,
using cudaMemcpy) to control timing and minimize overhead.

16
Dept. of ISE
Parallel Computing-BCS702 Module 5

17
Dept. of ISE
Parallel Computing-BCS702 Module 5

5.8.4 Other functions called from main

The Serial_vec_add function (Program6.4) just adds x and y on the host using a for loop. It
stores the result in the host array cz.
The Two_norm_diff function computes the ―distance‖ between the vector z computed by the
kernel and the vector cz computed by Serial_vec_add.

5.8.5 Explicit Memory Transfers in CUDA (No Unified

Memory)
 This approach modifies the vector addition program for systems without unified
memory.
 Host and device cannot share memory directly, so explicit memory allocation and
data transfer are required.

18
Dept. of ISE
Parallel Computing-BCS702 Module 5

2. Kernel
 Unchanged from the unified memory version.
 Receives x, y, z, and n as arguments.
 Each thread computes:
c
CopyEdit
if (my_elt < n)
z[my_elt] = x[my_elt] + y[my_elt];
3. Main Function Changes
 Structure is mostly the same, but memory handling differs:
o Pointers on the host are not valid on the device, and vice versa.
o Using a host pointer on the GPU can lead to errors or unintended behavior.
4. Separate Memory for Host and Device
 Host arrays: hx, hy, hz → Allocated with malloc()
 Device arrays: dx, dy, dz → Allocated with cudaMalloc()
 Memory is separated and must be manually transferred between host and device.
5. Allocation
 Memory allocation is handled in Allocate_vectors() (Program 6.10):
o Host arrays: malloc()
o Device arrays: cudaMalloc()
Without unified memory, developers must explicitly manage data transfers and allocations
using malloc and cudaMalloc. This provides more control but requires additional care to
ensure correctness and avoid memory access errors across host and device.

19
Dept. of ISE
Parallel Computing-BCS702 Module 5

20
Dept. of ISE
Parallel Computing-BCS702 Module 5

5.9 Returning Results from CUDA Kernels

1. Returning Values from Kernels
 CUDA kernels cannot return values like standard C functions.
 They must use pointers to memory shared between the host and device (or
explicitly copied back).

2. Pitfalls with Pointer-Based Returns

 Passing host pointers to kernels (e.g., via standard C pass-by-reference) often
causes errors because:
o The host and device have separate address spaces.
o Attempting to dereference a host pointer inside a device kernel will likely
crash or hang.

3. Correct Approaches
A. With Unified Memory
 Allocate result memory using cudaMallocManaged():

21
Dept. of ISE
Parallel Computing-BCS702 Module 5

 int *sum_p;
 cudaMallocManaged(&sum_p, sizeof(int));
 Kernel writes directly to sum_p, and the value becomes accessible from the host
after kernel execution.
B. Without Unified Memory
 Allocate:
o Device memory for result using cudaMalloc()
o Host memory for result using malloc()
 Copy result back using cudaMemcpy():
 cudaMemcpy(hsum_p, dsum_p, sizeof(int), cudaMemcpyDeviceToHost);

4. Alternate Option
 Instead of using a pointer, return the result in a global memory variable declared
with __device__.
 CUDA kernels don’t return values directly.
 Use shared or copied memory to communicate results from the device to the host.
 Unified memory simplifies the process, but manual copying is required without it.

5.10 CUDA trapezoidal rule I

5.10.1 The trapezoidal rule

1. Purpose
 Implement a CUDA version of the trapezoidal rule to approximate the area under a
curve between two points a and b.

2. Concept: The Trapezoidal Rule

 The interval [a,b][a, b] is divided into n subintervals.
 Each subinterval has width:
h=b−anh = \frac{b - a}{n}
 The left endpoint of the ithi^{th} subinterval:
xi=a+i⋅hx_i = a + i \cdot h
 The area of the ithi^{th} trapezoid:
h2[f(xi)+f(xi+1)]\frac{h}{2} [f(x_i) + f(x_{i+1})]

22
Dept. of ISE
Parallel Computing-BCS702 Module 5

3. Total Area Approximation

 Sum of all trapezoid areas:
h2[f(x0)+2f(x1)+2f(x2)+⋯+2f(xn−1)+f(xn)]\frac{h}{2} [f(x_0) + 2f(x_1) + 2f(x_2) + \cdots +
2f(x_{n-1}) + f(x_n)]
 Can be computed as:
h[12f(a)+f(x1)+f(x2)+⋯+f(xn−1)+12f(b)]h \left[ \frac{1}{2}f(a) + f(x_1) + f(x_2) + \cdots +
f(x_{n-1}) + \frac{1}{2}f(b) \right]

4. Serial Implementation (Program 6.11)

 Function Serial_trap() calculates the trapezoidal area on the CPU:
o Computes h
o Initializes result with 0.5⋅(f(a)+f(b))0.5 \cdot (f(a) + f(b))
o Adds f(xi)f(x_i) for each interior point
o Returns the total area

 The trapezoidal rule provides a simple way to numerically integrate functions.

 This section lays the foundation for converting the algorithm to a parallel CUDA
version.

5.10.2 A CUDA implementation

CUDA Implementation of the Trapezoidal Rule
Objective
 Parallelize the trapezoidal rule for numerical integration using CUDA.

23
Dept. of ISE
Parallel Computing-BCS702 Module 5

 The loop in the serial version is split into tasks:

1. Compute f(xi)f(x_i)f(xi)
2. Add result to the running total trap

Approach
 Each thread computes one iteration of the serial loop:
c
CopyEdit
int my_i = blockDim.x * blockIdx.x + threadIdx.x;
float my_x = a + my_i * h;
float my_trap = f(my_x);
trap += my_trap;
 Threads independently compute f(x_i) and contribute to the shared variable trap.

Problems Identified
1. Uninitialized Variables: h and trap are not initialized in the kernel.
2. Incorrect Index Range: my_i may go out of bounds; the loop should range from 1 to
n−1n-1n−1, but CUDA indices can start from 0.
3. Race Condition: trap is shared across threads. Multiple threads writing to it
simultaneously causes data races, leading to incorrect results.
4. Return Type Mismatch: CUDA kernels return void, so trap can’t be returned directly
like in the serial function.
5. Final Summation Needed: The values computed by threads must be summed
correctly on the host to get the total result.

Key Takeaway
While it's possible to parallelize the trapezoidal rule using CUDA, care must be taken to
handle indexing, memory access, and result aggregation to avoid errors and ensure
correctness.

5.10.3 Initialization, return value, and final update

Goal:
Handle initialization, returning results, and final updates in a CUDA version of the
trapezoidal rule.

24
Dept. of ISE
Parallel Computing-BCS702 Module 5

Proposed Approach:
 Use thread 0 (in block 0) to:
1. Initialize h and trap:
2. if (my_i == 0) {
3. h = (b - a)/n;
4. trap = 0.5 * (f(a) + f(b));
5. }
6. Perform final multiplication:
7. if (my_i == 0)
8. trap = trap * h;

Problems with This Approach:

1. Formal arguments (like h and trap) are private to each thread:
o Each CUDA thread has its own stack.
o Changes made by one thread to h or trap are not visible to other threads.
2. Thread synchronization is missing:
o No guarantee other threads will wait for initialization to finish before using h.
o Leads to incorrect or inconsistent results.

Solution Directions:
 For h:
o Compute on the host before kernel launch and pass as a kernel argument.
 For trap:
o Use shared memory or dynamically allocate device memory.
o Each thread writes its result to an element of a global array (e.g., trap[]),
which can be reduced on the host after kernel execution.

Key Takeaway:
To correctly parallelize the trapezoidal rule in CUDA:
 Avoid thread-local variables for shared computation.
 Use global memory for shared data like trap.
 Ensure synchronization or properly designed memory access to avoid race
conditions.

25
Dept. of ISE
Parallel Computing-BCS702 Module 5

A wrapper function is a function whose main purpose is to call another function. It can
perform any preparation needed for the call. It can also perform any additional work needed
after the call.

5.10.4 Using the Correct Threads

 The serial version’s loop runs from i = 1 to n - 1.
 So, in the CUDA kernel:
o We must ensure 0 < my_i < n before computing f(x) and adding to the sum.
o This avoids running the computation for thread 0 or out-of-bounds indices.

26
Dept. of ISE
Parallel Computing-BCS702 Module 5

o Correct thread check:

c
CopyEdit
if (0 < my_i && my_i < n) {
// Compute and add
}

5.10.5 Updating the Return Value and atomicAdd

 The shared memory location *trap_p is updated by multiple threads.
 A simple *trap_p += my_trap causes race conditions and unpredictable results.
 Solution: use CUDA’s atomicAdd() function, which performs an atomic (indivisible)
addition.
Why Atomic?
 Prevents multiple threads from interfering during read-modify-write operations.
 Ensures that every update to trap_p happens in isolation and is completed fully.
Function Signature:
c
CopyEdit
__device__ float atomicAdd(float* float_p, float val);
 Adds val to *float_p safely and atomically.
 Returns the previous value stored at float_p.

Key Takeaway:
 Use conditional execution to avoid invalid threads.
 Use atomicAdd() for safe parallel accumulation into shared variables like trap_p.

5.10.6 Performance of the CUDA trapezoidal rule

 The runtime is measured by timing the CUDA trapezoidal function call which
includes:
o Initialization of variables,
o Execution of the CUDA kernel,
o Final update of the result.

27
Dept. of ISE
Parallel Computing-BCS702 Module 5

 Timing is done using macros like GET_TIME to record start and finish times, and
then computing elapsed time.
 The same approach can be used to time the serial trapezoidal rule for comparison.
 Timings can vary due to system load, so it’s best to run the program multiple times
and report the mean or median.

Table 6.5 – Mean runtimes (seconds):

ARM Cortex- Nvidia Intel Core Nvidia GeForce GTX Titan

System
A15 GK20A i7 X

Clock 2.3 GHz 852 MHz 3.5 GHz 1.08 GHz

SMs, SM clock 4, 1.92 GHz - - 24, 3072

Run-time
33.6 20.7 4.48 3.08
(serial)

 The GPU runs are much faster than serial CPU runs.
 CUDA on the Nvidia GK20A GPU and GTX Titan X offers significant speedups over
ARM and Intel CPUs.
 Performance gains come from the massive parallelism on GPUs (thousands of
threads).

Key points:
 The CUDA trapezoidal rule can drastically reduce runtime compared to serial code.
 Multiple runs help average out variability in timing measurements.
 Hardware and GPU architecture greatly affect performance.

5.11 CUDA trapezoidal rule II: improving performance

5.11.1 Tree-structured communication:
CUDA trapezoidal rule II: improving performance
 To speed up the CUDA trapezoidal rule, use a tree-structured global sum instead
of a simple linear sum.

28
Dept. of ISE
Parallel Computing-BCS702 Module 5

 This method reduces serialization caused by multiple threads updating the same
variable trap_p with atomicAdd.
 A tree-structured sum aggregates partial sums in a hierarchical way, improving
performance by decreasing the number of sequential steps.

Tree-structured communication
 The global sum is done in multiple stages, where pairs of threads add their values,
then pairs of pairs, and so on, halving the number of active threads each step.
 Table 6.6 shows an example with 8 threads adding values in a tree-like manner,
gradually combining partial sums until one final result is produced.
 This approach reduces the time complexity from linear (n additions) to logarithmic
(log₂ n additions).
 The tree-structured approach reduces thread serialization and improves efficiency,
especially with many threads.

Key points:
 Tree-structured summation is more efficient for parallel reduction.
 It avoids bottlenecks caused by many threads updating the same variable.
 This method is common in parallel programming for summing or reducing values.

29
Dept. of ISE
Parallel Computing-BCS702 Module 5

5.11.2 Local variables, registers, shared and global

memory
CUDA uses a three-level memory hierarchy, each with different size and speed
characteristics:
1. Registers (fastest, smallest):
o Private to each thread.
o ~1 clock cycle to access.
o Best for storing local variables if enough space is available.
2. Shared Memory (medium speed and size):
o Shared among threads in the same block.
o Faster than global memory, slower than registers.
3. Global Memory (largest, slowest):
o Accessible by all threads on the GPU.

30
Dept. of ISE
Parallel Computing-BCS702 Module 5

o Much slower access (100–1000x slower than registers).

Local Variables:
 Local variables are ideally stored in registers.
 If there aren't enough registers, they are "spilled" to global memory (but in a
thread-private space).
 This hurts performance, since access to global memory is much slower.

Performance Tip:
 Use registers as much as possible to improve kernel performance.
 Minimize usage of shared and especially global memory unless necessary.
 Be aware of the limited register size, which can cause spilling if overused.

31
Dept. of ISE
Parallel Computing-BCS702 Module 5

5.11.3 Warps and warp shuffles

Warp shuffle functions allow threads within the same warp to exchange values through
registers, which is faster than using shared or global memory. This improves performance
significantly.

Key Concepts:
 Warp: A group of 32 threads (usually) in CUDA that execute together in SIMD
(Single Instruction, Multiple Data) fashion.
 _shfl_down_sync: A warp shuffle function that allows a thread to access the value
of a variable held by another thread in the same warp at a lower lane index.

Common Usage:
Used to implement tree-structured reductions (e.g., summing values across threads).
cpp
CopyEdit
val = __shfl_down_sync(mask, val, diff);
 mask: Bitmask indicating participating threads (usually set to 0xFFFFFFFF).
 val: Value to shuffle.
 diff: The offset by which to shuffle (e.g., 1, 2, 4...).
If a thread's lane + diff exceeds warpSize, the result is undefined.

Edge Cases & Notes:

 Be cautious of undefined values if lane + diff ≥ warpSize.
 The function only works within a single warp.
 CUDA 9.0 introduced the _sync versions (e.g., __shfl_down_sync), replacing earlier
non-sync versions.
 Use threadIdx.x % warpSize to determine a thread’s lane within the warp.

Use register-level communication via warp shuffle when reducing or aggregating data
within a warp — it’s faster and avoids race conditions common in shared/global memory.

32
Dept. of ISE
Parallel Computing-BCS702 Module 5

5.11.4 Implementing tree-structured global sum with a warp

shuffle
🧠 Goal:

To implement an efficient tree-structured global sum within a warp using the fast warp
shuffle instruction.

🔁 Code Overview:

__device__ float Warp_sum(float var) {

unsigned mask = 0xFFFFFFFF;
for (int diff = warpSize/2; diff > 0; diff = diff/2)
var += __shfl_down_sync(mask, var, diff);
return var;
}
 Performs a reduction within a warp (typically 32 threads).
 Uses __shfl_down_sync to pass and sum values across lanes in a binary tree
pattern.
 Works in log2(warpSize) steps.

🧮 How It Works:

 At each step, threads with lower lane IDs add values from higher lanes offset by
diff.
 diff starts at warpSize / 2 and halves every loop.
 Thread 0 ends up with the full sum of values from all threads in the warp.

📊 Figure 6.5 Explanation (assuming warpSize = 8):

 Shows how values propagate and are summed across lanes using a tree structure.
 E.g., in the first iteration (diff = 4), thread 0 adds value from thread 4, thread 1 from
thread 5, etc.
 The result is progressively accumulated in thread 0.

⚠️ Important Notes:

 This function only returns the full sum on lane 0 (i.e., thread with lane ID 0).

33
Dept. of ISE
Parallel Computing-BCS702 Module 5

 If all threads need the result, use __shfl_xor_sync instead (discussed in Exercise
6.6).

✅ Summary:

 Efficient, low-latency warp-wide reduction.

 Suitable for small group reductions where only one thread (usually lane 0) needs the
final result.
 Avoids shared/global memory, making it faster than conventional methods.

5.11.5 Shared memory and an alternative to the warp

shuffle
🧠 Purpose:

To implement a warp-level sum without using warp shuffle functions, which aren’t
available on GPUs with compute capability < 3.0.

💡 Key Idea:

Use shared memory instead of warp shuffle for inter-thread communication within a warp.
Threads in the same warp execute in lockstep (SIMD fashion), so updates to shared
memory can be synchronized implicitly.

📜 Function: Shared_mem_sum

__device__ float Shared_mem_sum(float shared_vals[]) {

int my_lane = threadIdx.x % warpSize;

34
Dept. of ISE
Parallel Computing-BCS702 Module 5

for (int diff = warpSize/2; diff > 0; diff = diff/2) {

int source = (my_lane + diff) % warpSize;
shared_vals[my_lane] += shared_vals[source];
}

return shared_vals[my_lane];
}
 Works like a tree-based reduction, but uses shared memory.
 Each thread reads from shared_vals and updates its own lane’s value.
 No race condition: all threads are executing in parallel and synchronously.

📊 Figure 6.6 Explanation:

 Demonstrates how values propagate across threads in a warp using shared

memory.
 Each thread adds values from higher-indexed threads in log2 steps, similar to tree-
structured sum.

⚠️ Notes:

 Although not technically a ―tree-structured‖ sum, this pattern is called a

dissemination sum.
 All threads perform the same updates, but only lane 0 ends up with the final total.
 Other threads may have intermediate or incorrect values unless corrected.

✅ Summary:

 A viable alternative to warp shuffle for older GPUs.

 Slightly slower than register-based shuffle, but still efficient.
 Relies on fast shared memory and warp-level synchronization.

35
Dept. of ISE
Parallel Computing-BCS702 Module 5

 Because all threads in a warp execute in lockstep (SIMD), there’s no performance

penalty for having every thread execute the same instruction.
 Hence, shared memory access within a warp is efficient.
 shared_vals can live in either global or shared memory, but using shared memory is
faster.
 You can allocate it statically or dynamically depending on your needs.
 Efficient usage of shared memory improves performance of parallel reductions like
global sums.


5.12 Implementation of trapezoidal rule with warpSize

thread blocks
This section puts together earlier concepts like warps, warp shuffles, and shared memory to
create optimized CUDA implementations of the trapezoidal rule.

5.12.1 Host Code

 The host code is nearly identical to the earlier CUDA versions.
 The only difference: there's no th_per_blk variable anymore.
 Assumes each thread block has exactly warpSize threads.

5.12.2 Kernel with Warp Shuffle

 The kernel uses warp shuffle to sum values across threads within a warp.
 Each thread computes a value (f(my_x)) and accumulates it into my_trap.

36
Dept. of ISE
Parallel Computing-BCS702 Module 5

 Warp_sum(my_trap) computes the sum of all my_trap values in the warp.

 Thread with lane ID 0 (i.e. threadIdx.x == 0) uses atomicAdd to add the warp’s sum
into the global total (trap_p).
Benefits:
 Reduces contention and synchronization overhead.
 Uses CUDA’s fast intra-warp communication via __shfl_down_sync.
 Only the first thread in each warp performs the atomic update to avoid race
conditions.

This optimized version uses warp-level parallelism and minimizes global memory updates.
It’s both faster and more scalable than the earlier naive implementations.

5.12.3Kernel with shared memory

It is almost identical to the version that uses the warp shuffle. The main differences are that it
declares an array of shared memory in Line7; it initializes this array in Lines11and14; and, of
course, the call to Shared mem_sum is passed this array rather than a scalar register.

37
Dept. of ISE
Parallel Computing-BCS702 Module 5

5.12.4 Performance

This section evaluates the performance of optimized trapezoidal rule implementations

(warp shuffle and shared memory) compared to the original CUDA and CPU versions.
The experiment uses:
 f(x)=x2+1f(x) = x^2 + 1
 Interval: [−3,3][-3, 3]
 220=1,048,5762^{20} = 1,048,576 trapezoids
 32 threads per block → 32,768 thread blocks
Table 6.8: Mean Run-Times (ms)

System Original Warp Shuffle Shared Memory

ARM Cortex-A15 33.6 – –

Nvidia GK20A 20.7 14.4 15.0

Intel Core i7 4.48 – –

Nvidia GTX Titan X 3.08 0.210 0.206

Key Takeaways:

38
Dept. of ISE
Parallel Computing-BCS702 Module 5

 Warp Shuffle and Shared Memory implementations outperform the original

CUDA version significantly:
o Nvidia GK20A:
 Warp Shuffle = ~70% of original time
 Shared Memory = ~72% of original time
o GTX Titan X:
 Both new versions run in < 7% of the original time.
 Warp shuffle is slightly slower than shared memory here.
 Conclusion:
o Both optimizations drastically reduce run-time.
o Performance gains are system-dependent, but shared memory can
match or exceed warp shuffle performance even though it is generally
slower than registers.

5.13 CUDA trapezoidal rule III: blocks with more than one
warp

This section explores extending CUDA implementations of the trapezoidal rule by using
larger thread blocks—beyond a single warp (i.e., more than 32 threads, up to 1024). This
approach boosts flexibility and performance by better utilizing GPU resources.

5.13.1: Using __syncthreads() to Coordinate Warps

Problem:
When summing results across multiple warps in a block, there's a race condition:
 Warp 0 may begin summing partial results before other warps (e.g., warp 1) finish
their calculations.
 This leads to incorrect final sums if warps do not synchronize.
Solution:
Use CUDA’s barrier synchronization function __syncthreads():
c
CopyEdit
__syncthreads(); // Ensures all threads reach this point before proceeding
 Forces all threads in a block to wait until every other thread reaches the same point
in code.

39
Dept. of ISE
Parallel Computing-BCS702 Module 5

 Ensures all warp-level sums are computed before a designated warp (e.g., warp 0)
begins the final summation.

Key Takeaways:
 Larger blocks (e.g., 1024 threads = 32 warps) allow parallel summation across
warps.
 __syncthreads() is critical to coordinate these warps and prevent race conditions.
 Synchronization ensures correctness when aggregating results across multiple warps
in a block.
 Developers must be careful using __syncthreads() with conditional execution, as
divergent threads may never reach the sync point, causing deadlock.
This approach enables more scalable, performant implementations by combining intra-warp
optimizations (like warp shuffle) with inter-warp coordination using __syncthreads().
The second important caveat about __syncthreads() is that it only synchronizes threads
within the same block, not across different blocks. Even if every thread in the entire grid
calls __syncthreads(), threads in different blocks will still operate independently.
Therefore, __syncthreads() cannot be used to synchronize threads across a whole
grid.

5.13.2 More shared memory

 Instead of relying on warp shuffles, shared memory can be used to allow threads in
warp 0 to access registers from other warps. A shared array sized for the number of
warps stores warp sums, enabling efficient communication.

5.13.3 Shared memory warp sums:

 When using shared memory for warp sums, each warp stores its threads’ partial
sums in a shared subarray. This subarray holds contributions from all threads in the
warp. This method requires declaring a sufficiently large shared array to hold all
thread contributions across warps.
Overall, the use of shared memory enables flexible and efficient aggregation of sums across
warps, overcoming limitations of warp shuffle operations.

5.13.4 Shared memory banks

This section explains how Nvidia's shared memory is organized into "banks" to allow
simultaneous access by threads in a warp:
 Shared memory is divided into 32 banks (16 for older GPUs with compute capability
< 2.0).
 Each bank stores consecutive elements of arrays like thread_calcs.
 When threads access different banks simultaneously, access is parallel and fast.
 If multiple threads access the same bank but different locations, accesses are
serialized, causing delays.

40
Dept. of ISE
Parallel Computing-BCS702 Module 5

 If multiple threads access the same location in a bank, the value is broadcast
simultaneously, which is efficient.
 Proper memory layout is crucial to avoid serialization and achieve maximum speed.
 The CUDA programming guide provides full details on bank conflicts and
optimization.
In short, understanding and optimizing shared memory bank usage is key for high-
performance CUDA programs.

41
Dept. of ISE
Parallel Computing-BCS702 Module 5

42
Dept. of ISE