Add GPU programming patterns tutorials (#1918)
Update projects/hip/docs/tutorial/programming-patterns/atomic_operations_histogram.rst WIP Co-authored-by: Julia Jiang <56359287+jujiang-del@users.noreply.github.com>
Этот коммит содержится в:
Istvan Kiss
GitHub
+310
@@ -0,0 +1,310 @@
|
||||
.. meta::
|
||||
:description: HIP atomic operations histogram tutorial
|
||||
:keywords: AMD, ROCm, HIP, atomic operations, GPU programming, histogram, synchronization primitives
|
||||
|
||||
*******************************************************************************
|
||||
Atomic operations: Histogram tutorial
|
||||
*******************************************************************************
|
||||
|
||||
In GPU programming, a core design principle is to **avoid simultaneous writes to
|
||||
the same memory address by multiple threads**. When multiple threads write to
|
||||
the same location without proper synchronization, this creates a
|
||||
**race condition**, where the final result depends on unpredictable thread
|
||||
execution order.
|
||||
|
||||
Unlike CPUs, GPUs are designed for high-throughput parallel execution with
|
||||
relaxed memory consistency models and limited cache coherence mechanisms. This
|
||||
architectural choice maximizes bandwidth and scalability but introduces
|
||||
challenges when multiple threads need to safely update shared state.
|
||||
|
||||
This tutorial demonstrates how to safely handle **concurrent memory updates**
|
||||
using **atomic operations**, illustrated through the practical example of
|
||||
computing an image brightness histogram on the GPU.
|
||||
|
||||
.. include:: ../prerequisites.rst
|
||||
|
||||
Race condition
|
||||
==============
|
||||
|
||||
A **race condition** occurs when two or more threads attempt to
|
||||
read-modify-write the same memory location concurrently without proper
|
||||
synchronization. Because GPU threads execute asynchronously across multiple
|
||||
cores (compute units), concurrent writes can interleave unpredictably,
|
||||
leading to incorrect results.
|
||||
|
||||
For example, if two threads simultaneously attempt:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
histogram[bin] = histogram[bin] + 1;
|
||||
|
||||
both may read the same old value before either writes back,
|
||||
resulting in only one increment being reflected. This results in **lost updates**
|
||||
and **nondeterministic output**, which must be avoided.
|
||||
|
||||
Histogram
|
||||
=========
|
||||
|
||||
A **histogram** partitions continuous data into discrete intervals called
|
||||
**bins** and counts how many data points fall into each bin. In image processing,
|
||||
a histogram typically represents the **distribution of pixel intensities** for
|
||||
example brightness or color channel values.
|
||||
|
||||
The histogram algorithm can be expressed as:
|
||||
|
||||
.. math::
|
||||
|
||||
H[b] = \sum_{i=1}^{N} \delta(b - \lfloor f(x_i) \rfloor)
|
||||
|
||||
where :math:`f(x_i)` maps each data value to its corresponding bin index
|
||||
:math:`b`, and :math:`\delta()` is 1 when the value belongs to bin :math:`b` and
|
||||
0 otherwise.
|
||||
|
||||
The basic computational steps are:
|
||||
|
||||
1. Iterate through all pixels (or data points).
|
||||
2. Determine the appropriate bin for each value.
|
||||
3. Increment that bin’s count.
|
||||
|
||||
In a serial CPU program, this is straightforward. On a GPU, thousands of threads
|
||||
may attempt to increment the same bin concurrently, leading to **race
|
||||
conditions** unless atomic synchronization is used.
|
||||
|
||||
The Challenge in parallel context
|
||||
---------------------------------
|
||||
|
||||
When multiple threads attempt to increment the same bin:
|
||||
|
||||
* One thread’s update can overwrite another’s pending increment.
|
||||
|
||||
* Memory coherence cannot guarantee ordered visibility across thread blocks.
|
||||
|
||||
* The final result may be inconsistent or incorrect.
|
||||
|
||||
This necessitates synchronization mechanisms to ensure that updates occur in a
|
||||
**mutually exclusive** manner without introducing high overhead.
|
||||
|
||||
Atomic operations
|
||||
=================
|
||||
|
||||
An **atomic operation** ensures that a compound operation — typically a
|
||||
read-modify-write sequence — executes as an **indivisible unit**. From the
|
||||
programmer’s perspective, atomicity guarantees that no other thread can observe
|
||||
a partially completed operation.
|
||||
|
||||
Formally, an operation :math:`O(x)` on shared variable :math:`x` is **atomic**
|
||||
if its execution satisfies:
|
||||
|
||||
.. math::
|
||||
|
||||
\forall T_i, T_j, \text{ the effects of } O(x) \text{ appear serializable.}
|
||||
|
||||
That is, all threads observe results as if operations occurred in a single,
|
||||
sequential order.
|
||||
|
||||
Mechanics
|
||||
---------
|
||||
|
||||
Atomic operations on GPUs are implemented in hardware through a **memory
|
||||
arbitration unit** that locks a cache line, performs the modification, and
|
||||
releases the lock. This ensures correctness even under massive parallelism.
|
||||
|
||||
When a thread performs an atomic operation:
|
||||
|
||||
1. The target memory location is temporarily locked.
|
||||
2. The value is fetched and updated.
|
||||
3. The update is written back, and the lock is released.
|
||||
|
||||
No other thread can modify the same memory location during this sequence.
|
||||
|
||||
Atomic functions
|
||||
----------------
|
||||
|
||||
HIP provides a wide set of atomic primitives to synchronize updates to shared
|
||||
memory or global memory locations:
|
||||
|
||||
.. list-table::
|
||||
:header-rows: 1
|
||||
:widths: 20 80
|
||||
|
||||
* - Operation
|
||||
- Description
|
||||
* - ``atomicAdd``
|
||||
- Atomically adds a value to a memory location and returns the old value.
|
||||
* - ``atomicSub``
|
||||
- Atomically subtracts a value.
|
||||
* - ``atomicExch``
|
||||
- Atomically exchanges values between a register and memory.
|
||||
* - ``atomicCAS``
|
||||
- Performs an atomic compare-and-swap; fundamental for implementing locks.
|
||||
* - ``atomicMax`` / ``atomicMin``
|
||||
- Updates to the maximum or minimum of two values.
|
||||
* - ``atomicInc`` / ``atomicDec``
|
||||
- Atomically increments or decrements a counter, wrapping at a boundary.
|
||||
|
||||
Atomic operations in kernels can operate on block scope (shared memory),
|
||||
device scope (global memory), or system scope (system memory), depending on
|
||||
:doc:`hardware support <rocm:reference/gpu-atomics-operation>`.
|
||||
|
||||
For more information, please check :ref:`atomic functions <atomic functions>`.
|
||||
|
||||
Image brightness histogram
|
||||
==========================
|
||||
|
||||
We will compute a histogram that captures the **distribution of pixel
|
||||
brightness** in an RGB image. The algorithm:
|
||||
|
||||
1. Reads image data in **channel-height-width** format.
|
||||
2. Converts RGB values to grayscale brightness.
|
||||
3. Maps brightness to a histogram bin.
|
||||
4. Atomically increments the corresponding bin counter.
|
||||
|
||||
Kernel implementation
|
||||
---------------------
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
__global__ void calculateHistogram(float* imageData, int* histogram,
|
||||
int width, int height,
|
||||
int channels, int numBins)
|
||||
{
|
||||
int x = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
int y = blockIdx.y * blockDim.y + threadIdx.y;
|
||||
|
||||
if (x >= width || y >= height)
|
||||
return;
|
||||
|
||||
int idx = (y * width + x) * channels;
|
||||
float brightness = 0.0f;
|
||||
|
||||
for (int c = 0; c < channels; ++c)
|
||||
brightness += imageData[idx + c];
|
||||
|
||||
brightness /= channels; // Normalize to [0, 1]
|
||||
int bin = static_cast<int>(brightness * numBins);
|
||||
|
||||
// Atomic increment to avoid race conditions
|
||||
atomicAdd(&histogram[bin], 1);
|
||||
}
|
||||
|
||||
Thread identification
|
||||
~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
Each thread computes one pixel’s contribution using its 2D thread and block
|
||||
indices:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
int x = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
int y = blockIdx.y * blockDim.y + threadIdx.y;
|
||||
|
||||
This mapping provides a 1:1 correspondence between threads and pixels, making
|
||||
the computation naturally parallel.
|
||||
|
||||
Brightness computation
|
||||
~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
Each pixel’s brightness is computed as the arithmetic mean of its RGB channels:
|
||||
|
||||
.. math::
|
||||
:nowrap:
|
||||
|
||||
\[
|
||||
I'(x, y) = \frac{R + G + B}{3}
|
||||
\]
|
||||
|
||||
This value is then normalized to [0, 1] and mapped to one of `numBins`
|
||||
histogram intervals.
|
||||
|
||||
Safe histogram update
|
||||
~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
The key step is:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
atomicAdd(&histogram[bin], 1);
|
||||
|
||||
This ensures that even if thousands of threads map to the same bin, each
|
||||
increment is serialized correctly, maintaining an accurate bin count.
|
||||
|
||||
Performance characteristics
|
||||
===========================
|
||||
|
||||
Benefits
|
||||
--------
|
||||
|
||||
* **Correctness under parallel updates:** Ensures race-free accumulation.
|
||||
|
||||
* **Simplified synchronization:** No explicit locks or barriers needed.
|
||||
|
||||
* **Hardware-level efficiency:** Implemented directly in the GPU memory
|
||||
subsystem.
|
||||
|
||||
Limitations
|
||||
-----------
|
||||
|
||||
While atomic operations guarantee correctness, they can **serialize execution**
|
||||
when multiple threads target the same memory address. This causes contention and
|
||||
reduces effective parallelism.
|
||||
|
||||
Typical performance degradation sources include:
|
||||
|
||||
* **Hot bins:** When many pixels fall into a small subset of bins.
|
||||
|
||||
* **Global memory atomics:** Global memory atomics are slower than shared memory
|
||||
atomics due to higher access latency.
|
||||
|
||||
* **Warp serialization:** Threads within a warp waiting for the same atomic
|
||||
target serialize.
|
||||
|
||||
Best practices
|
||||
==============
|
||||
|
||||
1. **Apply atomic operations only where necessary**
|
||||
|
||||
Atomic instructions serialize access to a memory location and use can
|
||||
diminish SIMT parallel efficiency and increase warp stalls. Restrict atomic
|
||||
usage to code paths where data races cannot be eliminated through algorithmic
|
||||
restructuring.
|
||||
|
||||
2. **Minimize contention**
|
||||
|
||||
High contention on a single address or a small set of addresses leads to
|
||||
serialization. Distribute writes across independent memory locations.
|
||||
|
||||
3. **Leverage shared memory**
|
||||
|
||||
Use fast, low-latency shared memory to aggregate partial results within a
|
||||
block before issuing a single atomic update to global memory.
|
||||
|
||||
4. **Validate correctness**
|
||||
|
||||
Validate the numerical and logical correctness of GPU kernels by comparing
|
||||
against single-threaded or deterministic multi-threaded CPU baselines.
|
||||
|
||||
5. **Profile regularly**
|
||||
|
||||
GPU performance is highly sensitive to thread divergence, memory-access
|
||||
patterns, and workload distribution. Regularly use profiling tools such as
|
||||
:doc:`rocprofv3<rocprofiler-sdk:how-to/using-rocprofv3>` or
|
||||
:doc:`ROCm compute profiler<rocprofiler-compute:how-to/profile/mode>` to
|
||||
examine warp-level execution efficiency, memory-coalescing behavior,
|
||||
occupancy, and atomic throughput bottlenecks.
|
||||
|
||||
Conclusion
|
||||
==========
|
||||
|
||||
Atomic operations provide a low-level synchronization mechanism that allows
|
||||
correct and deterministic parallel updates to shared data structures. In the
|
||||
histogram example, :cpp:func:`atomicAdd` ensures that all threads safely
|
||||
contribute to their corresponding bins, preventing race conditions.
|
||||
|
||||
While atomics incur some serialization overhead, they are indispensable for
|
||||
algorithms that require concurrent accumulation or counting. By applying
|
||||
techniques like privatization and reduction, developers can achieve both
|
||||
**correctness** and **high performance** on modern GPUs.
|
||||
|
||||
Atomic operations form the foundation for more advanced synchronization
|
||||
patterns, including parallel reductions, prefix sums, and graph traversal, and
|
||||
are essential for developing scalable, data-parallel GPU algorithms.
|
||||
@@ -0,0 +1,646 @@
|
||||
.. meta::
|
||||
:description: CPU-GPU cooperative computing at K-means clustering
|
||||
:keywords: AMD, ROCm, HIP, CPU-GPU cooperative computing, K-means, clustering, K-means clustering
|
||||
|
||||
*******************************************************************************
|
||||
CPU-GPU cooperative computing: K-means clustering tutorial
|
||||
*******************************************************************************
|
||||
|
||||
Modern heterogeneous systems combine CPUs and GPUs to maximize computational
|
||||
throughput. GPUs provide massive data-parallel performance, whereas CPUs excel
|
||||
at complex control logic, and latency-sensitive serial tasks. Many real-world
|
||||
algorithms—including unsupervised clustering—contain both parallelizable and
|
||||
inherently sequential components.
|
||||
|
||||
This tutorial demonstrates a hybrid CPU–GPU cooperative execution model using
|
||||
the K-means clustering algorithm, showcasing partitioned workload distribution,
|
||||
memory management, and performance optimization strategies in heterogeneous
|
||||
architectures.
|
||||
|
||||
.. include:: ../prerequisites.rst
|
||||
|
||||
CPU–GPU cooperative computing
|
||||
=============================
|
||||
|
||||
Modern computing platforms integrate **central processing units (CPUs)** and
|
||||
**graphics processing units (GPUs)** into unified heterogeneous systems.
|
||||
Each processor type offers distinct architectural strengths:
|
||||
|
||||
* **CPUs** feature a small number of complex cores optimized for sequential
|
||||
execution, branching, and latency-sensitive control flow.
|
||||
|
||||
* **GPUs** consist of thousands of simpler cores optimized for throughput and
|
||||
massive data-level parallelism.
|
||||
|
||||
**Cooperative computing** refers to a programming paradigm that combines these
|
||||
capabilities in a coordinated execution model, where both CPU and GPU perform
|
||||
complementary portions of a workload. Rather than treating the GPU as a passive
|
||||
accelerator, this approach distributes computation dynamically between devices
|
||||
to maximize total system utilization.
|
||||
|
||||
K-means clustering
|
||||
==================
|
||||
|
||||
K-means is an unsupervised machine learning algorithm that partitions a dataset
|
||||
into :math:`k` clusters (groups of similar data points) by minimizing
|
||||
intra-cluster variance. It iteratively refines cluster assignments and centroid
|
||||
(center point of a cluster) locations until convergence.
|
||||
|
||||
The optimization objective is defined as:
|
||||
|
||||
.. math::
|
||||
|
||||
\min_{C_1, \dots, C_k} \sum_{i=1}^{k} \sum_{x_j \in C_i} \|x_j - \mu_i\|^2
|
||||
|
||||
where :math:`\mu_i` is the centroid of cluster :math:`C_i`.
|
||||
|
||||
K-means is frequently used in:
|
||||
|
||||
* **Customer segmentation**: Grouping customers by behavior patterns
|
||||
|
||||
* **Image compression**: Reducing color palette by clustering similar colors
|
||||
|
||||
* **Anomaly detection**: Identifying outliers that don't fit clusters
|
||||
|
||||
* **Document clustering**: Organizing similar documents together
|
||||
|
||||
* **Feature engineering**: Creating new features based on cluster membership
|
||||
|
||||
Algorithm
|
||||
=========
|
||||
|
||||
The K-means algorithm iteratively refines a partition of a dataset into
|
||||
:math:`k` clusters by alternating between two primary computational phases:
|
||||
**assignment** and **update**. Each iteration minimizes the total within-cluster
|
||||
variance, driving the system toward convergence where centroid movement or
|
||||
membership changes fall below a defined threshold.
|
||||
|
||||
Iterative procedure
|
||||
-------------------
|
||||
|
||||
1. **Initialization**: Select initial centroids (random or via K-means++).
|
||||
|
||||
2. **Assignment**: Assign each data point to the nearest centroid.
|
||||
|
||||
3. **Update**: Recalculate each centroid as the mean of its assigned members.
|
||||
|
||||
4. **Convergence**: Repeat until centroids stabilize or a maximum iteration
|
||||
limit is reached.
|
||||
|
||||
|
||||
These phases alternate until the solution converges or a maximum iteration count
|
||||
is reached.
|
||||
|
||||
Initialization
|
||||
~~~~~~~~~~~~~~
|
||||
|
||||
The algorithm begins by selecting initial centroid positions. This step
|
||||
significantly influences both convergence rate and clustering quality.
|
||||
|
||||
* **Random initialization**: centroids are randomly chosen from the dataset.
|
||||
* **K-means++**: probabilistic seeding to spread centroids across data space.
|
||||
* **Domain-specific heuristics**: custom initializations for structured data.
|
||||
|
||||
Initial centroid diversity reduces the likelihood of poor local minima and
|
||||
improves overall stability.
|
||||
|
||||
Assignment
|
||||
~~~~~~~~~~
|
||||
|
||||
For each data point :math:`x_i`, the algorithm computes the Euclidean distance
|
||||
to each centroid :math:`\mu_j` and assigns the point to the nearest cluster:
|
||||
|
||||
.. math::
|
||||
|
||||
C_i = \arg \min_j \|x_i - \mu_j\|^2
|
||||
|
||||
* Each point–centroid distance calculation is independent.
|
||||
|
||||
* Ideal for SIMD/SIMT architectures (GPU execution).
|
||||
|
||||
* Dominated by dense floating-point arithmetic and memory bandwidth utilization.
|
||||
|
||||
This phase represents the **embarrassingly parallel** portion of the algorithm
|
||||
and provides the largest opportunity for GPU acceleration.
|
||||
|
||||
Update
|
||||
~~~~~~
|
||||
|
||||
Once all data points are assigned, new centroid positions are computed as the
|
||||
mean of their corresponding cluster members:
|
||||
|
||||
.. math::
|
||||
|
||||
\mu_j = \frac{1}{|C_j|} \sum_{x_i \in C_j} x_i
|
||||
|
||||
* Requires aggregation and division per cluster.
|
||||
|
||||
* Involves variable membership counts across clusters.
|
||||
|
||||
* Reduction-heavy and branch-divergent.
|
||||
|
||||
Although reduction can be implemented on GPUs, small :math:`k` values and
|
||||
irregular membership sizes typically make the CPU more efficient for this phase,
|
||||
given its superior cache hierarchy and flexible control flow.
|
||||
|
||||
Convergence
|
||||
~~~~~~~~~~~
|
||||
|
||||
After the update phase, convergence is evaluated using one of the following
|
||||
criteria:
|
||||
|
||||
* No change in point memberships.
|
||||
|
||||
* Centroid displacement below a user-defined threshold.
|
||||
|
||||
* Iteration count exceeds maximum limit.
|
||||
|
||||
If the convergence condition is not met, the assignment and update phases
|
||||
repeat.
|
||||
|
||||
Summary of cooperative execution
|
||||
--------------------------------
|
||||
|
||||
1. **GPU:** performs parallel distance evaluations and membership assignments.
|
||||
2. **CPU:** executes centroid averaging and convergence checks.
|
||||
3. **Synchronization:** only essential data (centroids and membership
|
||||
arrays) is exchanged per iteration.
|
||||
|
||||
This division of labor minimizes data movement and maximizes hardware
|
||||
utilization. The resulting hybrid implementation combines GPU throughput for
|
||||
massive data-parallel operations with CPU efficiency for aggregation and control
|
||||
tasks, enabling scalable clustering performance on heterogeneous systems.
|
||||
|
||||
Implementation
|
||||
==============
|
||||
|
||||
This implementation follows a hybrid CPU/GPU design to accelerate the most
|
||||
computationally expensive phase of K-means: assigning each data point to its
|
||||
nearest centroid. The GPU performs distance calculations in parallel, while the
|
||||
CPU handles the centroid recomputation, where sequential reductions are
|
||||
efficient and data sizes are small. The algorithm iterates between GPU-based
|
||||
membership updates and CPU-based centroid averaging until convergence or a
|
||||
maximum iteration count is reached.
|
||||
|
||||
Data Structures
|
||||
---------------
|
||||
|
||||
The implementation stores all data in simple, contiguous arrays for efficient
|
||||
memory access on both CPU and GPU. Data points and centroids are represented as
|
||||
flattened ``std::vector<float>`` arrays, while cluster assignments are stored as
|
||||
``std::vector<int>``.
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
// length: Number of data points to cluster
|
||||
// dimension: Number of features per data point
|
||||
// k: Number of clusters
|
||||
|
||||
std::vector<float> data; // Data points (length * dimension)
|
||||
std::vector<float> centroids; // Centroid positions (k * dimension)
|
||||
std::vector<int> memberships; // Cluster assignments (length)
|
||||
|
||||
Main Loop
|
||||
---------
|
||||
|
||||
The core K-means iteration:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
// length is an integer for the number of entries to be clustered.
|
||||
// dimension is an integer for the number of properties of each entry.
|
||||
// k is an integer that determines the number of clusters.
|
||||
|
||||
std::vector<float> centroids = initializeCentroids(length * dimension, k);
|
||||
std::vector<int> memberships(length, 0);
|
||||
|
||||
for (int iteration = 0; iteration < maxIterations; ++iteration) {
|
||||
// Determine the cluster that each entry belongs to.
|
||||
// The function returns how many entries changed membership.
|
||||
int membershipChanges = updateMembership(data, centroids, memberships);
|
||||
|
||||
// Converge checking.
|
||||
if (membershipChanges == 0) {
|
||||
break;
|
||||
}
|
||||
|
||||
// Calculate new centroids.
|
||||
std::vector<Point> newCentroids = centroids;
|
||||
updateCentroid(data, newCentroids, memberships);
|
||||
centroids = newCentroids;
|
||||
}
|
||||
|
||||
**Loop structure:**
|
||||
|
||||
1. Update memberships using GPU
|
||||
2. Check for convergence (no changes)
|
||||
3. Update centroids using CPU
|
||||
4. Repeat until convergence or max iterations
|
||||
|
||||
GPU membership update function
|
||||
------------------------------
|
||||
|
||||
The ``updateMembership`` function serves as the interface between CPU and GPU:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
int updateMembership(float* data, float* centroids, int* membership,
|
||||
int dataSize, int dimension, int k) {
|
||||
// gpuData is allocated and copied to the GPU earlier.
|
||||
|
||||
float *gpuCentroids, *gpuMembership;
|
||||
|
||||
// Allocate GPU memory
|
||||
hipMalloc(&gpuCentroids, k * dimension * sizeof(float));
|
||||
hipMalloc(&gpuMembership, dataSize * sizeof(int));
|
||||
|
||||
// Copy data from CPU to GPU
|
||||
hipMemcpy(gpuCentroids, centroids, k * dimension * sizeof(float),
|
||||
hipMemcpyHostToDevice);
|
||||
|
||||
// Calculate the sizes for the kernel launch
|
||||
int localSize = 256;
|
||||
int globalSize = (dataSize + localSize - 1) / localSize;
|
||||
|
||||
// Launch the kernel
|
||||
updateMembershipGPU<<<globalSize, localSize>>>(
|
||||
gpuData, gpuCentroids, gpuMembership, dataSize, dimension, k);
|
||||
hipDeviceSynchronize();
|
||||
|
||||
// Create CPU Data to hold the results
|
||||
std::vector<int> cpuNewMembership(dataSize);
|
||||
|
||||
// Copy GPU data back to CPU
|
||||
hipMemcpy(cpuNewMembership.data(), gpuMembership,
|
||||
dataSize * sizeof(int), hipMemcpyDeviceToHost);
|
||||
|
||||
// Count membership updates
|
||||
int membershipUpdate = 0;
|
||||
for (int i = 0; i < dataSize; ++i) {
|
||||
if (membership[i] != cpuNewMembership[i]) {
|
||||
membershipUpdate++;
|
||||
membership[i] = cpuNewMembership[i]; // Update the original
|
||||
}
|
||||
}
|
||||
|
||||
// Free GPU memory
|
||||
hipFree(gpuCentroids);
|
||||
hipFree(gpuMembership);
|
||||
|
||||
return membershipUpdate;
|
||||
}
|
||||
|
||||
Step 1: GPU memory allocation
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
float *gpuCentroids, *gpuMembership;
|
||||
hipMalloc(&gpuCentroids, k * dimension * sizeof(float));
|
||||
hipMalloc(&gpuMembership, dataSize * sizeof(int));
|
||||
|
||||
Allocate space on the GPU for:
|
||||
|
||||
* **Centroids**: :math:`k` centroids, each with dimension features
|
||||
* **Membership assignments**: One cluster ID per data point
|
||||
|
||||
Step 2: Data transfer to GPU
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
hipMemcpy(gpuCentroids, centroids, k * dimension * sizeof(float),
|
||||
hipMemcpyHostToDevice);
|
||||
|
||||
Transfer the current centroid positions from CPU to GPU. Note that the data
|
||||
points (``gpuData``) are already on the GPU from earlier initialization.
|
||||
|
||||
Step 3: Kernel configuration
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
int localSize = 256;
|
||||
int globalSize = (dataSize + localSize - 1) / localSize;
|
||||
|
||||
* ``localSize``: 256 threads per block (common choice)
|
||||
* ``globalSize``: Number of blocks needed to cover all data points
|
||||
* Rounding up ensures we process all data points
|
||||
|
||||
Step 4: Kernel launch
|
||||
~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
updateMembershipGPU<<<globalSize, localSize>>>(
|
||||
gpuData, gpuCentroids, gpuMembership, dataSize, dimension, k);
|
||||
hipDeviceSynchronize();
|
||||
|
||||
Launch the GPU kernel to compute cluster assignments in parallel, then wait for
|
||||
completion.
|
||||
|
||||
Step 5: Retrieve results
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
std::vector<int> cpuNewMembership(dataSize);
|
||||
hipMemcpy(cpuNewMembership.data(), gpuMembership,
|
||||
dataSize * sizeof(int), hipMemcpyDeviceToHost);
|
||||
|
||||
Copy the new membership assignments back from GPU to CPU.
|
||||
|
||||
Step 6: Count changes
|
||||
~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
int membershipUpdate = 0;
|
||||
for (int i = 0; i < dataSize; ++i) {
|
||||
if (membership[i] != cpuNewMembership[i]) {
|
||||
membershipUpdate++;
|
||||
membership[i] = cpuNewMembership[i];
|
||||
}
|
||||
}
|
||||
|
||||
Count how many data points changed clusters. This value determines if the algorithm has converged.
|
||||
|
||||
Step 7: Cleanup
|
||||
~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
hipFree(gpuCentroids);
|
||||
hipFree(gpuMembership);
|
||||
return membershipUpdate;
|
||||
|
||||
Free temporary GPU memory and return the change count.
|
||||
|
||||
Kernel implementation
|
||||
~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
__global__ void updateMembershipGPU(
|
||||
float* data, float* centroids, int* membership,
|
||||
int dataSize, int dimension, int k)
|
||||
{
|
||||
int tid = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
|
||||
if (tid < dataSize) {
|
||||
float minDistance = INFINITY;
|
||||
int bestCluster = 0;
|
||||
|
||||
// Find nearest centroid
|
||||
for (int cluster = 0; cluster < k; ++cluster) {
|
||||
float distance = 0.0f;
|
||||
|
||||
// Calculate Euclidean distance
|
||||
for (int d = 0; d < dimension; ++d) {
|
||||
float diff = data[tid * dimension + d] -
|
||||
centroids[cluster * dimension + d];
|
||||
distance += diff * diff;
|
||||
}
|
||||
|
||||
if (distance < minDistance) {
|
||||
minDistance = distance;
|
||||
bestCluster = cluster;
|
||||
}
|
||||
}
|
||||
|
||||
membership[tid] = bestCluster;
|
||||
}
|
||||
}
|
||||
|
||||
**Kernel operation:**
|
||||
|
||||
1. Each thread processes one data point
|
||||
2. Calculates distance to all :math:`k` centroids
|
||||
3. Assigns point to nearest centroid
|
||||
4. Stores result in membership array
|
||||
|
||||
CPU centroid update
|
||||
-------------------
|
||||
|
||||
The CPU handles the averaging operation:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
void updateCentroid(float* data, float* centroids, int* membership,
|
||||
int dataSize, int dimension, int k)
|
||||
{
|
||||
std::vector<int> counts(k, 0);
|
||||
std::vector<float> sums(k * dimension, 0.0f);
|
||||
|
||||
// Accumulate sums for each cluster
|
||||
for (int i = 0; i < dataSize; ++i) {
|
||||
int cluster = membership[i];
|
||||
counts[cluster]++;
|
||||
|
||||
for (int d = 0; d < dimension; ++d) {
|
||||
sums[cluster * dimension + d] += data[i * dimension + d];
|
||||
}
|
||||
}
|
||||
|
||||
// Calculate averages (new centroids)
|
||||
for (int cluster = 0; cluster < k; ++cluster) {
|
||||
if (counts[cluster] > 0) {
|
||||
for (int d = 0; d < dimension; ++d) {
|
||||
centroids[cluster * dimension + d] =
|
||||
sums[cluster * dimension + d] / counts[cluster];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
The centroid update requires:
|
||||
|
||||
1. **Reduction operation**: Sum all points in each cluster.
|
||||
2. **Variable-sized groups**: Clusters have different numbers of points.
|
||||
3. **Division**: Calculate average (not easily parallelized).
|
||||
|
||||
While GPUs can perform reductions, for this workload:
|
||||
|
||||
* The number of clusters (:math:`k`) is typically small (< 100).
|
||||
* The CPU can efficiently handle this sequential aggregation.
|
||||
* Avoiding GPU complexity keeps code simpler.
|
||||
|
||||
Data transfer considerations
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
**Data already on CPU:**
|
||||
|
||||
* Membership array was just copied back.
|
||||
* Data points can be kept on CPU for small datasets.
|
||||
* Centroids are small (k × dimension values).
|
||||
|
||||
**Transfer costs:**
|
||||
|
||||
* Only centroids need to copy back to GPU.
|
||||
* Much smaller than the full dataset.
|
||||
* Overhead justified by parallel speedup in assignment phase.
|
||||
|
||||
Best practices
|
||||
==============
|
||||
|
||||
This section outlines recommended practices for implementing an efficient
|
||||
GPU-accelerated Breadth-First Search (BFS). It highlights design principles,
|
||||
memory-management strategies, and debugging techniques that help ensure
|
||||
correctness, maintainability, and high performance when mapping BFS onto modern
|
||||
GPU architectures.
|
||||
|
||||
Design principles
|
||||
-----------------
|
||||
|
||||
1. **Identify parallelism**
|
||||
|
||||
Decompose the K-means workflow into independent, data-parallel kernels such
|
||||
as distance computation. Minimize sequential dependencies in assignment and
|
||||
reduction steps.
|
||||
|
||||
2. **Minimize host–device data movement**
|
||||
|
||||
Maintain dataset and intermediate buffers resident on the GPU across
|
||||
iterations. Transfer only convergence metrics or centroids when absolutely
|
||||
necessary.
|
||||
|
||||
3. **Batch and fuse operations**
|
||||
|
||||
Combine smaller kernels such as distance computation and assignment, or
|
||||
process multiple mini-batches per kernel launch to reduce kernel invocation
|
||||
and PCIe transfer overhead.
|
||||
|
||||
4. **Overlap communication with computation**
|
||||
|
||||
Utilize asynchronous HIP streams to overlap host–device memory transfers with
|
||||
GPU kernel execution. Employ pinned (page-locked) memory for faster DMA
|
||||
transfers.
|
||||
|
||||
Memory strategy
|
||||
---------------
|
||||
|
||||
1. **Persistent device allocations**
|
||||
|
||||
Allocate GPU memory once before iterative processing. Reuse buffers such as
|
||||
those for centroids, assignments, and temporary reductions to avoid the
|
||||
overhead of repeated :cpp:func:`hipMalloc` and :cpp:func:`hipFree` calls.
|
||||
|
||||
2. **Data locality optimization**
|
||||
|
||||
Co-locate data structures according to access frequency:
|
||||
|
||||
* Store high-throughput arrays (points, centroids) in global memory with
|
||||
coalesced access.
|
||||
|
||||
* Cache small, frequently reused values in shared memory or registers.
|
||||
|
||||
3. **Transfer scheduling**
|
||||
|
||||
Schedule host–device transfers asynchronously while the GPU executes kernels.
|
||||
Use HIP events or streams to synchronize only when necessary.
|
||||
|
||||
4. **Memory pooling and reuse**
|
||||
|
||||
Use memory pools such as :cpp:func:`hipMallocAsync`,
|
||||
:cpp:func:`hipMallocFromPoolAsync`, or custom allocators to mitigate
|
||||
fragmentation and reduce allocation latency, especially in iterative or
|
||||
batched workloads.
|
||||
|
||||
Performance Considerations
|
||||
--------------------------
|
||||
|
||||
.. list-table::
|
||||
:header-rows: 1
|
||||
:widths: 30 70
|
||||
|
||||
* - Aspect
|
||||
|
||||
- Recommendation
|
||||
|
||||
* - **Large datasets**
|
||||
|
||||
- GPU acceleration scales linearly with dataset size. Prioritize keeping
|
||||
data on GPU to avoid PCIe bottlenecks.
|
||||
|
||||
* - **Many iterations**
|
||||
|
||||
- Use persistent buffers and asynchronous reduction to minimize
|
||||
per-iteration synchronization overhead.
|
||||
|
||||
* - **Small number of clusters (k)**
|
||||
|
||||
- Kernel execution may be underutilized. CPU execution or hybrid scheduling
|
||||
may be more efficient.
|
||||
|
||||
* - **High-dimensional data**
|
||||
|
||||
- Distance computation dominates cost. Leverage GPU shared memory for
|
||||
centroid caching and ensure memory coalescing for point vectors.
|
||||
|
||||
When to use CPU-GPU cooperation
|
||||
===============================
|
||||
|
||||
* **Hybrid computational structure**
|
||||
|
||||
The algorithm exhibits a clear division between compute-intensive,
|
||||
data-parallel kernels such as distance computation in K-means and lightweight,
|
||||
control-heavy serial stages such as centroid updates or convergence checks.
|
||||
|
||||
* **High arithmetic intensity in parallel regions**
|
||||
|
||||
The GPU-executed portion performs substantial arithmetic operations per memory
|
||||
access, ensuring efficient utilization of GPU cores and minimizing the impact
|
||||
of memory latency.
|
||||
|
||||
* **Low-complexity reduction or synchronization on CPU**
|
||||
|
||||
The CPU phase primarily aggregates results or performs decision logic that
|
||||
does not justify GPU kernel execution overhead.
|
||||
|
||||
* **Large-scale data processing**
|
||||
|
||||
Dataset size exceeds the threshold where PCIe transfer costs can be amortized,
|
||||
enabling sustained GPU throughput and effective memory reuse.
|
||||
|
||||
* **Iterative or streaming workloads**
|
||||
|
||||
Algorithms involving multiple iterations benefit from persistent GPU memory
|
||||
allocations and overlapped data transfer, significantly reducing
|
||||
per-iteration latency.
|
||||
|
||||
When to avoid CPU-GPU cooperation
|
||||
=================================
|
||||
|
||||
* **Fully parallelizable workloads**
|
||||
|
||||
Algorithms with uniform, independent computations across all data points are
|
||||
better suited for GPU-only execution without CPU coordination overhead.
|
||||
|
||||
* **Low computational density per element**
|
||||
|
||||
If each data element requires few operations relative to transfer cost,
|
||||
CPU–GPU communication latency will dominate, leading to suboptimal
|
||||
performance.
|
||||
|
||||
* **High inter-phase data dependency**
|
||||
|
||||
Frequent synchronization or data exchange between CPU and GPU phases prevents
|
||||
effective pipeline overlap and leads to idle compute units.
|
||||
|
||||
* **Small problem sizes**
|
||||
|
||||
When datasets fit comfortably in CPU cache or system memory, the GPU launch
|
||||
overhead and transfer latency outweigh any computational gains from
|
||||
offloading.
|
||||
|
||||
Conclusion
|
||||
==========
|
||||
|
||||
The K-means implementation illustrates **heterogeneous workload partitioning**
|
||||
between CPU and GPU. By mapping compute-intensive, data-parallel operations to
|
||||
the GPU and reduction-heavy serial logic to the CPU, total runtime is minimized
|
||||
while code complexity remains manageable.
|
||||
|
||||
This CPU–GPU cooperative execution paradigm generalizes to many algorithms
|
||||
combining reduction, aggregation, and distance-based computation—enabling
|
||||
scalable, efficient utilization of modern heterogeneous hardware.
|
||||
@@ -0,0 +1,475 @@
|
||||
.. meta::
|
||||
:description: Matrix multiplication tutorial
|
||||
:keywords: AMD, ROCm, HIP, two-dimensional kernels, matrix multiplication tutorial
|
||||
|
||||
*******************************************************************************
|
||||
Two-dimensional kernels: Matrix multiplication tutorial
|
||||
*******************************************************************************
|
||||
|
||||
GPUs provide a massively parallel architecture consisting of thousands of cores,
|
||||
making them exceptionally well-suited for data-parallel computations.
|
||||
Two-dimensional kernel patterns are commonly data-parallel, enabling us to
|
||||
leverage GPU capabilities to exploit this inherent parallelism.
|
||||
|
||||
Tasks involving large matrices, which are common in image processing and machine
|
||||
learning applications, can be significantly accelerated by distributing
|
||||
computations across GPU cores. This tutorial explores how to implement matrix
|
||||
multiplication using two-dimensional GPU kernels with :doc:`HIP <hip:index>`.
|
||||
|
||||
.. include:: ../prerequisites.rst
|
||||
|
||||
Characteristics of 2D computational problems
|
||||
============================================
|
||||
|
||||
* Spatial locality and data dependencies: Adjacent elements in the grid often
|
||||
exhibit strong spatial correlations, making memory access patterns and cache
|
||||
utilization critical for performance.
|
||||
|
||||
* Natural 2D data representation: Many datasets—such as images, numerical
|
||||
matrices, and discretized physical fields—map directly onto a two-dimensional
|
||||
coordinate space.
|
||||
|
||||
* Prevalence in simulation and modeling: Numerous scientific and engineering
|
||||
workloads such as finite difference methods, fluid dynamics, heat transfer,
|
||||
and image processing, are inherently two-dimensional.
|
||||
|
||||
Modern GPU architectures are engineered to exploit the parallelism of 2D
|
||||
computational grids. By designing kernels that operate on two-dimensional thread
|
||||
blocks and memory layouts, developers can optimize global memory access,
|
||||
minimize latency, and maximize throughput. Leveraging 2D kernel configurations
|
||||
not only aligns the computation with the GPU’s hardware topology but also
|
||||
enables substantial performance improvements for domain-specific applications.
|
||||
|
||||
Matrix multiplication
|
||||
=====================
|
||||
|
||||
Let **A** and **B** be two matrices defined as follows,
|
||||
:math:`A \in \mathbb{R}^{m \times n}` and :math:`B \in \mathbb{R}^{n \times k}`.
|
||||
The matrix product :math:`C = A \cdot B` is defined only when the number of
|
||||
columns of :math:`A` equals the number of rows of :math:`B`. The resulting
|
||||
matrix :math:`C \in \mathbb{R}^{m \times k}` has elements given by:
|
||||
|
||||
.. math::
|
||||
|
||||
C_{ij} = \sum_{r=1}^{n} A_{ir} B_{rj}
|
||||
|
||||
for all :math:`i = 1, \dots, m` and :math:`j = 1, \dots, k`.
|
||||
|
||||
In other words, each element :math:`C_{ij}` is computed as the dot product of
|
||||
the :math:`i`-th row vector of :math:`A` and the *j*-th column vector of
|
||||
:math:`B`. This operation is repeated for all valid pairs of :math:`(i, j)` to
|
||||
construct the complete matrix :math:`C`.
|
||||
|
||||
For two square matrices of size :math:`N \times N`, the computational cost of
|
||||
classical matrix multiplication is :math:`O(N^3)`. As an example, consider
|
||||
:math:`N = 32`:
|
||||
|
||||
* **Multiplication operations**: :math:`32^3 = 32{,}768`
|
||||
* **Addition operations**: :math:`32^2 \times (32 - 1) = 31{,}744`
|
||||
|
||||
Each element in the resulting matrix :math:`C` is computed independently of the
|
||||
others, since it depends only on a single row of :math:`A` and a single column
|
||||
of :math:`B`. This property makes matrix multiplication
|
||||
**highly parallelizable** and well-suited for execution on GPUs, multi-core
|
||||
CPUs, or distributed computing architectures.
|
||||
|
||||
CPU implementation
|
||||
==================
|
||||
|
||||
A baseline CPU implementation provides a clear understanding of the classical
|
||||
matrix multiplication algorithm before exploring parallel GPU execution.
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
#include <iostream>
|
||||
#include <cstdlib>
|
||||
#define N 32
|
||||
|
||||
void cpu_matrix_multiplication(float *a, float *b, float *c, int n) {
|
||||
for (int i = 0; i < n; ++i) {
|
||||
for (int j = 0; j < n; ++j) {
|
||||
float sum = 0.0f;
|
||||
for (int k = 0; k < n; ++k) {
|
||||
sum += a[i * n + k] * b[k * n + j];
|
||||
}
|
||||
c[i * n + j] = sum;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int main() {
|
||||
float *a, *b, *c;
|
||||
a = (float*)malloc(sizeof(float) * N * N);
|
||||
b = (float*)malloc(sizeof(float) * N * N);
|
||||
c = (float*)malloc(sizeof(float) * N * N);
|
||||
|
||||
// Initialize matrix A
|
||||
for (int i = 0; i < N; ++i) {
|
||||
for (int j = 0; j < N; ++j) {
|
||||
a[i * N + j] = static_cast<float>(rand()) / RAND_MAX;
|
||||
}
|
||||
}
|
||||
|
||||
// Initialize matrix B
|
||||
for (int i = 0; i < N; ++i) {
|
||||
for (int j = 0; j < N; ++j) {
|
||||
b[i * N + j] = static_cast<float>(rand()) / RAND_MAX;
|
||||
}
|
||||
}
|
||||
|
||||
cpu_matrix_multiplication(a, b, c, N);
|
||||
free(a);
|
||||
free(b);
|
||||
free(c);
|
||||
return 0;
|
||||
}
|
||||
|
||||
The ``cpu_matrix_multiplication`` function performs the classical :math:`O(N^3)`
|
||||
matrix multiplication algorithm using three nested loops. The implementation
|
||||
proceeds as follows:
|
||||
|
||||
* **Input parameters:** Three pointers to contiguous memory blocks representing
|
||||
matrices :math:`A`, :math:`B`, and the output matrix :math:`C`.
|
||||
|
||||
* **Outer and middle loops:** The indices :math:`i` and :math:`j` iterate over
|
||||
the rows and columns of the output matrix :math:`C`, respectively.
|
||||
|
||||
* **Innermost loop:** For each element :math:`C_{ij}`, the loop over :math:`k`
|
||||
performs a dot product between the *i*-th row of :math:`A` and the :math:`j`-th
|
||||
column of :math:`B`:
|
||||
|
||||
.. math::
|
||||
|
||||
C_{ij} = \sum_{k=0}^{n-1} A_{ik} \cdot B_{kj}
|
||||
|
||||
* **Temporary accumulation:** A local scalar :code:`sum` accumulates the
|
||||
intermediate sum before being written to :code:`c[i * n + j]`.
|
||||
|
||||
This implementation has a computational complexity of :math:`O(N^3)` and poor
|
||||
cache locality for large matrices, but it serves as a reference for
|
||||
understanding sequential computation before introducing GPU parallelization.
|
||||
|
||||
GPU implementation
|
||||
==================
|
||||
|
||||
The following example demonstrates a complete HIP implementation of matrix
|
||||
multiplication, including host and device memory management, kernel
|
||||
implementation, configuration, and synchronization.
|
||||
|
||||
GPU kernel
|
||||
----------
|
||||
|
||||
The core computation is performed by the GPU kernel, where each thread computes
|
||||
one element of the output matrix. For comparison, the CPU implementation is
|
||||
also provided.
|
||||
|
||||
.. list-table::
|
||||
:header-rows: 1
|
||||
|
||||
* - **GPU version**
|
||||
|
||||
* - .. code-block:: c++
|
||||
|
||||
__global__ void gpu_matrix_multiplication(float *a, float *b, float *c, int n) {
|
||||
|
||||
int row = blockIdx.y * blockDim.y + threadIdx.y;
|
||||
int col = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
float sum = 0.0f;
|
||||
|
||||
if (row < n && col < n) {
|
||||
for (int k = 0; k < n; ++k) {
|
||||
sum += a[row * n + k] * b[k * n + col];
|
||||
}
|
||||
c[row * n + col] = sum;
|
||||
}
|
||||
}
|
||||
|
||||
.. list-table::
|
||||
:header-rows: 1
|
||||
|
||||
* - **CPU version**
|
||||
|
||||
* - .. code-block:: c++
|
||||
|
||||
void cpu_matrix_multiplication(float *a, float *b, float *c, int n) {
|
||||
|
||||
for (int i = 0; i < n; ++i) {
|
||||
for (int j = 0; j < n; ++j) {
|
||||
float sum = 0.0f;
|
||||
for (int k = 0; k < n; ++k) {
|
||||
sum += a[i * n + k] * b[k * n + j];
|
||||
}
|
||||
c[i * n + j] = sum;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
The outer and middle loops of the CPU implementation are replaced by
|
||||
the parallel execution of the GPU implementation. Each GPU thread computes the
|
||||
**single element** :math:`C_{ij}` of the output matrix corresponding to one dot
|
||||
product between a row of :math:`A` and a column of :math:`B`. This decomposition
|
||||
exposes massive parallelism, as all elements of :math:`C` can be computed
|
||||
independently and concurrently.
|
||||
|
||||
Thread and block identification
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
Each thread’s global position within the grid is determined by:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
int row = blockIdx.y * blockDim.y + threadIdx.y;
|
||||
int col = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
|
||||
Here:
|
||||
|
||||
* ``threadIdx.(x|y)``: Local thread indices within a block.
|
||||
* ``blockIdx.(x|y)``: Block indices within the grid.
|
||||
* ``blockDim.(x|y)``: Dimensions of each block (used to scale offsets).
|
||||
|
||||
Boundary checking
|
||||
~~~~~~~~~~~~~~~~~
|
||||
|
||||
Since the total number of threads launched may exceed :math:`N^2`, boundary
|
||||
checking ensures that threads outside the matrix domain do not perform invalid
|
||||
memory accesses:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
if (row < n && col < n) {
|
||||
// Safe computation region
|
||||
}
|
||||
|
||||
Dot product computation
|
||||
~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
Within the valid region, each thread executes a dot product over :math:`k`:
|
||||
|
||||
* Loads one element from row :code:`row` of matrix :math:`A`.
|
||||
|
||||
* Loads one element from column :code:`col` of matrix :math:`B`.
|
||||
|
||||
* Multiplies and accumulates these values into :code:`sum`.
|
||||
|
||||
* Writes the final scalar result to :code:`c[row * n + col]`.
|
||||
|
||||
This kernel performs the same :math:`O(N^3)` arithmetic operations as the CPU
|
||||
version but distributes them across thousands of concurrent GPU threads,
|
||||
achieving significant acceleration through parallel execution and memory
|
||||
throughput optimization.
|
||||
|
||||
Step 1: Host memory allocation and initialization
|
||||
--------------------------------------------------
|
||||
|
||||
Host memory is allocated for matrices and initialized with random floating-point
|
||||
values.
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
int main() {
|
||||
float *h_a, *h_b, *h_c;
|
||||
h_a = (float*)malloc(sizeof(float) * N * N);
|
||||
h_b = (float*)malloc(sizeof(float) * N * N);
|
||||
h_c = (float*)malloc(sizeof(float) * N * N);
|
||||
|
||||
// Initialize matrix A
|
||||
for (int i = 0; i < N; ++i) {
|
||||
for (int j = 0; j < N; ++j) {
|
||||
h_a[i * N + j] = static_cast<float>(rand()) / RAND_MAX;
|
||||
}
|
||||
}
|
||||
|
||||
// Initialize matrix B
|
||||
for (int i = 0; i < N; ++i) {
|
||||
for (int j = 0; j < N; ++j) {
|
||||
h_b[i * N + j] = static_cast<float>(rand()) / RAND_MAX;
|
||||
}
|
||||
}
|
||||
|
||||
**Description:**
|
||||
|
||||
* Declare host (CPU) pointers ``h_a``, ``h_b``, and ``h_c``.
|
||||
|
||||
* Allocate contiguous memory for each :math:`N \times N` matrix.
|
||||
|
||||
* Initialize input matrices :math:`A` and :math:`B` with pseudo-random
|
||||
floating-point values in the range [0, 1).
|
||||
|
||||
Step 2: Device memory allocation and data transfer
|
||||
---------------------------------------------------
|
||||
|
||||
Memory is allocated on the GPU, and input matrices are transferred from host to
|
||||
device.
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
float *d_a, *d_b, *d_c;
|
||||
hipMalloc((void**)&d_a, sizeof(float) * N * N);
|
||||
hipMalloc((void**)&d_b, sizeof(float) * N * N);
|
||||
hipMalloc((void**)&d_c, sizeof(float) * N * N);
|
||||
|
||||
hipMemcpy(d_a, h_a, sizeof(float) * N * N, hipMemcpyHostToDevice);
|
||||
hipMemcpy(d_b, h_b, sizeof(float) * N * N, hipMemcpyHostToDevice);
|
||||
|
||||
**Operations:**
|
||||
|
||||
1. Allocate GPU (device) memory for matrices ``d_a``, ``d_b``, and ``d_c``.
|
||||
2. Transfer data from host to device using :cpp:func:`hipMemcpy` with direction
|
||||
:cpp:enum:`hipMemcpyHostToDevice`.
|
||||
|
||||
Step 3: Configure and launch kernel
|
||||
------------------------------------
|
||||
|
||||
Kernel launch parameters define how threads are organized across blocks and the
|
||||
grid.
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
dim3 threadsPerBlock(BLOCK_SIZE, BLOCK_SIZE);
|
||||
int n_blocks = static_cast<int>(ceil(static_cast<float>(N) / BLOCK_SIZE));
|
||||
dim3 blocksPerGrid(n_blocks, n_blocks);
|
||||
|
||||
hipLaunchKernelGGL(gpu_matrix_multiplication,
|
||||
blocksPerGrid,
|
||||
threadsPerBlock,
|
||||
0, 0,
|
||||
d_a, d_b, d_c, N);
|
||||
|
||||
hipDeviceSynchronize();
|
||||
|
||||
Configuration details
|
||||
~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
The :code:`dim3` type defines thread and block dimensions:
|
||||
|
||||
* :code:`threadsPerBlock`: Number of threads per block. A 16 × 16 block
|
||||
(256 threads total).
|
||||
|
||||
* :code:`n_blocks`: Number of blocks per dimension. Computed as
|
||||
:math:`\lceil N / \mathrm{BLOCK\_SIZE} \rceil`.
|
||||
|
||||
* :code:`blocksPerGrid`: A grid of blocks covering the entire :math:`N \times N`
|
||||
matrix.
|
||||
|
||||
For example :math:`N = 256` and ``BLOCK_SIZE = 16``:
|
||||
|
||||
* :code:`n_blocks`: :math:`\lceil 256 / 16 \rceil = 16`
|
||||
|
||||
* :code:`blocksPerGrid`: :math:`16 \times 16 = 256`
|
||||
|
||||
* Total threads: :math:`256 \text{ blocks} \times 256 \text{ threads/block} = 65{,}536`
|
||||
threads.
|
||||
|
||||
Rounding up ensures full coverage of the matrix even when :math:`N` is not an
|
||||
exact multiple of ``BLOCK_SIZE``. The boundary check in the kernel
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
if (row < n && col < n)
|
||||
|
||||
prevents out-of-bounds memory access for extra threads.
|
||||
|
||||
Synchronization
|
||||
~~~~~~~~~~~~~~~~
|
||||
|
||||
The call to :cpp:func:`hipDeviceSynchronize()` ensures that all GPU computations
|
||||
complete before the CPU accesses results or proceeds to subsequent operations.
|
||||
This is essential for correctness and debugging.
|
||||
|
||||
Step 4: Copy results back and cleanup
|
||||
--------------------------------------
|
||||
|
||||
After the kernel execution, results are transferred back to host memory, and
|
||||
all allocated resources are released.
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
hipMemcpy(h_c, d_c, sizeof(float) * N * N, hipMemcpyDeviceToHost);
|
||||
|
||||
hipFree(d_a);
|
||||
hipFree(d_b);
|
||||
hipFree(d_c);
|
||||
free(h_a);
|
||||
free(h_b);
|
||||
free(h_c);
|
||||
return 0;
|
||||
}
|
||||
|
||||
**Summary of final steps:**
|
||||
|
||||
1. Copy the result matrix from device to host using ``hipMemcpy``.
|
||||
2. Free all device memory allocations with ``hipFree``.
|
||||
3. Release host memory with ``free``.
|
||||
4. Return control to the operating system, indicating successful program termination.
|
||||
|
||||
Parallelization benefits
|
||||
========================
|
||||
|
||||
For a 256 × 256 matrix multiplication:
|
||||
|
||||
* **Sequential CPU version:** Computes 65,536 output elements serially.
|
||||
|
||||
* **Parallel GPU version:** Executes up to 65,536 independent threads concurrently.
|
||||
|
||||
This results in a theoretical performance gain proportional to the number of
|
||||
active GPU threads and the device’s compute throughput.
|
||||
|
||||
Each output element :math:`C_{ij}` is computed independently from others, since
|
||||
it depends solely on row :math:`i` of matrix :math:`A` and column :math:`j` of
|
||||
matrix :math:`B`. This independence allows full utilization of GPU streaming
|
||||
multiprocessors and makes the algorithm highly scalable.
|
||||
|
||||
Best practices
|
||||
==============
|
||||
|
||||
1. **Choose optimal block sizes**
|
||||
|
||||
Powers of two (e.g., 16 or 32) often yield better occupancy and memory
|
||||
alignment.
|
||||
|
||||
2. **Handle boundary conditions**
|
||||
|
||||
Always include thread boundary checks.
|
||||
|
||||
3. **Synchronize appropriately**
|
||||
|
||||
Use :cpp:func:`hipDeviceSynchronize()` after kernel launches to ensure data
|
||||
consistency.
|
||||
|
||||
4. **Memory coalescing**
|
||||
|
||||
Arrange data access patterns so consecutive threads access contiguous memory
|
||||
locations, maximizing bandwidth utilization.
|
||||
|
||||
5. **Use shared memory**
|
||||
|
||||
Use shared memory to cache sub-blocks of matrices, significantly
|
||||
reducing global memory latency.
|
||||
|
||||
6. **Profile and tune**
|
||||
|
||||
Use tools such as :doc:`rocprofv3<rocprofiler-sdk:how-to/using-rocprofv3>`
|
||||
or :doc:`ROCm compute profiler<rocprofiler-compute:how-to/profile/mode>`
|
||||
to identify bottlenecks and fine-tune kernel launch configurations.
|
||||
|
||||
Conclusion
|
||||
==========
|
||||
|
||||
Two-dimensional GPU kernels provide an efficient mechanism to accelerate dense
|
||||
linear algebra computations such as matrix multiplication by exploiting
|
||||
fine-grained data parallelism. This example demonstrates:
|
||||
|
||||
* Structuring GPU kernels for 2D problems.
|
||||
|
||||
* Managing memory transfers between host and device.
|
||||
|
||||
* Configuring thread and block hierarchies.
|
||||
|
||||
* Achieving substantial speedups via massive parallel execution.
|
||||
|
||||
Understanding these concepts enables developers to implement optimized GPU
|
||||
solutions for computationally intensive workloads, including scientific
|
||||
simulations, numerical linear algebra, and machine learning. Because each output
|
||||
element is computed independently, matrix multiplication serves as an ideal
|
||||
introductory example for mastering GPU programming paradigms applicable to a
|
||||
wide range of data-parallel applications.
|
||||
@@ -0,0 +1,484 @@
|
||||
.. meta::
|
||||
:description: Multi-kernel programming with breadth-first search (BFS)
|
||||
:keywords: AMD, ROCm, HIP, GPU programming, multi-kernel, BFS, breadth-first search
|
||||
|
||||
*******************************************************************************
|
||||
Multi-kernel programming: breadth-first search tutorial
|
||||
*******************************************************************************
|
||||
|
||||
Many real-world GPU workloads involve multiple kernels cooperating to solve a
|
||||
single problem. This tutorial explores **multi-kernel GPU programming** using
|
||||
the breadth-first search (BFS) algorithm, a foundational graph traversal
|
||||
method widely used in networking, path-finding, and social network analysis.
|
||||
|
||||
The implementation is adapted from the **Rodinia benchmark suite**, a
|
||||
well-known collection of heterogeneous computing workloads that demonstrate
|
||||
different parallel programming strategies.
|
||||
|
||||
.. include:: ../prerequisites.rst
|
||||
|
||||
Multi-kernel GPU programming
|
||||
============================
|
||||
|
||||
In GPU computing, some algorithms cannot be efficiently expressed using a
|
||||
single kernel due to synchronization or dependency constraints. Instead, they
|
||||
are decomposed into multiple kernels that execute sequentially, with each
|
||||
kernel responsible for a specific computation phase.
|
||||
|
||||
This approach, called **multi-kernel programming**, is essential when:
|
||||
|
||||
* Results from one kernel determine the input for the next.
|
||||
|
||||
* Global synchronization between thread blocks is required.
|
||||
|
||||
* Control flow depends on runtime conditions.
|
||||
|
||||
* The algorithm involves iterative or level-wise processing.
|
||||
|
||||
Breadth-first search (BFS)
|
||||
==========================
|
||||
|
||||
Breadth-first search (BFS) is a **layered graph traversal algorithm** that
|
||||
explores nodes level by level, starting from a root node. It guarantees finding
|
||||
the shortest path (in edge count) to all reachable nodes in an unweighted
|
||||
graph.
|
||||
|
||||
Applications of BFS include:
|
||||
|
||||
* **Path-finding**: Finding shortest paths between nodes.
|
||||
* **Peer-to-peer networking**: Network topology discovery.
|
||||
* **GPS navigation**: Route planning and optimization.
|
||||
* **Social networks**: Friend recommendations and connection analysis.
|
||||
* **Web crawling**: Systematic website exploration.
|
||||
|
||||
Algorithm characteristics
|
||||
-------------------------
|
||||
|
||||
BFS is structured as a level-synchronous algorithm:
|
||||
|
||||
* Nodes in the same graph level are processed concurrently.
|
||||
* A queue (or "frontier") tracks which nodes to explore next.
|
||||
* Each node is visited once to prevent redundant processing.
|
||||
|
||||
Sequential BFS is straightforward but inherently serial due to the
|
||||
level-by-level dependency between nodes. GPU parallelization requires
|
||||
restructuring the traversal to exploit data parallelism across nodes
|
||||
within the same frontier.
|
||||
|
||||
Sequential BFS algorithm
|
||||
=========================
|
||||
|
||||
Let's first understand how BFS works sequentially before parallelizing it.
|
||||
|
||||
Example Graph
|
||||
~~~~~~~~~~~~~
|
||||
|
||||
Consider a simple graph with four nodes:
|
||||
|
||||
.. code-block:: text
|
||||
|
||||
R (root)
|
||||
/ \
|
||||
A B
|
||||
\ /
|
||||
C
|
||||
|
||||
Step-by-step execution
|
||||
~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
**Step 1**: Start at the root node ``R``
|
||||
|
||||
* Mark ``R`` as visited
|
||||
* Enqueue ``R``
|
||||
* Queue: [R]
|
||||
|
||||
**Step 2**: Process ``R``
|
||||
|
||||
* Dequeue ``R``
|
||||
* Discover neighbors: ``A`` and ``B``
|
||||
* Enqueue both, mark as visited
|
||||
* Queue: [A, B]
|
||||
|
||||
**Step 3**: Process ``A``
|
||||
|
||||
* Dequeue ``A``
|
||||
* Neighbors: ``R`` (visited) and ``C`` (new)
|
||||
* Enqueue ``C``
|
||||
* Queue: [B, C]
|
||||
|
||||
**Step 4**: Process ``B``
|
||||
|
||||
* Dequeue ``B``
|
||||
* Neighbors: ``R`` (visited) and ``C`` (visited)
|
||||
* Queue: [C]
|
||||
|
||||
**Step 5**: Process ``C``
|
||||
|
||||
* Dequeue ``C``
|
||||
* All neighbors visited
|
||||
* Queue becomes empty — traversal complete
|
||||
|
||||
Parallel BFS on GPU
|
||||
===================
|
||||
|
||||
Unlike dense linear algebra, BFS is an **irregular** algorithm. The amount of
|
||||
work per node varies, and the connectivity pattern of the graph drives
|
||||
execution. The main challenges are:
|
||||
|
||||
1. **Data dependencies**: nodes in the next level depend on the previous level.
|
||||
|
||||
2. **Irregular parallelism**: each frontier may contain a very different number of nodes.
|
||||
|
||||
3. **Dynamic workload**: the size of the next frontier is unknown at runtime.
|
||||
|
||||
4. **Synchronization**: all nodes in one frontier must complete before the next begins.
|
||||
|
||||
The **frontier** is the set of nodes being processed at a given BFS level.
|
||||
Parallel BFS executes all frontier nodes simultaneously, using one thread per
|
||||
node to discover new neighbors and mark them for the next iteration.
|
||||
|
||||
Implementation strategy
|
||||
-----------------------
|
||||
|
||||
The GPU implementation performs BFS using **two cooperating kernels**:
|
||||
|
||||
1. **Kernel 1**: processes all nodes in the current frontier.
|
||||
2. **Kernel 2**: updates the next frontier and checks if work remains.
|
||||
|
||||
This design provides **implicit synchronization** between levels while avoiding
|
||||
race conditions. The host (CPU) manages the iterative control loop, launching
|
||||
kernels repeatedly until no more frontier nodes exist.
|
||||
|
||||
Data structures
|
||||
===============
|
||||
|
||||
The graph is represented using adjacency lists stored in arrays:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
struct Node {
|
||||
int starting; // starting index in the edge list
|
||||
int no_of_edges; // number of outgoing edges
|
||||
};
|
||||
|
||||
**Main arrays:**
|
||||
|
||||
* ``g_graph_nodes``: node array storing offsets into the edge list.
|
||||
* ``g_graph_edges``: flattened list of edge destinations.
|
||||
* ``g_graph_mask``: boolean array indicating active frontier nodes.
|
||||
* ``g_updating_graph_mask``: marks nodes to be added to the next frontier.
|
||||
* ``g_graph_visited``: tracks which nodes were visited.
|
||||
* ``g_graph_cost``: stores the distance (edge count) from the source node.
|
||||
|
||||
**Control flow flags:**
|
||||
|
||||
* ``g_over``: device-side flag indicating whether another iteration is needed.
|
||||
* The host resets this flag each iteration and checks it after kernel execution.
|
||||
|
||||
The two-kernel approach
|
||||
=======================
|
||||
|
||||
The two-kernel structure ensures correctness and efficient synchronization:
|
||||
|
||||
* **Exploration kernel (Kernel 1)** discovers new nodes.
|
||||
* **Update kernel (Kernel 2)** finalizes state for the next iteration.
|
||||
|
||||
This separation:
|
||||
|
||||
* Avoids race conditions between threads of different levels.
|
||||
|
||||
* Provides synchronization between BFS levels.
|
||||
|
||||
* Keeps control logic simple on the host side.
|
||||
|
||||
Kernel 1: process current frontier
|
||||
----------------------------------
|
||||
|
||||
Each thread processes one node from the current frontier, examining all of its
|
||||
outgoing edges:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
__global__ void Kernel1(
|
||||
Node* g_graph_nodes,
|
||||
int* g_graph_edges,
|
||||
bool* g_graph_mask,
|
||||
bool* g_updating_graph_mask,
|
||||
bool* g_graph_visited,
|
||||
int* g_graph_cost,
|
||||
int no_of_nodes)
|
||||
{
|
||||
int tid = hipBlockIdx_x * MAX_THREADS_PER_BLOCK + hipThreadIdx_x;
|
||||
|
||||
if (tid < no_of_nodes && g_graph_mask[tid]) {
|
||||
g_graph_mask[tid] = false;
|
||||
|
||||
for (int i = g_graph_nodes[tid].starting;
|
||||
i < g_graph_nodes[tid].starting + g_graph_nodes[tid].no_of_edges;
|
||||
i++) {
|
||||
int id = g_graph_edges[i];
|
||||
if (!g_graph_visited[id]) {
|
||||
g_graph_cost[id] = g_graph_cost[tid] + 1;
|
||||
g_updating_graph_mask[id] = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
**Kernel 1 responsibilities:**
|
||||
|
||||
* Clear the node’s mask (mark processed).
|
||||
|
||||
* Explore all edges.
|
||||
|
||||
* For each unvisited neighbor:
|
||||
|
||||
* Compute cost (distance).
|
||||
|
||||
* Add to the next frontier.
|
||||
|
||||
Kernel 2: update frontier
|
||||
-------------------------
|
||||
|
||||
This kernel finalizes the next frontier:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
__global__ void Kernel2(
|
||||
bool* g_graph_mask,
|
||||
bool* g_updating_graph_mask,
|
||||
bool* g_graph_visited,
|
||||
bool* g_over,
|
||||
int no_of_nodes)
|
||||
{
|
||||
int tid = hipBlockIdx_x * MAX_THREADS_PER_BLOCK + hipThreadIdx_x;
|
||||
|
||||
if (tid < no_of_nodes && g_updating_graph_mask[tid]) {
|
||||
g_graph_mask[tid] = true;
|
||||
g_graph_visited[tid] = true;
|
||||
*g_over = true;
|
||||
g_updating_graph_mask[tid] = false;
|
||||
}
|
||||
}
|
||||
|
||||
**Kernel 2 responsibilities:**
|
||||
|
||||
* Move newly discovered nodes into the active frontier.
|
||||
|
||||
* Mark them as visited.
|
||||
|
||||
* Signal continuation via ``*g_over``.
|
||||
|
||||
Host-side control loop
|
||||
======================
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
do {
|
||||
h_over = false;
|
||||
hipMemcpy(d_over, &h_over, sizeof(bool), hipMemcpyHostToDevice);
|
||||
|
||||
Kernel1<<<num_blocks, MAX_THREADS_PER_BLOCK>>>(
|
||||
d_graph_nodes, d_graph_edges, d_graph_mask,
|
||||
d_graph_updating_graph_mask, d_graph_visited,
|
||||
d_graph_cost, no_of_nodes);
|
||||
hipDeviceSynchronize();
|
||||
|
||||
Kernel2<<<num_blocks, MAX_THREADS_PER_BLOCK>>>(
|
||||
d_graph_mask, d_graph_updating_graph_mask,
|
||||
d_graph_visited, d_over, no_of_nodes);
|
||||
hipDeviceSynchronize();
|
||||
|
||||
hipMemcpy(&h_over, d_over, sizeof(bool), hipMemcpyDeviceToHost);
|
||||
} while (h_over);
|
||||
|
||||
|
||||
The loop exits when no new nodes are discovered. ``g_over`` or ``h_over`` on
|
||||
host side remains ``false`` after one full iteration.
|
||||
|
||||
Performance Characteristics
|
||||
===========================
|
||||
|
||||
Parallelism Patterns
|
||||
--------------------
|
||||
|
||||
**Within each iteration:**
|
||||
|
||||
- High parallelism: All frontier nodes processed simultaneously
|
||||
- Work distribution: One thread per node
|
||||
|
||||
**Across iterations:**
|
||||
|
||||
- Sequential: Must complete one level before starting the next
|
||||
- Variable parallelism: Different levels may have different numbers of nodes
|
||||
|
||||
Workload Characteristics
|
||||
------------------------
|
||||
|
||||
.. list-table::
|
||||
:header-rows: 1
|
||||
:widths: 30 70
|
||||
|
||||
* - Characteristic
|
||||
- Description
|
||||
* - **Irregular**
|
||||
- Frontier size varies dramatically across levels
|
||||
* - **Data-dependent**
|
||||
- Graph structure determines parallel work available
|
||||
* - **Dynamic**
|
||||
- Cannot predict workload statically
|
||||
* - **Memory-bound**
|
||||
- Many memory accesses per computation
|
||||
|
||||
Best practices
|
||||
==============
|
||||
|
||||
This section outlines recommended practices for implementing an efficient
|
||||
GPU-accelerated Breadth-First Search (BFS). It highlights design principles,
|
||||
memory-management strategies, and debugging techniques that help ensure
|
||||
correctness, maintainability, and high performance when mapping BFS onto modern
|
||||
GPU architectures.
|
||||
|
||||
Design principles
|
||||
-----------------
|
||||
|
||||
1. **Define clear kernel roles**
|
||||
|
||||
Decompose BFS into well-defined GPU kernels, each responsible for a specific
|
||||
phase of computation. For example:
|
||||
|
||||
* **Kernel 1**: frontier expansion (discovering new nodes)
|
||||
* **Kernel 2**: frontier update (marking next-level nodes)
|
||||
|
||||
This separation simplifies synchronization and ensures that each kernel
|
||||
operates on independent data regions.
|
||||
|
||||
2. **Minimize host–device communication**
|
||||
|
||||
Keep graph data structures (nodes, edges, masks) resident on the GPU across
|
||||
iterations. Only transfer lightweight control flags such as ``g_over`` to the
|
||||
host each loop iteration to check termination conditions.
|
||||
|
||||
3. **Kernel boundaries as synchronization points**
|
||||
|
||||
Kernel launch boundaries on the same stream naturally enforce global
|
||||
synchronization across all threads on the GPU. Each kernel invocation
|
||||
completes before the next begins, ensuring that:
|
||||
|
||||
* All nodes in the current frontier are fully processed before updating the
|
||||
next frontier.
|
||||
|
||||
* Memory updates to arrays like ``g_graph_cost`` or ``g_graph_mask`` are
|
||||
visible to all threads in subsequent kernels.
|
||||
|
||||
This avoids the need for costly device-wide barriers or explicit
|
||||
synchronization primitives within a single kernel. Leverage kernel sequencing
|
||||
to structure iterative algorithms cleanly—each kernel represents one
|
||||
computation phase per BFS level.
|
||||
|
||||
4. **Flag-based control**
|
||||
|
||||
Use device-side flags for dynamic termination and conditional control flow.
|
||||
In BFS, the Boolean flag ``g_over`` serves as a device-to-host signal
|
||||
indicating whether new nodes were discovered during the current iteration.
|
||||
|
||||
* Initialize ``g_over`` to ``false`` on the host at the start of each
|
||||
iteration.
|
||||
|
||||
* Allow GPU threads in **Kernel 2** to set ``*g_over = true`` when adding new
|
||||
nodes to the next frontier.
|
||||
|
||||
* After kernel completion, copy the flag back to the host using
|
||||
:cpp:func:`hipMemcpy`. If ``g_over`` remains false, the traversal is
|
||||
complete.
|
||||
|
||||
This mechanism avoids repeated host intervention and enables a tight
|
||||
CPU–GPU control loop that dynamically adapts to workload size without
|
||||
transferring large data structures.
|
||||
|
||||
Memory strategy
|
||||
---------------
|
||||
|
||||
1. **Persistent device allocations**
|
||||
|
||||
Allocate all required device buffers once prior to traversal. Reuse these
|
||||
allocations across multiple BFS runs or multiple source nodes to minimize
|
||||
the overhead of repeated :cpp:func:`hipMalloc` and :cpp:func:`hipFree`
|
||||
calls.
|
||||
|
||||
2. **Minimize host–device communication**
|
||||
|
||||
Keep graph data structures (nodes, edges, masks) resident on the GPU across
|
||||
iterations. Only transfer lightweight control flags such as ``g_over`` to the
|
||||
host each loop iteration to check termination conditions.
|
||||
|
||||
3. **Use pinned host memory for control flags**
|
||||
|
||||
When copying ``g_over`` or other control signals between host and device,
|
||||
allocate host memory using pinned (page-locked) buffers to accelerate DMA
|
||||
transfers.
|
||||
|
||||
Debugging and validation
|
||||
------------------------
|
||||
|
||||
1. **Frontier validation**
|
||||
|
||||
After each iteration, verify the number of nodes marked in
|
||||
``g_graph_mask``. Unexpected empty or overfull frontiers often indicate
|
||||
incorrect synchronization or uninitialized masks.
|
||||
|
||||
2. **Termination condition check**
|
||||
|
||||
Confirm that the host-side loop terminates when ``g_over`` remains false
|
||||
for one iteration. If the loop never ends, ensure ``g_over`` is reset on the
|
||||
host before each kernel launch.
|
||||
|
||||
3. **Result verification**
|
||||
|
||||
Compare computed distances in ``g_graph_cost`` against a CPU reference
|
||||
implementation for small graphs to validate correctness.
|
||||
|
||||
4. **Profiling and bottleneck detection**
|
||||
|
||||
Use tools such as :doc:`rocprofv3<rocprofiler-sdk:how-to/using-rocprofv3>`
|
||||
or :doc:`ROCm compute profiler<rocprofiler-compute:how-to/profile/mode>`
|
||||
to measure per-kernel execution times, memory throughput, and
|
||||
synchronization overhead.
|
||||
|
||||
5. **Logging and debug builds**
|
||||
|
||||
Enable optional logging for iteration counts, frontier sizes, and
|
||||
synchronization states during development. Disable logging in production
|
||||
builds to avoid performance impact.
|
||||
|
||||
Conclusion
|
||||
==========
|
||||
|
||||
Multi-kernel GPU programming is essential for complex algorithms that require:
|
||||
|
||||
* Multiple phases of computation.
|
||||
|
||||
* Data dependencies between phases.
|
||||
|
||||
* Dynamic control flow based on intermediate results.
|
||||
|
||||
The BFS example demonstrates:
|
||||
|
||||
* How to decompose algorithms into multiple cooperating kernels.
|
||||
|
||||
* Techniques for managing frontiers and iterative processing.
|
||||
|
||||
* Strategies for handling irregular and dynamic parallelism.
|
||||
|
||||
* Proper synchronization between kernel launches.
|
||||
|
||||
Key takeaways:
|
||||
|
||||
1. **Kernel boundaries provide synchronization**: Use them strategically to ensure correctness.
|
||||
2. **Separate exploration from update**: Prevents race conditions in level-based algorithms.
|
||||
3. **Host controls iteration**: CPU manages the overall loop while GPU does heavy lifting.
|
||||
4. **Flags enable dynamic control**: Device-side flags allow work-dependent termination.
|
||||
|
||||
Understanding multi-kernel patterns enables developers to implement
|
||||
sophisticated algorithms like graph processing, dynamic programming, and
|
||||
iterative refinement methods efficiently on GPUs.
|
||||
@@ -0,0 +1,444 @@
|
||||
.. meta::
|
||||
:description: Image convolution tutorial using HIP on AMD GPUs
|
||||
:keywords: AMD, ROCm, HIP, stencil operations, image convolution, GPU programming, data parallelism
|
||||
|
||||
*******************************************************************************
|
||||
Stencil operations: Image convolution tutorial
|
||||
*******************************************************************************
|
||||
|
||||
Stencil operations represent an important class of **embarrassingly parallel**
|
||||
algorithms that are ideally suited for GPU acceleration. A stencil algorithm
|
||||
iteratively updates data in an array based on a data item's adjacent cells,
|
||||
making it a fundamental technique in various computational domains.
|
||||
|
||||
.. include:: ../prerequisites.rst
|
||||
|
||||
Applications of stencil operations
|
||||
===================================
|
||||
|
||||
Stencil algorithms are commonly used in:
|
||||
|
||||
* **Physics simulations**: Modeling heat transfer, fluid dynamics, and wave
|
||||
propagation
|
||||
|
||||
* **Partial differential equations**: Numerical solutions to scientific
|
||||
computing problems
|
||||
|
||||
* **Image processing**: Convolutional operations for smoothing, sharpening, and
|
||||
edge detection
|
||||
|
||||
* **Convolutional Neural Networks**: A major building block in modern
|
||||
deep learning
|
||||
|
||||
By applying different image convolution kernels (not to be confused with GPU
|
||||
kernels), these algorithms can smooth and sharpen image features and detect
|
||||
edges effectively.
|
||||
|
||||
Image convolution
|
||||
=================
|
||||
|
||||
An image convolution applies a small matrix (the **mask** or **filter kernel**)
|
||||
to an input image. For each pixel :math:`(x, y)`, the output is computed as:
|
||||
|
||||
.. math::
|
||||
|
||||
I'(x, y) = \sum_{i=-r}^{r} \sum_{j=-r}^{r} M(i, j) \cdot I(x + i, y + j)
|
||||
|
||||
where :math:`M(i, j)` is the mask coefficient, and :math:`r` is half the mask
|
||||
width (assuming a square kernel).
|
||||
|
||||
Step by step description of the equation:
|
||||
|
||||
1. Center the mask over the current pixel
|
||||
2. Multiply each mask value by the corresponding image pixel
|
||||
3. Sum all the products
|
||||
4. Store the result as the new pixel value
|
||||
|
||||
Dimensionality of stencils
|
||||
--------------------------
|
||||
|
||||
Stencil computations extend beyond 2D image grids:
|
||||
|
||||
* **1D:** Signal filtering, time series processing
|
||||
|
||||
* **2D:** Image processing, texture analysis
|
||||
|
||||
* **3D:** Volume data, fluid flow, and physical field simulation
|
||||
|
||||
This tutorial focuses on **2D image convolution**, the most common stencil
|
||||
operation in visual and scientific computing.
|
||||
|
||||
The smoothing operation
|
||||
-----------------------
|
||||
|
||||
The tutorial implements a **box blur** (uniform smoothing filter). Each
|
||||
pixel’s new value is the average of its local neighborhood. This operation:
|
||||
|
||||
* Reduces noise by averaging local intensity variations.
|
||||
|
||||
* Acts as a low-pass filter, attenuating high-frequency components.
|
||||
|
||||
* Provides an ideal example of a stencil computation with uniform weights.
|
||||
|
||||
Two-dimensional grid architecture
|
||||
==================================
|
||||
|
||||
The tutorial uses a two-dimensional grid that maps to the shape of the image, significantly
|
||||
simplifying the implementation. This approach:
|
||||
|
||||
* Maps threads directly to pixel positions
|
||||
|
||||
* Simplifies coordinate calculations
|
||||
|
||||
* Enables intuitive spatial reasoning
|
||||
|
||||
* Aligns with the natural structure of images
|
||||
|
||||
Grid configuration
|
||||
------------------
|
||||
|
||||
Rather than using a single integer to represent the size of the grid, the tutorial uses a
|
||||
``dim3`` object containing three values to represent the number of block-based
|
||||
work items per dimension:
|
||||
|
||||
* **x dimension**: Width of the image
|
||||
|
||||
* **y dimension**: Height of the image
|
||||
|
||||
* **z dimension**: Set to 1 for 2D problems
|
||||
|
||||
Complete implementation
|
||||
=======================
|
||||
|
||||
Header and setup
|
||||
----------------
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
#include <hip/hip_runtime.h>
|
||||
#include <vector>
|
||||
#include "image.h"
|
||||
|
||||
The convolution kernel
|
||||
----------------------
|
||||
|
||||
Here's the complete 2D convolution kernel for image smoothing:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
__global__ void conv2d(uint8_t *image, float *mask, int image_width,
|
||||
int image_height, int mask_width, int mask_height)
|
||||
{
|
||||
int x = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
int y = blockIdx.y * blockDim.y + threadIdx.y;
|
||||
|
||||
if (x >= image_width || y >= image_height) {
|
||||
return;
|
||||
}
|
||||
|
||||
float sum = 0;
|
||||
for (int i = 0; i < mask_width; i++) {
|
||||
for (int j = 0; j < mask_height; j++) {
|
||||
// Calculate the coordinate of the pixel to read.
|
||||
int image_x = x + i - mask_width / 2;
|
||||
int image_y = y + j - mask_height / 2;
|
||||
|
||||
// Do not read outside the image.
|
||||
if (image_x < 0 || image_x >= image_width ||
|
||||
image_y < 0 || image_y >= image_height) {
|
||||
continue;
|
||||
}
|
||||
|
||||
// Accumulate the value of the pixel.
|
||||
int image_index = image_y * image_width + image_x;
|
||||
int mask_index = j * mask_width + i;
|
||||
sum += image[image_index] / 255.0f * mask[mask_index];
|
||||
}
|
||||
}
|
||||
|
||||
int image_index = y * image_width + x;
|
||||
image[image_index] = sum * 255;
|
||||
}
|
||||
|
||||
Thread identification in 2D
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
To obtain thread IDs in both x and y dimensions:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
int x = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
int y = blockIdx.y * blockDim.y + threadIdx.y;
|
||||
|
||||
This calculation combines:
|
||||
|
||||
* ``threadIdx.x`` and ``threadIdx.y``: Local thread coordinates within a block
|
||||
|
||||
* ``blockIdx.x`` and ``blockIdx.y``: Block coordinates within the grid
|
||||
|
||||
* ``blockDim.x`` and ``blockDim.y``: Block dimensions
|
||||
|
||||
Boundary checking
|
||||
~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
if (x >= image_width || y >= image_height) {
|
||||
return;
|
||||
}
|
||||
|
||||
This ensures threads don't process pixels outside the image bounds. Threads that
|
||||
exceed the image dimensions simply return without doing work.
|
||||
|
||||
Mask application loop
|
||||
~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
The nested loops iterate over the mask dimensions:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
for (int i = 0; i < mask_width; i++) {
|
||||
for (int j = 0; j < mask_height; j++) {
|
||||
// Process each mask element
|
||||
}
|
||||
}
|
||||
|
||||
Coordinate calculation
|
||||
~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
For each position in the mask, calculate the corresponding image coordinate:
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
int image_x = x + i - mask_width / 2;
|
||||
int image_y = y + j - mask_height / 2;
|
||||
|
||||
The ``- mask_width / 2`` and ``- mask_height / 2`` center the mask over the
|
||||
current pixel.
|
||||
|
||||
Edge handling
|
||||
~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
if (image_x < 0 || image_x >= image_width ||
|
||||
image_y < 0 || image_y >= image_height) {
|
||||
continue;
|
||||
}
|
||||
|
||||
This prevents reading outside the image boundaries. When the mask extends beyond
|
||||
the image edge, the code simply skips those pixels (continue to the next
|
||||
iteration).
|
||||
|
||||
Accumulation
|
||||
~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
int image_index = image_y * image_width + image_x;
|
||||
int mask_index = j * mask_width + i;
|
||||
sum += image[image_index] / 255.0f * mask[mask_index];
|
||||
|
||||
For each mask position:
|
||||
|
||||
1. Calculate the flattened array index for the image pixel
|
||||
2. Calculate the flattened array index for the mask value
|
||||
3. Normalize the pixel value (divide by 255 to get 0-1 range)
|
||||
4. Multiply by the mask weight and accumulate
|
||||
|
||||
Writing the result
|
||||
~~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
int image_index = y * image_width + x;
|
||||
image[image_index] = sum * 255;
|
||||
|
||||
After processing all mask positions, write the accumulated result back to the
|
||||
output image, scaling back to the 0-255 range.
|
||||
|
||||
Host code implementation
|
||||
------------------------
|
||||
|
||||
Main function setup
|
||||
~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
int main() {
|
||||
int width, height, channels;
|
||||
static const int maskWidth = 200;
|
||||
static const int maskHeight = 200;
|
||||
std::vector<float> mask(maskWidth * maskHeight * channels);
|
||||
|
||||
// Initialize mask with uniform averaging weights
|
||||
for (int i = 0; i < maskWidth * maskHeight; ++i) {
|
||||
mask[i] = 1.0f / maskWidth / maskHeight / channels;
|
||||
}
|
||||
|
||||
// Load an image from disk (implementation not shown)
|
||||
|
||||
Mask initialization
|
||||
~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
The mask is initialized with uniform weights that sum to 1.0:
|
||||
|
||||
* Each element is ``1.0 / (maskWidth * maskHeight * channels)``
|
||||
|
||||
* This creates an averaging filter
|
||||
|
||||
* When applied, it produces a smoothing (blurring) effect
|
||||
|
||||
Memory allocation and data transfer
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
// Allocate GPU memory and copy data to the GPU.
|
||||
uint8_t *d_image;
|
||||
float *d_mask;
|
||||
hipMalloc(&d_image, width * height * channels * sizeof(uint8_t));
|
||||
hipMalloc(&d_mask, maskWidth * maskHeight * channels * sizeof(float));
|
||||
|
||||
hipMemcpy(d_image, image, width * height * channels * sizeof(uint8_t),
|
||||
hipMemcpyHostToDevice);
|
||||
hipMemcpy(d_mask, mask.data(),
|
||||
maskWidth * maskHeight * channels * sizeof(float),
|
||||
hipMemcpyHostToDevice);
|
||||
|
||||
Grid configuration and kernel launch
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
// Calculate grid size and launch the kernel.
|
||||
dim3 block_size = {16, 16, 1};
|
||||
dim3 grid_size = {(width + block_size.x - 1) / block_size.x,
|
||||
(height + block_size.y - 1) / block_size.y, 1};
|
||||
|
||||
conv2d<<<grid_size, block_size>>>(d_image, d_mask, width, height,
|
||||
maskWidth, maskHeight);
|
||||
hipDeviceSynchronize();
|
||||
|
||||
**Grid size calculation:**
|
||||
|
||||
* :code:`block_size`: 16 × 16 threads per block (256 threads total).
|
||||
|
||||
* :code:`grid_size`: Calculated to cover the entire image.
|
||||
|
||||
* The ``(width + block_size.x - 1) / block_size.x`` formula ensures there are
|
||||
enough blocks to cover all pixels, rounding up.
|
||||
|
||||
Retrieving results and cleanup
|
||||
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
.. code-block:: c++
|
||||
|
||||
// Copy the data back to the host.
|
||||
hipMemcpy(image, d_image, width * height * channels * sizeof(uint8_t),
|
||||
hipMemcpyDeviceToHost);
|
||||
|
||||
// Store the image to disk (implementation not shown)
|
||||
|
||||
hipFree(d_image);
|
||||
hipFree(d_mask);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
Performance considerations
|
||||
==========================
|
||||
|
||||
For high-resolution images for example 4096×4096 pixels, the GPU acceleration
|
||||
provides:
|
||||
|
||||
* **Massive parallelism:** Tens of thousands of concurrent threads.
|
||||
|
||||
* **High throughput:** Leveraging GPU memory bandwidth and computational density.
|
||||
|
||||
* **Scalability:** Linear speedup with increased SM occupancy.
|
||||
|
||||
Typical speedups over CPU implementations range from **10× to 100×**, depending
|
||||
on image size, mask complexity, and GPU architecture.
|
||||
|
||||
Memory access patterns
|
||||
-----------------------
|
||||
|
||||
Image convolution requires repeated access to neighboring pixels, leading to
|
||||
non-coalesced memory transactions. Optimizing memory access is essential:
|
||||
|
||||
* Non-coalesced memory accesses (not all accesses are contiguous).
|
||||
|
||||
* Repeated reads of the same pixels by adjacent threads.
|
||||
|
||||
* Potential for optimization using shared memory (advanced technique).
|
||||
|
||||
Best practices
|
||||
==============
|
||||
|
||||
1. **Handle boundaries carefully**
|
||||
|
||||
Conditional branches in edge regions can cause wavefronts or warps to execute
|
||||
serially. Prefer approaches that avoid per-thread branching, including:
|
||||
|
||||
* Pre-clamping coordinates to valid index ranges
|
||||
|
||||
* Padding or haloing input images so all threads operate on valid data
|
||||
|
||||
* Using hardware boundary modes (such as texture sampling modes on RDNA) to
|
||||
offload boundary handling
|
||||
|
||||
These techniques help maintain high SIMD lane utilization across the entire
|
||||
grid.
|
||||
|
||||
2. **Center your stencil properly**
|
||||
|
||||
Compute the stencil origin using mask_width / 2 (or the equivalent for
|
||||
rectangular masks) to ensure correct alignment between the input data and the
|
||||
mask coefficients. This prevents off-by-one misalignment that can propagate
|
||||
as spatial artifacts.
|
||||
|
||||
3. **Select mask sizes based on compute–memory tradeoffs**
|
||||
|
||||
Larger kernels increase arithmetic intensity but also expand the set of
|
||||
neighbor loads per output element. Balance mask dimensions with available
|
||||
bandwidth, register pressure, and shared-memory capacity, particularly when
|
||||
implementing separable or multi-pass stencils.
|
||||
|
||||
4. **Normalize properly**
|
||||
|
||||
Ensure mask weights sum to the intended normalization constant, commonly 1.0
|
||||
for averaging operations. When using integer or half-precision paths, verify
|
||||
scaling behavior to avoid overflow or unintended bias.
|
||||
|
||||
5. **Consider edge strategies**
|
||||
|
||||
Adopt a clear policy for pixels whose neighborhoods extend outside the valid
|
||||
domain. Options include skipping output generation, clamping to the nearest
|
||||
valid coordinate, wrapping coordinates, or mirroring.
|
||||
|
||||
Conclusion
|
||||
==========
|
||||
|
||||
Stencil operations are a fundamental pattern in GPU computing that enables
|
||||
efficient parallel processing of spatially-dependent data. The 2D convolution
|
||||
example demonstrates:
|
||||
|
||||
* How to structure kernels for stencil patterns
|
||||
|
||||
* Proper boundary handling for neighborhood operations
|
||||
|
||||
* Effective use of 2D thread grids that map naturally to image structure
|
||||
|
||||
* Memory access patterns for adjacent data elements
|
||||
|
||||
By understanding stencil operations, developers can implement a wide range of
|
||||
image processing algorithms, scientific simulations, and deep learning
|
||||
operations on GPUs. The patterns demonstrated here extend beyond image
|
||||
processing to any computational problem involving spatial relationships in
|
||||
multi-dimensional data.
|
||||
|
||||
The key to successful stencil implementations is carefully managing boundary
|
||||
conditions, ensuring correct coordinate calculations, and leveraging the GPU's
|
||||
parallel architecture to process many independent stencil operations
|
||||
simultaneously.
|
||||
Ссылка в новой задаче
Block a user