diff --git a/projects/hip/docs/index.md b/projects/hip/docs/index.md
index 7fc82eff26..8d5aa045af 100644
--- a/projects/hip/docs/index.md
+++ b/projects/hip/docs/index.md
@@ -63,6 +63,7 @@ The HIP documentation is organized into the following categories:
 * [HIP basic examples](https://github.com/ROCm/rocm-examples/tree/develop/HIP-Basic)
 * [HIP examples](https://github.com/ROCm/rocm-examples)
 * [SAXPY tutorial](./tutorial/saxpy)
+* [GPU programming patterns](./tutorial/programming-patterns)
 * [Reduction tutorial](./tutorial/reduction)
 * [Cooperative groups tutorial](./tutorial/cooperative_groups_tutorial)
 * [HIP Graph API tutorial](./tutorial/graph_api)
diff --git a/projects/hip/docs/sphinx/_toc.yml.in b/projects/hip/docs/sphinx/_toc.yml.in
index 3b76d5454c..cd66584f0a 100644
--- a/projects/hip/docs/sphinx/_toc.yml.in
+++ b/projects/hip/docs/sphinx/_toc.yml.in
@@ -124,6 +124,14 @@ subtrees:
   - url: https://github.com/ROCm/rocm-examples
     title: HIP examples
   - file: tutorial/saxpy
+  - file: tutorial/programming-patterns
+    subtrees:
+    - entries:
+      - file: tutorial/programming-patterns/matrix_multiplication
+      - file: tutorial/programming-patterns/atomic_operations_histogram
+      - file: tutorial/programming-patterns/cpu_gpu_kmeans
+      - file: tutorial/programming-patterns/stencil_operations
+      - file: tutorial/programming-patterns/multikernel_bfs
   - file: tutorial/reduction
   - file: tutorial/cooperative_groups_tutorial
   - file: tutorial/graph_api
diff --git a/projects/hip/docs/tutorial/prerequisites.rst b/projects/hip/docs/tutorial/prerequisites.rst
new file mode 100644
index 0000000000..e8285bab03
--- /dev/null
+++ b/projects/hip/docs/tutorial/prerequisites.rst
@@ -0,0 +1,8 @@
+Prerequisites
+=============
+
+To follow this tutorial, you'll need installed drivers and a HIP compiler
+toolchain to compile your code. HIP supports compiling and running on Linux and
+Windows with AMD GPUs, the combination of install instructions is more than
+worth covering as part of this tutorial. For more information about installing
+HIP development packages, see :doc:`/install/install`.
\ No newline at end of file
diff --git a/projects/hip/docs/tutorial/programming-patterns.rst b/projects/hip/docs/tutorial/programming-patterns.rst
new file mode 100644
index 0000000000..add0ae10a6
--- /dev/null
+++ b/projects/hip/docs/tutorial/programming-patterns.rst
@@ -0,0 +1,85 @@
+.. meta::
+  :description: GPU programming patterns and tutorials
+  :keywords: AMD, ROCm, HIP, GPU, programming patterns, parallel computing, tutorial
+
+.. _gpu_programming-patterns:
+
+********************************************************************************
+GPU programming patterns
+********************************************************************************
+
+GPU programming patterns are fundamental algorithmic structures that enable
+efficient parallel computation on GPUs. Understanding these
+patterns is essential for developers looking to effectively harness the massive parallel
+processing capabilities of modern GPUs for scientific computing, machine learning,
+image processing, and other computationally intensive applications.
+
+These tutorials describe core programming patterns demonstrating how to
+efficiently implement common parallel algorithms using the HIP runtime API and
+kernel extensions. Each pattern addresses a specific computational challenge and
+provides practical implementations with detailed explanations.
+
+Common GPU programming challenges
+==================================
+
+GPU programming introduces unique challenges not present in traditional CPU
+programming:
+
+* **Memory coherence**: GPUs lack robust cache coherence mechanisms, requiring
+  careful coordination when multiple threads access shared memory.
+
+* **Race conditions**: Concurrent memory access requires atomic operations or
+  careful algorithm design.
+
+* **Irregular parallelism**: Real-world algorithms often have varying amounts of
+  parallel work across iterations.
+
+* **CPU-GPU communication**: Data transfer overhead between host and device must
+  be minimized.
+
+Tutorial overview
+=================
+
+This collection provides comprehensive tutorials on essential GPU programming
+patterns:
+
+* :doc:`Two-dimensional kernels <./programming-patterns/matrix_multiplication>`:
+  Processing grid-structured data such as matrices and images.
+
+* :doc:`Stencil operations  <./programming-patterns/stencil_operations>`:
+  Updating array elements based on neighboring values.
+
+* :doc:`Atomic operations <./programming-patterns/atomic_operations_histogram>`:
+  Ensuring data integrity during concurrent memory access.
+
+* :doc:`Multi-kernel applications <./programming-patterns/multikernel_bfs>`:
+  Coordinating multiple GPU kernels to solve complex problems.
+
+* :doc:`CPU-GPU cooperation <./programming-patterns/cpu_gpu_kmeans>`: Strategic
+  work distribution between CPU and GPU.
+
+Prerequisites
+-------------
+
+To get the most from these tutorials, you should have:
+
+* Basic understanding of C/C++ programming.
+
+* Familiarity with parallel programming concepts.
+
+* HIP runtime environment installed (see :doc:`../install/install`).
+
+* Basic knowledge of GPU architecture (recommended).
+
+Getting started
+---------------
+
+Each tutorial is self-contained and can be studied independently, though we
+recommend following the order presented for a comprehensive understanding:
+
+1. **Start with Two-dimensional kernels** to understand basic GPU thread
+   organization and memory access patterns.
+2. **Progress to stencil operations** to learn about neighborhood dependencies.
+3. **Study atomic operations** to understand concurrent memory access.
+4. **Explore multi-kernel programming** for complex algorithmic patterns.
+5. **Check CPU-GPU cooperation** to handle mixed-parallelism workloads.
diff --git a/projects/hip/docs/tutorial/programming-patterns/atomic_operations_histogram.rst b/projects/hip/docs/tutorial/programming-patterns/atomic_operations_histogram.rst
new file mode 100644
index 0000000000..1e3c58bdde
--- /dev/null
+++ b/projects/hip/docs/tutorial/programming-patterns/atomic_operations_histogram.rst
@@ -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.
diff --git a/projects/hip/docs/tutorial/programming-patterns/cpu_gpu_kmeans.rst b/projects/hip/docs/tutorial/programming-patterns/cpu_gpu_kmeans.rst
new file mode 100644
index 0000000000..210d51c6f3
--- /dev/null
+++ b/projects/hip/docs/tutorial/programming-patterns/cpu_gpu_kmeans.rst
@@ -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.
diff --git a/projects/hip/docs/tutorial/programming-patterns/matrix_multiplication.rst b/projects/hip/docs/tutorial/programming-patterns/matrix_multiplication.rst
new file mode 100644
index 0000000000..203cecfc33
--- /dev/null
+++ b/projects/hip/docs/tutorial/programming-patterns/matrix_multiplication.rst
@@ -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.
diff --git a/projects/hip/docs/tutorial/programming-patterns/multikernel_bfs.rst b/projects/hip/docs/tutorial/programming-patterns/multikernel_bfs.rst
new file mode 100644
index 0000000000..47f5a1e19b
--- /dev/null
+++ b/projects/hip/docs/tutorial/programming-patterns/multikernel_bfs.rst
@@ -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.
diff --git a/projects/hip/docs/tutorial/programming-patterns/stencil_operations.rst b/projects/hip/docs/tutorial/programming-patterns/stencil_operations.rst
new file mode 100644
index 0000000000..195543c7d2
--- /dev/null
+++ b/projects/hip/docs/tutorial/programming-patterns/stencil_operations.rst
@@ -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.