Evolution of GPU Programming Paradigms: A 15-Year Retrospective

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, establishing the GPU as the primary vehicle for massively parallel computation.

This review paper examines the evolution of GPU programming paradigms through the specific lens of multi-level parallelism—the simultaneous exploitation of task-level, work-group-level, and work-item-level parallelism that we first systematically analyzed in our 2013 M.Tech thesis [2]. That original research introduced the K-G-I (Kernels-Groups-Items) classification model and proposed API extensions for OpenCL to enable hybrid parallel execution on NVIDIA's then-new Kepler architecture with its revolutionary Hyper-Q technology.

Revisiting these concepts thirteen years later reveals a fascinating trajectory: many ideas that required custom API extensions in 2013 are now mainstream features. The K-G-I model maps directly to SYCL's hierarchical execution model; concurrent kernel execution is formalized in CUDA Graphs; and unified memory architectures have eliminated entire categories of explicit data transfer overhead.

1.1 The 2013 Problem Statement

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, with limited support for task-level parallelism across independent kernels
• False dependencies: Pre-Kepler NVIDIA architectures used a single hardware queue, causing independent kernels submitted to different software streams to serialize unnecessarily
• Manual optimization burden: Exploiting concurrent kernel execution required intricate manual scheduling, explicit stream management, and deep understanding of hardware queue behavior
• Portability vs. performance trade-off: OpenCL promised "write once, run anywhere" but consistently underperformed CUDA by 5–15% on NVIDIA hardware due to abstraction overhead

1.2 The K-G-I Classification Model

To address these challenges, we proposed the K-G-I classification model for categorizing parallel execution patterns in GPU applications:

This classification enabled systematic analysis of how algorithms could exploit multiple parallelism levels simultaneously—what we termed hybrid parallelism. Our proposed clEnqueueNDRangeHyperKernel API extension automated the distribution of independent kernels across Hyper-Q's 32 hardware queues on Kepler GPUs.

2. Historical Evolution of GPU Programming

2.1 Phase 1: Shader-Based GPGPU (2001–2006)

Before dedicated compute APIs existed, pioneering researchers exploited GPUs by repurposing graphics shaders for general computation. This approach required encoding computational problems as rendering operations—input data packed into texture formats, computations expressed as fragment shader operations.

Notable early work included Stanford BrookGPU (2004), which provided a stream programming abstraction compiled to shader code, presaging CUDA's programming model [4].

2.2 Phase 2: The CUDA Revolution (2007–2011)

NVIDIA's release of CUDA 1.0 in February 2007 [1], alongside the GeForce 8800 GTX (G80 architecture), fundamentally transformed GPU programming. Key innovations included:

• C-like kernel syntax: Kernels written in extended C/C++, dramatically lowering the barrier to entry
• Hierarchical thread organization: Threads → Thread blocks → Grids
• Memory hierarchy: Registers, shared memory, global memory, constant/texture memory
• SIMT execution model: Single Instruction, Multiple Thread with warp-based divergence handling

The Fermi architecture (2010) brought GPU computing to maturity with true cache hierarchy, ECC memory support, and concurrent kernel execution (up to 16 kernels) [10].

2.3 Phase 3: OpenCL and the Standards War (2008–2013)

Apple proposed OpenCL to the Khronos Group in June 2008, with the OpenCL 1.0 specification ratified in December 2008. OpenCL aimed to provide a vendor-neutral alternative to CUDA with platform-agnostic abstractions [5].

Despite its cross-platform promise, OpenCL faced persistent challenges: 5–15% performance gap versus CUDA, vendor fragmentation, inconsistent OpenCL 2.0 adoption, and inferior tooling ecosystem.

2.4 Phase 4: Hyper-Q and Multi-Level Parallelism (2012–2017)

The Kepler architecture (2012) introduced Hyper-Q [9], which directly addressed the concurrent kernel execution limitations that motivated our 2013 research. Hyper-Q increased hardware work queues from 1 to 32, enabling true concurrent execution of independent kernels from different streams.

Our proposed clEnqueueNDRangeHyperKernel API extension automated exploitation of Hyper-Q for OpenCL applications, achieving 1.4× speedup on Strassen's Matrix Multiplication at recursion depth 1 [2].

2.5 Phase 5: SYCL, Metal, and Ecosystem Fragmentation (2017–Present)

The period from 2017 to present has seen significant fragmentation. Apple deprecated OpenCL in macOS 10.14 (2018), recommending developers transition to Metal. SYCL emerged from the Khronos Group as a higher-level C++ abstraction, with SYCL 2020 bringing C++17 support, Unified Shared Memory, and backend independence [6].

Today's ecosystem includes:

• CUDA 12.x: Dominant for NVIDIA hardware, with CUDA Graphs for dependency-based scheduling
• SYCL 2020: Primary portable standard with implementations from Intel (DPC++), Codeplay (AdaptiveCpp)
• Metal 4: Apple's proprietary API, optimized for unified memory on M-series chips
• HIP/ROCm: AMD's CUDA-compatible layer for Radeon GPUs
• OpenCL 3.0: Maintenance mode, primarily for legacy support

3. Architectural Analysis: Execution and Memory Models

3.1 Execution Model Comparison

The three frameworks differ significantly in their programming approaches:

3.2 Kernel Launch Patterns

The three frameworks differ significantly in how kernels are defined and launched:

OpenCL NDRange Kernel Launch

// OpenCL: Separate kernel string compiled at runtime
const char* kernelSource = R"(
    __kernel void vector_add(__global float* A,
                             __global float* B,
                             __global float* C,
                             int N) {
        int gid = get_global_id(0);
        if (gid < N) {
            C[gid] = A[gid] + B[gid];
        }
    })";

// Host code
cl_program program = clCreateProgramWithSource(context, 1, &kernelSource, NULL, &err);
clBuildProgram(program, 1, &device_id, NULL, NULL, NULL);
cl_kernel kernel = clCreateKernel(program, "vector_add", &err);

clSetKernelArg(kernel, 0, sizeof(cl_mem), &bufA);
clSetKernelArg(kernel, 1, sizeof(cl_mem), &bufB);
clSetKernelArg(kernel, 2, sizeof(cl_mem), &bufC);
clSetKernelArg(kernel, 3, sizeof(int), &N);

size_t globalSize = N;
size_t localSize = 256;
clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &globalSize, &localSize, 0, NULL, NULL);

SYCL Kernel Launch with Accessors

// SYCL: Single-source C++, buffers with automatic data movement
sycl::queue q;
sycl::buffer<float> bufA(A.data(), sycl::range<1>(N));
sycl::buffer<float> bufB(B.data(), sycl::range<1>(N));
sycl::buffer<float> bufC(C.data(), sycl::range<1>(N));

q.submit([&](sycl::handler& h) {
    auto accA = bufA.get_access<sycl::access::mode::read>(h);
    auto accB = bufB.get_access<sycl::access::mode::read>(h);
    auto accC = bufC.get_access<sycl::access::mode::write>(h);

    h.parallel_for<class VectorAdd>(sycl::range<1>(N), [=](sycl::id<1> i) {
        accC[i] = accA[i] + accB[i];
    });
}).wait();

SYCL with Unified Shared Memory (USM)

// SYCL USM: Pointer-based API, explicit memory management
sycl::queue q;
float* A = sycl::malloc_shared<float>(N, q);
float* B = sycl::malloc_shared<float>(N, q);
float* C = sycl::malloc_shared<float>(N, q);

// Initialize A and B...

q.parallel_for<class VectorAddUSM>(sycl::range<1>(N), [=](sycl::id<1> i) {
    C[i] = A[i] + B[i];
}).wait();

sycl::free(A, q);
sycl::free(B, q);
sycl::free(C, q);

Metal Shader File (MatrixMultiply.metal)

#include <metal_stdlib>
using namespace metal;

kernel void vector_add(device const float* A [[buffer(0)]],
                       device const float* B [[buffer(1)]],
                       device float* C [[buffer(2)]],
                       uint gid [[thread_position_in_grid]]) {
    C[gid] = A[gid] + B[gid];
}

Metal Host Code (Swift)

// Swift host code for Metal
let device = MTLCreateSystemDefaultDevice()!
let library = device.makeDefaultLibrary()!
let function = library.makeFunction(name: "vector_add")!
let pipeline = try! device.makeComputePipelineState(function: function)

let commandQueue = device.makeCommandQueue()!
let commandBuffer = commandQueue.makeCommandBuffer()!
let encoder = commandBuffer.makeComputeCommandEncoder()!

encoder.setComputePipelineState(pipeline)
encoder.setBuffer(bufferA, offset: 0, index: 0)
encoder.setBuffer(bufferB, offset: 0, index: 1)
encoder.setBuffer(bufferC, offset: 0, index: 2)

let gridSize = MTLSize(width: N, height: 1, depth: 1)
let threadgroupSize = MTLSize(width: min(256, pipeline.maxTotalThreadsPerThreadgroup), height: 1, depth: 1)
encoder.dispatchThreads(gridSize, threadsPerThreadgroup: threadgroupSize)

encoder.endEncoding()
commandBuffer.commit()
commandBuffer.waitUntilCompleted()

3.3 Memory Model Deep Dive

Apple Silicon's unified memory architecture represents a fundamentally different approach that eliminates explicit data transfers:

3.4 Synchronization Mechanisms

The K-G-I model identifies three synchronization levels, each with distinct mechanisms across frameworks.

K-Level Synchronization (Between Kernels)

Kernel-level synchronization ensures ordering between independent kernel invocations:

// OpenCL: Event-based dependencies
cl_event event1, event2;
clEnqueueNDRangeKernel(queue, kernel1, ..., 0, NULL, &event1);
clEnqueueNDRangeKernel(queue, kernel2, ..., 1, &event1, &event2);
// kernel2 waits for kernel1 to complete

// SYCL: Implicit dependencies via accessors
q.submit([&](handler& h) {
    auto acc = buf.get_access<access::mode::write>(h);
    h.parallel_for<class K1>(range<1>(N), [=](id<1> i) { acc[i] = i; });
});

q.submit([&](handler& h) {
    auto acc = buf.get_access<access::mode::read>(h);
    // Implicit dependency: waits for K1 to complete
    h.parallel_for<class K2>(range<1>(N), [=](id<1> i) { result[i] = acc[i] * 2; });
});

// Metal: Command buffer ordering with fences
let commandBuffer1 = commandQueue.makeCommandBuffer()!
// ... encode kernel1 ...
commandBuffer1.commit()

let commandBuffer2 = commandQueue.makeCommandBuffer()!
// ... encode kernel2 ...
commandBuffer2.commit()
// Sequential by default; use MTLSharedEvent for explicit sync

G-Level Synchronization (Within Work-Group)

I-Level Synchronization (Between Work-Items)

Atomic operations for fine-grained synchronization:

// OpenCL atomics
atomic_add(&shared_counter, 1);

// SYCL atomics (C++ style)
sycl::atomic_ref<int, memory_order::relaxed, memory_scope::device,
                 access::address_space::global_space> ref(counter[0]);
ref.fetch_add(1);

// Metal atomics
atomic_fetch_add_explicit(&counter, 1, memory_order_relaxed);

4. Systematic Framework Comparison

4.1 Platform Support

4.2 Advanced Features

4.3 Unified Memory Comparison: USM vs SVM

4.4 ML Framework Backend Support

4.5 Development Ecosystem

5. Multi-Level Parallelism: K-G-I Model Application

This section applies our K-G-I classification model to analyze multi-level parallelism support across modern frameworks.

5.1 K-G-I Mapping to Modern Frameworks

5.2 Evolution from Hyper-Q to CUDA Graphs

The concepts we explored in 2013 with Hyper-Q have evolved significantly. Modern CUDA Graphs [8] provide:

• Explicit dependency specification: DAG structure is clear and verifiable
• Reduced launch overhead: Graph instantiation amortizes setup cost
• Optimization opportunities: Runtime can optimize node placement
• Repeatability: Same graph can be launched repeatedly without reconstruction

2013 Approach: Manual Multi-Stream Scheduling

// Manual stream management for Strassen (2013)
cudaStream_t streams[7];
for (int i = 0; i < 7; i++) {
    cudaStreamCreate(&streams[i]);
}

// Launch P1-P7 on separate streams
strassen_multiply<<<grid, block, 0, streams[0]>>>(A1, B1, P1);
strassen_multiply<<<grid, block, 0, streams[1]>>>(A2, B2, P2);
strassen_multiply<<<grid, block, 0, streams[2]>>>(A3, B3, P3);
strassen_multiply<<<grid, block, 0, streams[3]>>>(A4, B4, P4);
strassen_multiply<<<grid, block, 0, streams[4]>>>(A5, B5, P5);
strassen_multiply<<<grid, block, 0, streams[5]>>>(A6, B6, P6);
strassen_multiply<<<grid, block, 0, streams[6]>>>(A7, B7, P7);

// Synchronize all streams before recomposition
for (int i = 0; i < 7; i++) {
    cudaStreamSynchronize(streams[i]);
}

2024 Approach: CUDA Graphs with Explicit Dependencies

// CUDA Graphs for Strassen (modern approach)
cudaGraph_t graph;
cudaGraphCreate(&graph, 0);

// Add kernel nodes with dependencies
cudaGraphNode_t nodes[11]; // 7 multiplies + 4 recompose

// P1-P7 are independent (no dependencies)
for (int i = 0; i < 7; i++) {
    cudaKernelNodeParams params = {...};
    cudaGraphAddKernelNode(&nodes[i], graph, NULL, 0, &params);
}

// Recomposition depends on all P nodes
cudaGraphNode_t deps[] = {nodes[0], nodes[1], nodes[2], nodes[3],
                          nodes[4], nodes[5], nodes[6]};
cudaGraphAddKernelNode(&nodes[7], graph, deps, 7, &recomposeC11);
cudaGraphAddKernelNode(&nodes[8], graph, deps, 7, &recomposeC12);
cudaGraphAddKernelNode(&nodes[9], graph, deps, 7, &recomposeC21);
cudaGraphAddKernelNode(&nodes[10], graph, deps, 7, &recomposeC22);

// Instantiate and launch
cudaGraphExec_t instance;
cudaGraphInstantiate(&instance, graph, NULL, NULL, 0);
cudaGraphLaunch(instance, stream);

5.3 SYCL's Hierarchical Parallelism

SYCL 2020 provides native support for multi-level parallelism through its hierarchical execution model with parallel_for_work_group and parallel_for_work_item constructs that map directly to our K-G-I model.

// SYCL Hierarchical Parallelism Example
q.submit([&](handler& cgh) {
    auto acc = buf.get_access<access::mode::read_write>(cgh);

    cgh.parallel_for_work_group<class Hierarchical>(
        range<1>(numGroups),    // Number of work-groups
        range<1>(groupSize),    // Work-group size
        [=](group<1> g) {
            // G-level: Code here runs once per work-group
            int groupId = g.get_id(0);

            // Shared memory allocation at G-level
            int localSum = 0;

            g.parallel_for_work_item([&](h_item<1> item) {
                // I-level: Code here runs for each work-item
                int localId = item.get_local_id(0);
                int globalId = groupId * groupSize + localId;

                // Each work-item contributes to local sum
                localSum += acc[globalId];
            });

            // Implicit barrier after parallel_for_work_item

            // G-level: Write group result
            if (groupId < numGroups) {
                groupResults[groupId] = localSum;
            }
        });
});

Metal Multi-Level Parallelism

// Metal hierarchical parallelism with SIMD groups
kernel void hierarchical_reduction(
    device float* input [[buffer(0)]],
    device float* output [[buffer(1)]],
    threadgroup float* shared [[threadgroup(0)]],
    uint tid [[thread_index_in_threadgroup]],
    uint gid [[threadgroup_position_in_grid]],
    uint simd_lane [[thread_index_in_simdgroup]],
    uint simd_group [[simdgroup_index_in_threadgroup]])
{
    // I-level: Each thread loads one element
    shared[tid] = input[gid * 256 + tid];

    threadgroup_barrier(mem_flags::mem_threadgroup);

    // SIMD-level reduction (within 32-thread SIMD group)
    float val = shared[tid];
    val += simd_shuffle_down(val, 16);
    val += simd_shuffle_down(val, 8);
    val += simd_shuffle_down(val, 4);
    val += simd_shuffle_down(val, 2);
    val += simd_shuffle_down(val, 1);

    // First lane of each SIMD group writes to shared memory
    if (simd_lane == 0) {
        shared[simd_group] = val;
    }

    threadgroup_barrier(mem_flags::mem_threadgroup);

    // G-level: First SIMD group reduces across groups
    if (simd_group == 0 && simd_lane < 8) {
        val = shared[simd_lane];
        val += simd_shuffle_down(val, 4);
        val += simd_shuffle_down(val, 2);
        val += simd_shuffle_down(val, 1);

        if (simd_lane == 0) {
            output[gid] = val;
        }
    }
}

5.4 Case Study: Strassen's Algorithm Revisited

6. Performance Analysis

6.1 SYCL vs CUDA Performance

On NVIDIA hardware, SYCL (via CUDA backend) achieves near-native performance [11]:

Values normalized to native CUDA performance (1.00×)

6.2 Apple Silicon Efficiency

While raw performance favors discrete GPUs, power efficiency tells a different story:

Apple Silicon achieves 2.5–3.5× better GFLOPS/W than discrete GPUs, making it compelling for power-constrained deployments.

6.3 Multi-Level Parallelism Impact

Comparing our 2013 measurements on Kepler with modern equivalents:

7. Future Directions

7.1 Emerging Hardware Trends

Chiplet-Based GPU Architectures: AMD's MI300 and future NVIDIA designs use chiplet architectures, introducing new programming challenges:

• Non-uniform memory access between chiplets
• Inter-chiplet communication requiring explicit management
• Heterogeneous configurations mixing compute and memory chiplets

The K-G-I model may need extension to K-C-G-I (Chiplets) for work placement across NUMA-aware topologies.

7.2 Framework Evolution Predictions

• SYCL: Positioned to become the default portable GPU programming standard, with potential Metal backend for Apple Silicon
• Metal: Will remain Apple-exclusive but continue evolving with deeper MLX integration
• OpenCL: Primarily maintenance mode; gradual migration to SYCL for new development

7.3 Framework Selection Guide

8. Conclusion

The GPU programming landscape has undergone remarkable transformation over the fifteen years since CUDA's introduction. What began as a two-framework competition between CUDA and OpenCL has evolved into a fragmented but maturing ecosystem with distinct solutions for different needs.

8.1 Key Findings

1. CUDA maintains dominance through ecosystem depth, tooling quality, and continuous innovation. CUDA Graphs formalized the concurrent kernel execution concepts we explored in 2013.
2. SYCL has emerged as the credible portable alternative, achieving 90–102% of native CUDA performance through backend flexibility while providing modern C++17 programming model.
3. Metal demonstrates the power of vertical integration, with unified memory eliminating entire categories of programming complexity on Apple Silicon.
4. OpenCL is effectively deprecated for new development, surviving primarily for legacy support and embedded applications.
5. The K-G-I classification model remains relevant, though extension to A-K-G-I (Accelerators) is needed for tensor cores and neural engines.

8.2 Retrospective on 2013 Research

Revisiting our 2013 work on OpenCL API extensions reveals that:

• The clEnqueueNDRangeHyperKernel concept is now embodied in CUDA Graphs' dependency-aware execution model
• K-G-I classification maps directly to SYCL's hierarchical parallel execution
• Unified memory (then experimental SVM) is now a core feature in SYCL USM and native on Apple Silicon
• Multi-level parallelism is no longer an advanced technique requiring API extensions—it is the default execution model

References

1. NVIDIA Corporation. (2007). CUDA Programming Guide 1.0.
2. Kasiviswanathan, S. (2013). OpenCL API Extensions to achieve Multi-level Parallelism for Efficient Implementation of Strassen's Matrix Multiplication on GPUs. M.Tech Thesis, Indian Institute of Science, Bangalore.
3. Owens, J. D., et al. (2007). A Survey of General-Purpose Computation on Graphics Hardware. Computer Graphics Forum, 26(1), 80–113.
4. Buck, I., et al. (2004). Brook for GPUs: Stream Computing on Graphics Hardware. ACM SIGGRAPH.
5. Khronos Group. (2020). OpenCL Specification, Version 3.0.
6. Khronos Group. (2023). SYCL 2020 Specification (revision 9).
7. Apple Inc. (2024). Metal Programming Guide.
8. NVIDIA Corporation. (2024). CUDA Graphs. CUDA C++ Programming Guide.
9. NVIDIA Corporation. (2012). NVIDIA's Next Generation CUDA Compute Architecture: Kepler GK110. Whitepaper.
10. NVIDIA Corporation. (2009). NVIDIA's Next Generation CUDA Compute Architecture: Fermi. Whitepaper.
11. Alpay, A., Heß, B., & Heuveline, V. (2024). SYCL-Bench 2020: Benchmarking SYCL 2020 on AMD, Intel, and NVIDIA GPUs. IWOCL.
12. Apple Machine Learning Research. (2023). MLX: An Array Framework for Apple Silicon.
13. Strassen, V. (1969). Gaussian elimination is not optimal. Numerische Mathematik, 13, 354–356.