OpenCL API Extensions to achieve Multi-level Parallelism for Efficient Implementation of Strassen's Matrix Multiplication on GPUs

Abstract

1. Introduction

The transformation of graphics processing units (GPUs) from fixed-function rendering pipelines to fully programmable parallel processors represents one of the most significant shifts in computing architecture of the 21st century. When NVIDIA released CUDA 1.0 in February 2007, it initiated a revolution that would fundamentally reshape domains from scientific simulation to artificial intelligence.

However, GPU programming models have historically focused on single-level parallelism: expressing data parallelism at the work-item level with limited support for task-level parallelism across independent kernels. This limitation becomes problematic for algorithms that exhibit natural hierarchical parallelism.

1.1 The Problem

In 2013, the GPU programming landscape was characterized by:

• Single-level parallelism focus: Both CUDA and OpenCL primarily expressed data parallelism at the work-item level
• False dependencies: Pre-Kepler architectures used a single hardware queue, causing independent kernels to serialize unnecessarily
• Manual optimization burden: Exploiting concurrent kernel execution required intricate manual scheduling and deep hardware knowledge
• Portability vs. performance trade-off: OpenCL consistently underperformed CUDA by 5–15% on NVIDIA hardware

1.2 NDRange Limitations

The OpenCL NDRange execution model, while powerful for regular data parallelism, imposes significant constraints that limit exploitation of complex parallel patterns:

• Homogeneous work-items: All work-items in an NDRange execute the same kernel code, making heterogeneous computations difficult to express
• Uniform work-group sizes: All work-groups must have identical dimensions, preventing dynamic load balancing
• Single kernel granularity: One NDRange corresponds to one kernel invocation; independent computations require separate enqueues
• Limited inter-kernel parallelism: While command queues support out-of-order execution, hardware limitations (single queue pre-Kepler) caused serialization

1.3 Our Contribution

This thesis makes three main contributions:

1. K-G-I Classification Model: A systematic framework for categorizing parallel execution patterns at three levels: Kernels (K), work-Groups (G), and work-Items (I)
2. API Extension Proposal: The clEnqueueNDRangeHyperKernel function that automates multi-stream execution on Hyper-Q hardware
3. Strassen Case Study: Demonstration of multi-level parallelism benefits using Strassen's Matrix Multiplication algorithm

2. Related Work

Research on exploiting parallelism in matrix algorithms and GPU programming has a rich history. We categorize prior work into algorithmic optimizations and system-level approaches.

2.1 Algorithmic Approaches

Cohen [9] analyzed the overhead of Strassen's recursion, establishing that the crossover point where Strassen outperforms conventional multiplication depends critically on the constant factors in implementation—a finding that remains relevant for GPU implementations where kernel launch overhead must be balanced against reduced arithmetic operations.

West et al. [10] developed the Mentat Programming Language, introducing the concept of "medium-grained" parallelism through their MPL construct for expressing parallel object invocations. Their work influenced our thinking on K-level parallelism as a distinct abstraction layer.

2.2 Hierarchical and Tiled Approaches

Bikshandi et al. [11] introduced hierarchically tiled arrays (HTAs), providing a high-level abstraction for expressing tiled computations. HTAs capture G-level parallelism naturally but lack explicit support for K-level concurrent kernel execution.

Subhlok et al. [12] explored the integration of task and data parallelism in scientific codes, demonstrating speedups when both levels are exploited. Their work directly motivates our K-G-I classification.

2.3 Language and Runtime Systems

BenHassen et al. [13] developed the Orca language for expressing multi-level parallelism on distributed memory systems. Kambadur et al. [14] introduced the PFunc library for parallel task creation with work-stealing schedulers.

2.4 GPU-Specific Work

Chao-Tung et al. [15] examined multi-GPU cluster configurations for linear algebra, but focused on inter-GPU parallelism rather than intra-GPU multi-level parallelism.

Vejarano et al. [16] developed analytical cost models for GPU execution, providing foundation techniques that inform our cost model development.

Most recently, Elangovan et al. [17] explored the OmpSs task-based programming model for GPU clusters, demonstrating task parallelism at the cluster level but not exploiting intra-GPU concurrent kernel execution.

3. GPU Architecture Background

3.1 Hardware Architecture

Modern GPUs consist of an array of Streaming Multiprocessors (SMXs), each containing multiple Streaming Processors (SPs) or CUDA cores. On Kepler K20, the architecture comprises 13 SMXs, each with 192 SPs, yielding 2496 total CUDA cores.

3.2 CUDA Programming Model

NVIDIA's CUDA introduced foundational concepts that remain relevant today:

• Threads: Individual execution units, each with unique thread ID
• Thread blocks: Groups of threads that can synchronize via shared memory
• Grids: Collections of thread blocks executing the same kernel
• SIMT execution: Single Instruction, Multiple Thread model where threads in a warp (32 threads) execute the same instruction

3.3 OpenCL Correspondence

OpenCL provides equivalent abstractions with different terminology:

3.4 Kepler Architecture and Hyper-Q

The Kepler architecture (2012) introduced Hyper-Q, which directly addressed the concurrent kernel execution limitations that motivated this research. Prior to Kepler, NVIDIA GPUs used a single hardware work queue. Even with multiple CUDA streams, kernels from different streams could serialize due to false dependencies at the hardware queue level.

Hyper-Q increased the number of hardware work queues from 1 to 32, enabling:

• True concurrent execution of independent kernels from different streams
• MPI applications with multiple ranks per GPU without serialization
• Reduced latency for small kernel launches through queue parallelism

4. Strassen's Matrix Multiplication Algorithm

Strassen's algorithm [1] reduces the complexity of matrix multiplication from O(n³) to O(n^2.807) by replacing one multiplication with additions at each recursive level.

4.1 Matrix Decomposition

For n×n matrices A and B, we partition each into four (n/2)×(n/2) submatrices:

A = ┌ A₁₁  A₁₂ ┐    B = ┌ B₁₁  B₁₂ ┐    C = ┌ C₁₁  C₁₂ ┐
    └ A₂₁  A₂₂ ┘        └ B₂₁  B₂₂ ┘        └ C₂₁  C₂₂ ┘

4.2 Intermediate Matrices

We first compute auxiliary matrices from A and B:

4.3 Seven Independent Products

The core insight: compute only 7 matrix products instead of 8:

4.4 Result Recomposition

The final result C is assembled from the seven products:

4.5 Parallelism Opportunities

Strassen's algorithm exhibits natural multi-level parallelism:

• K-level (Kernels): The 7 intermediate products P₁–P₇ are completely independent at each recursion level
• G-level (Groups): Each sub-matrix multiplication can be tiled, with tiles processed by independent work-groups using shared memory
• I-level (Items): Element-wise additions and the innermost multiplication loop can be parallelized across work-items

4.6 Recursive Structure

At recursion depth d, the algorithm generates:

5. The K-G-I Classification Model

We propose a systematic framework for categorizing parallel execution patterns in GPU applications. The model captures three orthogonal dimensions of parallelism.

5.1 Parallelism Levels

5.2 Classification Modes

We classify parallel modes using notation [K]-[G]-[I] where each position indicates homogeneous (M) or heterogeneous (T) behavior:

5.3 Kernel Implementations

Each K-G-I mode requires different kernel structures. We present representative OpenCL kernel code for each mode.

Algorithm 1: S-M-T Mode Kernel (Classic Element-wise)

__kernel void smt_kernel(
    __global float* A,
    __global float* B,
    __global float* C,
    int N)
{
    int gid = get_global_id(0);

    // All work-groups execute same code
    // Each work-item operates on different data element
    if (gid < N) {
        C[gid] = A[gid] + B[gid];
    }
}

Algorithm 2: S-T-M Mode Kernel (Group-Heterogeneous)

__kernel void stm_kernel(
    __global float* data,
    __global float* result,
    int N)
{
    int group_id = get_group_id(0);
    int lid = get_local_id(0);
    __local float scratch[256];

    // Different groups perform different operations
    switch(group_id % 4) {
        case 0: // Sum reduction
            scratch[lid] = data[group_id * 256 + lid];
            break;
        case 1: // Max reduction
            scratch[lid] = data[group_id * 256 + lid];
            break;
        case 2: // Min reduction
            scratch[lid] = data[group_id * 256 + lid];
            break;
        case 3: // Product reduction
            scratch[lid] = data[group_id * 256 + lid];
            break;
    }
    barrier(CLK_LOCAL_MEM_FENCE);

    // Parallel reduction within work-group
    for(int s = 128; s > 0; s >>= 1) {
        if(lid < s) {
            // Operation depends on group_id
            scratch[lid] += scratch[lid + s]; // simplified
        }
        barrier(CLK_LOCAL_MEM_FENCE);
    }
}

Algorithm 3: S-M-M Mode Kernel (Embarrassingly Parallel)

__kernel void smm_kernel(
    __global float* output,
    unsigned int seed,
    int N)
{
    int gid = get_global_id(0);

    // All work-items execute identical code
    // (e.g., Monte Carlo simulation)
    unsigned int state = seed + gid;
    float sum = 0.0f;

    for(int i = 0; i < 1000; i++) {
        // PRNG iteration
        state = state * 1103515245 + 12345;
        float r = (float)(state & 0x7FFFFFFF) / (float)0x7FFFFFFF;
        sum += r * r;
    }

    output[gid] = sum / 1000.0f;
}

5.4 Strassen in K-G-I Terms

Strassen's algorithm maps to the K-G-I model as follows:

• Recursion Level 0 (root): Single kernel launch (S-*-*)
• Recursion Level 1: 7 independent kernels → M-M-M mode
• Recursion Level 2: 49 independent kernels → M-M-M mode
• Base case (tiled GEMM): S-M-T mode with tiled shared memory access

6. Theoretical Cost Model

We develop a cost model to predict execution time under different parallelism configurations.

6.1 Architecture Parameters

Our analysis targets the NVIDIA Kepler K20 architecture, compared against the previous-generation Fermi:

6.2 Model Parameters

6.3 Pilot Program Methodology

To calibrate our cost model, we use pilot programs that measure fundamental timing parameters:

6.4 Execution Time Analysis

For Strassen at recursion depth d with n concurrent streams:

The model predicts optimal recursion depth where:

• d = 1: 7 concurrent kernels, theoretical speedup up to 7×
• d = 2: 49 concurrent kernels, limited by H = 32 queues
• d ≥ 3: Launch overhead begins to dominate

7. Proposed API Extension

We propose an OpenCL API extension that automates the exploitation of Hyper-Q for applications with multi-kernel parallelism.

7.1 M-M-M Mode: Multi-Kernel Parallelism

The M-M-M (Multiple-hoMogeneous-hoMogeneous) mode represents true concurrent kernel execution, exploiting Kepler's Hyper-Q capability. This mode is essential for algorithms like Strassen where multiple independent kernels can execute simultaneously.

7.2 API Design

cl_int clEnqueueNDRangeHyperKernel(
    cl_command_queue queue,
    cl_uint num_kernels,
    const cl_kernel* kernels,
    const size_t** global_work_sizes,
    const size_t** local_work_sizes,
    const cl_kernel_dependency* dependencies,
    cl_event* event
);

// Dependency specification
typedef struct {
    cl_uint kernel_index;      // Dependent kernel
    cl_uint depends_on_index;  // Dependency target
    cl_dependency_type type;   // DATA, CONTROL, or NONE
} cl_kernel_dependency;

7.3 Functionality

The proposed API:

1. Accepts multiple kernels with explicit dependency specifications
2. Analyzes the dependency DAG to identify independent kernel sets
3. Automatically creates and manages CUDA streams (up to H = 32)
4. Distributes independent kernels across available hardware queues
5. Respects data and control dependencies through appropriate synchronization
6. Returns a single event representing completion of all kernels

8. Experimental Results

We evaluated our approach on a Tesla K20 (Kepler GK110) with 2496 CUDA cores and 5 GB GDDR5 memory.

8.1 Performance Results

8.2 Theoretical vs. Practical Results

Our cost model predictions compared against measured execution times:

The gap between theoretical and practical speedup is attributed to:

• Memory bandwidth saturation: Concurrent kernels compete for limited global memory bandwidth
• SMX occupancy: Not all SMXs can be fully utilized when kernels have varying resource requirements
• Launch overhead accumulation: While Hyper-Q reduces serialization, launch costs still accumulate

8.3 Analysis

Key findings:

• Speedup increases with matrix size due to better compute/launch ratio
• Depth 2 outperforms depth 1 for large matrices despite exceeding 32 queues
• Diminishing returns beyond depth 2 due to kernel launch overhead
• OpenCL overhead vs. CUDA is 8–12% for this workload

9. Conclusion

This thesis has addressed the challenge of exploiting multi-level parallelism in GPU applications. Our contributions include:

1. The K-G-I classification model provides a systematic framework for analyzing parallel execution patterns
2. The proposed clEnqueueNDRangeHyperKernel API automates multi-stream execution on Hyper-Q hardware
3. Experimental validation using Strassen's algorithm demonstrates 1.4× speedup through concurrent kernel execution

9.1 Future Work

Several directions warrant further investigation:

• Automatic parallelism extraction from sequential kernel code
• Dynamic load balancing across heterogeneous accelerators
• Integration with emerging GPU architectures beyond Kepler
• Application to other algorithms with natural multi-level parallelism

References

1. Strassen, V. (1969). Gaussian elimination is not optimal. Numerische Mathematik, 13, 354–356.
2. NVIDIA Corporation. (2012). NVIDIA's Next Generation CUDA Compute Architecture: Kepler GK110. Whitepaper.
3. Khronos Group. (2012). OpenCL Specification, Version 1.2.
4. Kirk, D. B., & Hwu, W. W. (2010). Programming Massively Parallel Processors. Morgan Kaufmann.
5. Volkov, V., & Demmel, J. W. (2008). Benchmarking GPUs to tune dense linear algebra. SC '08: Proceedings of the ACM/IEEE Conference on Supercomputing.
6. Nandy, S. K., et al. (2009). REDEFINE: Runtime Reconfigurable Polymorphic ASIC. ACM TRETS, 2(2).
7. Buck, I., et al. (2004). Brook for GPUs: Stream Computing on Graphics Hardware. ACM SIGGRAPH.
8. Owens, J. D., et al. (2007). A Survey of General-Purpose Computation on Graphics Hardware. Computer Graphics Forum, 26(1), 80–113.
9. Cohen, H. (1976). On the Overhead of Strassen's Algorithm. SIAM Journal on Computing, 5(3), 403–413.
10. West, E. A., & Grimshaw, A. S. (1995). Braid: Integrating Task and Data Parallelism. SIGPLAN Symposium on Principles and Practice of Parallel Programming.
11. Bikshandi, G., et al. (2006). Programming for Parallelism and Locality with Hierarchically Tiled Arrays. PPoPP '06.
12. Subhlok, J., & Vondran, G. (1997). Optimal Latency-Throughput Tradeoffs for Data Parallel Pipelines. SPAA '97.
13. BenHassen, S., et al. (1996). The Orca Parallel Programming Model. IEEE TPDS, 7(10).
14. Kambadur, P., et al. (2007). PFunc: Modern Task Parallelism for Modern High Performance Computing. SC '07.
15. Yang, C. T., et al. (2011). Implementation of Parallel Linear Algebra Functions on Multi-GPU Clusters. IPDPS Workshops.
16. Vejarano, G., & Gao, Y. (2011). GPU Kernel Optimization for Finite Element Computations. Parallel Computing, 37(8), 366–379.
17. Elangovan, V. K., et al. (2013). Hybrid Task and Data Parallel Programming Using OmpSs and GPUs. Euro-Par 2013.