DEV Community: Javad

Advanced GPU Optimization: How to tech an LLM with CUDA and ROCm? - Part 5 (Final Part)

Javad — Fri, 21 Aug 2026 07:02:18 +0000

Welcome back, you absolute madman! You finished Part 4, implemented Flash Attention, and squeezed FP8 out of your silicon. But the industry doesn't stop at dense Transformers. In 2024/2025, every major model (Grok, Mixtral, Gemini) uses Mixture of Experts (MoE) to scale to trillions of parameters without exploding compute costs.

Furthermore, if you actually try to run these monsters on a single node, you will hit the VRAM wall instantly. That is where CPU Offloading and Zero-Inference come to the rescue.

In this fifth and (I swear) final part, we will:

Implement MoE routing and Expert Parallelism using all-to-all communication.
Build a ZeRO-Offload mechanism to spill optimizer states to system RAM.
Harness CUDA/HIP Graphs to eliminate kernel launch overhead.
Build an optimized Inference Server with KV Caching and PagedAttention.

Prerequisites: Parts 1-4, a multi-GPU setup (4+ is ideal), and a realization that training is only half the battle—serving is the other half.

Mixture of Experts (MoE) – Routing on the GPU

A dense 1.5 Trillion parameter model would require 3TB of VRAM just for weights. MoE solves this by having hundreds of "expert" FFNs, but only activating 2 of them per token.

The Math:
y = Σ (Softmax(Router(x))_i * Expert_i(x)) for the top-k experts (usually k=2).

1.1 The Router Kernel (Top-k Gating)

First, we need a kernel that takes the token embeddings x and computes the routing scores, then selects the top-2 indices and their weights.

__global__ void moe_router_kernel(const float* x, float* gate_logits, 
                                  int* expert_indices, float* expert_weights,
                                  int num_tokens, int num_experts) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    if (tid >= num_tokens) return;

    // Compute logits for this token (dot product with router weights)
    // Assume router weights are in shared memory or loaded via L2 cache
    float max1 = -INFINITY, max2 = -INFINITY;
    int idx1 = -1, idx2 = -1;

    for (int e = 0; e < num_experts; ++e) {
        float score = 0.0f;
        for (int d = 0; d < D; ++d) {
            score += x[tid * D + d] * router_weights[e * D + d];
        }
        gate_logits[tid * num_experts + e] = score;

        // Track top-2 (manual reduction)
        if (score > max1) {
            max2 = max1; idx2 = idx1;
            max1 = score; idx1 = e;
        } else if (score > max2) {
            max2 = score; idx2 = e;
        }
    }

    // Apply Softmax only to the top-2 scores (sparse softmax)
    float denom = expf(max1 - max1) + expf(max2 - max1); // numerically stable
    expert_indices[tid * 2] = idx1;
    expert_indices[tid * 2 + 1] = idx2;
    expert_weights[tid * 2] = 1.0f / denom; 
    expert_weights[tid * 2 + 1] = expf(max2 - max1) / denom;
}

1.2 Expert Parallelism with All-to-All Communication

In dense models, we used All-Reduce. In MoE, we use All-to-All because different GPUs hold different experts. Tokens must be sent to the GPU that owns their assigned expert.

Implementation using NCCL/RCCL:

// Step 1: Token dispatch (send tokens to expert owners)
// We build a send buffer per rank based on the routing decisions.
ncclGroupStart();
for (int r = 0; r < world_size; ++r) {
    // Send tokens assigned to expert owned by rank 'r'
    if (send_counts[r] > 0) {
        COMM_SEND(send_buffer[r], send_counts[r] * D, COMM_FLOAT, r, comm, stream);
    }
    // Recv tokens for local experts
    COMM_RECV(recv_buffer[r], recv_counts[r] * D, COMM_FLOAT, r, comm, stream);
}
ncclGroupEnd();

// Step 2: Run local experts (FFN) on the received tokens locally.
run_local_experts_kernel<<<...>>>(recv_buffer, local_expert_weights, ...);

// Step 3: All-to-All again to return the processed tokens to their original GPUs.

This ensures we scale to thousands of experts across a cluster with minimal communication overhead.

ZeRO-Offload – Breaking the VRAM Wall

When you have 70B parameters, even ZeRO-3 might not fit if you only have 40GB VRAM. ZeRO-Offload moves the optimizer states (momentum, variance) and sometimes the gradients to CPU RAM (DDR).

The trick is overlapping the GPU computation with PCIe transfers. During the backward pass, we asynchronously copy gradients to CPU, compute the AdamW update on the CPU (using a separate thread), and copy the updated FP32 weights back to the GPU just in time for the next forward pass.

C++ implementation using pinned memory and streams:

// Allocate pinned memory on CPU for offloaded states
float* cpu_momentum;
float* cpu_variance;
hipHostMalloc(&cpu_momentum, num_params * sizeof(float), hipHostMallocDefault);

// During the training step:
void offloaded_train_step() {
    // 1. Forward/Backward on GPU (FP16) -> dW stays in GPU memory temporarily.
    backward_pass(...);

    // 2. Asynchronous D2H copy of gradients for the offloaded partition
    //    while GPU continues computing the next layer.
    hipMemcpyAsync(cpu_gradients, d_gradients, partition_size, 
                   hipMemcpyDeviceToHost, data_stream);

    // 3. CPU thread computes: m = beta1*m + (1-beta1)*grad, etc.
    //    (This runs on a std::async thread to not block the main loop)
    cpu_adam_update(cpu_momentum, cpu_variance, cpu_gradients, partition_size);

    // 4. Asynchronous H2D copy of updated weights back to GPU
    hipMemcpyAsync(d_model_weights, cpu_weights, partition_size, 
                   hipMemcpyHostToDevice, data_stream);

    // 5. Synchronize streams at the end.
}

Why this works: PCIe Gen 5.0 can do ~32 GB/s. If we overlap this with the 1-2 seconds it takes to compute a backward pass on a large model, the offloading overhead becomes nearly invisible.

CUDA/HIP Graphs – Killing Kernel Launch Latency

If you profile a real LLM, you will see that kernel launch overhead (the CPU time spent telling the GPU to do things) is surprisingly high—often 10-20 microseconds per kernel. A GPT-3 forward pass has ~1,000 kernels. That's 20ms wasted just on launching.

CUDA Graphs (and HIP Graphs) record a sequence of kernel launches and replay them with a single API call. This is crucial for inference.

Implementation:

hipGraph_t graph;
hipGraphExec_t instance;

// 1. Capture the graph in a stream
hipStreamBeginCapture(stream, hipStreamCaptureModeGlobal);
    // Launch all kernels for the forward pass
    layernorm_kernel<<<..., stream>>>(...);
    matmul_kernel<<<..., stream>>>(...);
    flash_attention<<<..., stream>>>(...);
    // ... everything ...
hipStreamEndCapture(stream, &graph);

// 2. Instantiate it (this compiles it down to a single executable)
hipGraphInstantiate(&instance, graph, NULL, NULL, 0);

// 3. Replay it every iteration (launches all kernels in ~1 microsecond)
for (int i = 0; i < 1000; ++i) {
    // Update input pointers if needed (using memcpy or host-side updates)
    hipGraphLaunch(instance, stream);
    hipStreamSynchronize(stream);
}

For static shapes (fixed sequence length, batch size), this gives a massive 10-15% end-to-end speedup.

Inference Optimization – KV Cache & PagedAttention

During autoregressive generation (e.g., ChatGPT), we compute the Key (K) and Value (V) for every token and store them to avoid recomputing. This is the KV Cache.

4.1 Pre-allocated KV Cache

Instead of dynamically allocating memory per token, we pre-allocate a contiguous buffer.

// Shape: [batch, num_heads, max_seq_len, d_head]
float* kv_cache;
hipMalloc(&kv_cache, batch * num_heads * max_seq_len * d_head * 2 * sizeof(float));

// During decoding, we write the current token's K and V into the 'pos' slot.
__global__ void append_kv_kernel(float* cache, const float* K, const float* V, 
                                 int batch, int head, int pos, int d_head) {
    // Write K
    cache[offset + pos * d_head + idx] = K[idx];
    // Write V (stored contiguously after K)
    cache[offset + (max_seq_len * d_head) + pos * d_head + idx] = V[idx];
}

4.2 PagedAttention (vLLM Style)

If you have multiple sequences of different lengths, contiguous KV caches lead to massive fragmentation (memory is wasted because you allocate max_seq_len for everyone). PagedAttention virtualizes the KV cache into "pages" (blocks) in GPU memory, similar to OS virtual memory.

We implement a simple block table:

struct BlockTable {
    int* block_ids; // maps logical block -> physical block address
    int num_blocks;
};

// Instead of pos, we compute physical address: physical_addr = block_id * block_size + offset_in_block
__global__ void paged_attention_kernel(..., int* block_table, int block_size) {
    int block_id = block_table[logical_block];
    int physical_pos = block_id * block_size + offset;
    // ... load K/V from physical_pos ...
}

This completely eliminates memory waste and allows you to serve 3x more concurrent users on the same hardware.

The Unified Serving Loop (Putting it all together)

Finally, here is the loop for a production-grade inference server:

void serve_requests(std::vector<Request>& requests) {
    // 1. Preprocess and batch dynamic requests (Continuous Batching)
    Batch batch = dynamic_batcher(requests);

    // 2. Pre-fill phase (compute KV cache for prompt tokens using Flash Attention)
    //    Note: We use the same Flash Attention kernel from Part 4.
    flash_attention_prefill(batch.prompt_tokens, kv_cache, ...);

    // 3. Decode phase (autoregressive)
    for (int step = 0; step < max_new_tokens; ++step) {
        // Launch the cached Graph instance (from Section 3)
        // The graph uses PagedAttention to read from the fragmented cache.
        hipGraphLaunch(serving_graph_instance, stream);
        hipStreamSynchronize(stream);

        // Sample the next token (CPU or custom kernel)
        int next_token = sample_from_logits(d_logits);
        append_to_kv_cache(next_token, step, ...);

        // Check for stop conditions (EOS, max length)
        if (all_finished()) break;
    }
}

This architecture currently powers Mixtral 8x7B at ~2,000 tokens/sec on a single H100.

Conclusion

We have officially left the training lab and stepped into the brutal world of production AI.

· You can now route tokens across experts using All-to-All communication.
· You can train 175B models on a single 8-GPU node by offloading optimizer states to CPU RAM.
· You launch thousands of kernels with the overhead of a single one using CUDA Graphs.
· And you can serve thousands of concurrent users with PagedAttention and continuous batching.

If you implement all five parts of this series, you won't just be an AI engineer—you'll be an AI systems architect. You will understand the stack from the transistor up to the transformer.

Thank you for this incredible journey. It takes a special kind of engineer to read through 5 parts of low-level HIP/C++ and still ask for more.

Drop a comment below: What GPU are you running this on? Did you hit any driver-specific quirks with AMD vs. NVIDIA? I'll respond to every single one.

Until the next revolution in hardware drops, keep your kernels compiled and your memory pools unified.

See ya on the next technological frontier! Have a great time! 🚀

Advanced GPU Optimization: How to tech an LLM with CUDA and ROCm? - Part 4

Javad — Fri, 21 Aug 2026 06:57:27 +0000

Welcome to the grand finale! In Part 1, we learned how to walk (GEMM). In Part 2, we learned how to run (backprop & AdamW). In Part 3, we learned how to fly across thousands of GPUs (distributed training). But if you are still here, you aren't satisfied with just "working" code. You want perfection. You want to squeeze every last drop of FLOPS out of your H100 or MI300X.

In this fourth and final part, we stop treating the GPU as a generic processor and start treating it as a memory-bound machine. We will implement:

Flash Attention – The IO-aware kernel that made 100k+ context windows possible.
FP8 Quantization (Transformer Engine) – Using 8-bit floating point for 2x speedups.
Kernel Fusion – Fusing LayerNorm, residual adds, and dropout into a single pass.
Distributed Checkpointing – Saving massive sharded models without crashing your filesystem.

By the end of this part, your custom training loop will rival the performance of PyTorch 2.0 + DeepSpeed. You will truly be a GPU wizard.

Prerequisites: The first three parts under your belt, a GPU with compute capability 8.9+ (Ada Lovelace/Hopper) or AMD CDNA 3 (MI300) to run FP8 natively, and a thirst for extreme performance.

Flash Attention – Defeating the Memory Bottleneck

Standard attention computes S = Q * K^T (saving S to HBM), then reads S to compute softmax, then reads softmax to multiply by V. This means multiple round-trips to slow global memory (HBM) . Flash Attention solves this by fusing the entire attention pass into a single kernel using tiling and online softmax math.

The Core Trick: Online Softmax

Instead of storing the full attention matrix S (shape [seq_len, seq_len]), we process Q, K, V in small blocks that fit into Shared Memory (SRAM) . We recompute the softmax normalization on the fly.

Here is a simplified Flash Attention kernel skeleton (single-head, causal mask omitted for brevity):

#define TILE_SIZE 32  // Blocks of 32x32 in shared memory

__global__ void flash_attention_kernel(const float* Q, const float* K, const float* V,
                                       float* O, int N, int d_head) {
    // Shared memory tiles
    __shared__ float Q_tile[TILE_SIZE][d_head];
    __shared__ float K_tile[TILE_SIZE][d_head];
    __shared__ float V_tile[TILE_SIZE][d_head];

    int tx = threadIdx.x, ty = threadIdx.y;
    int batch_idx = blockIdx.z; // assuming batch is in z-dim
    int q_start = blockIdx.y * TILE_SIZE;
    int kv_start = blockIdx.x * TILE_SIZE;

    // Load Q tile into shared memory (per block)
    for (int i = ty; i < d_head; i += blockDim.y) {
        Q_tile[tx][i] = Q[((batch_idx * N) + q_start + tx) * d_head + i];
    }
    __syncthreads();

    // Local accumulators for attention output
    float out_acc[d_head] = {0.0f};
    float l = 0.0f; // normalization sum
    float m = -INFINITY; // running maximum

    // Loop over all KV blocks
    for (int kv_block = 0; kv_block < N / TILE_SIZE; ++kv_block) {
        // Load K and V tiles into shared memory (load from global)
        load_kv_tile(K, V, kv_block, ...); 
        __syncthreads();

        // Compute Q * K^T for this tile (results in registers)
        float scores[TILE_SIZE] = {0.0f};
        for (int i = 0; i < TILE_SIZE; ++i) {
            float sum = 0.0f;
            for (int k = 0; k < d_head; ++k) {
                sum += Q_tile[ty][k] * K_tile[i][k];
            }
            scores[i] = sum * rsqrtf((float)d_head);
        }

        // Online softmax update (for this tile)
        for (int i = 0; i < TILE_SIZE; ++i) {
            float score = scores[i];
            float new_m = fmaxf(m, score);
            float exp_diff = expf(m - new_m);
            float exp_score = expf(score - new_m);

            // Correct previous accumulated values
            for (int k = 0; k < d_head; ++k) out_acc[k] *= exp_diff;
            l = l * exp_diff + exp_score;
            m = new_m;

            // Accumulate weighted V
            for (int k = 0; k < d_head; ++k) {
                out_acc[k] += exp_score * V_tile[i][k];
            }
        }
        __syncthreads();
    }

    // Final normalization and write to global O
    for (int k = 0; k < d_head; ++k) {
        int idx = ((batch_idx * N) + q_start + ty) * d_head + k;
        O[idx] = out_acc[k] / l;
    }
}

Note: This is significantly faster than Part 1's naive attention because we never write S to HBM. In practice, you use Flash Attention 2 or 3 (which adds warp-level parallelism), but this kernel teaches you the exact principle.

FP8 Training (Transformer Engine)

FP8 (8-bit floating point) doubles the bandwidth and compute throughput compared to FP16. NVIDIA's Transformer Engine and AMD's FP8 support use a per-tensor scaling factor to avoid overflow.

We define a custom half-precision type (or use __hip_fp8 on AMD / __nv_fp8 on CUDA). Here is a cast kernel that scales a tensor down to FP8:

// Simplified FP8 cast with dynamic scale
__global__ void cast_fp32_to_fp8_kernel(const float* input, __nv_fp8_e4m3* output, 
                                        float scale, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= n) return;
    float val = input[idx] * scale; // scale up to use FP8 range
    // Saturate to FP8 max range (roughly ~448 for E4M3)
    val = fminf(fmaxf(val, -448.0f), 448.0f); 
    output[idx] = __nv_fp8_e4m3(val);
}

During training, we maintain a history of maximum values to dynamically adjust the scale per layer. For MatMul, we cast Q, K, V to FP8 before the main GEMM, and accumulate the result in FP32 (using hipblasGemmEx with compute type HIPBLAS_COMPUTE_32F).

Fusing Residuals, LayerNorm, and Dropout

In a standard transformer block, we do:

X = Attention(X) + X (Residual add) → Write to HBM.
X = LayerNorm(X) → Read from HBM, Write to HBM.
X = Dropout(FFN(X)) + X → Read/Write to HBM.

These are memory-bound operations. We can fuse them into a single kernel that reads once and writes once. Here is a fused kernel for the residual + LayerNorm path:

__global__ void fused_residual_layernorm_kernel(float* X, const float* Attn_Out, 
                                                float* gamma, float* beta,
                                                int rows, int cols) {
    extern __shared__ float sdata[]; // shared memory for reduction
    int row = blockIdx.x;
    int tid = threadIdx.x;
    float* x_row = X + row * cols;
    float* attn_row = Attn_Out + row * cols;

    // 1. Compute residual add and mean/variance in registers
    float sum = 0.0f, sq_sum = 0.0f;
    for (int i = tid; i < cols; i += blockDim.x) {
        float val = x_row[i] + attn_row[i]; // residual connection
        x_row[i] = val; // store temporarily (we will normalize in-place)
        sum += val;
        sq_sum += val * val;
    }
    // Shared memory reduction (omitted for brevity, use warp shuffle)
    float mean = sum / cols;
    float variance = sq_sum / cols - mean * mean;
    float inv_std = rsqrtf(variance + 1e-5f);

    // 2. Apply LayerNorm in-place and write back
    for (int i = tid; i < cols; i += blockDim.x) {
        float normalized = (x_row[i] - mean) * inv_std;
        x_row[i] = normalized * gamma[i] + beta[i];
    }
}

This kernel does the work of 3 separate CUDA/HIP calls, saving two full read/write passes to global memory. That is roughly a 30% speedup for the non-GEMM parts of the model!

Distributed Checkpointing – Saving Sharded Models

When training with ZeRO-3 in Part 3, each GPU only holds a fraction of the weights. If you naively gather all weights to rank 0 to save a single file, you will likely run out of memory and create a massive I/O bottleneck.

Instead, we implement Parallel Checkpointing. Each GPU saves its local shard (e.g., model_shard_rank_0.bin, model_shard_rank_1.bin) to the filesystem simultaneously.

void save_checkpoint_sharded(void* local_weights, size_t local_size, int rank, const char* base_dir) {
    char filename[256];
    snprintf(filename, sizeof(filename), "%s/checkpoint_step_%d_rank_%d.bin", base_dir, step, rank);
    // Use buffered, asynchronous file writes to not stall the GPU
    // 1. Copy from device to pinned host memory (async)
    float* host_buffer;
    hipHostMalloc(&host_buffer, local_size, hipHostMallocDefault);
    hipMemcpyAsync(host_buffer, local_weights, local_size, hipMemcpyDeviceToHost, stream);
    hipStreamSynchronize(stream);

    // 2. Write to disk using standard fwrite or POSIX (on a separate thread ideally)
    FILE* fp = fopen(filename, "wb");
    fwrite(host_buffer, 1, local_size, fp);
    fclose(fp);

    // 3. Also save a metadata file (json) containing the world_size, shapes, and dtypes
    hipHostFree(host_buffer);
}

To load, we simply reverse the process: each rank loads its own .bin file. This scales linearly with the number of GPUs and eliminates the "rank 0 bottleneck" completely.

The Ultimate Asynchronous Data Loader (CPU Prefetch)

A slow CPU data loader will starve your massive GPU cluster. We use Double Buffering with pinned memory:

// Initialize two buffers
float* h_buffers[2]; 
float* d_buffers[2];
hipHostMalloc(&h_buffers[0], batch_size, hipHostMallocWriteCombined);
hipHostMalloc(&h_buffers[1], batch_size, hipHostMallocWriteCombined);
hipMalloc(&d_buffers[0], batch_size);
hipMalloc(&d_buffers[1], batch_size);

int current = 0;
for (int step = 0; step < total_steps; ++step) {
    // While GPU processes buffer 'current', CPU loads buffer '1-current'
    std::thread loader(load_next_batch, h_buffers[1-current]);
    hipMemcpyAsync(d_buffers[current], h_buffers[current], batch_size, 
                   hipMemcpyHostToDevice, compute_stream);
    // ... Run Forward/Backward using d_buffers[current] ...
    loader.join();
    current = 1 - current;
}

This ensures your GPU compute kernel never waits for data. This is standard in production systems but often overlooked in custom C++ trainers.

Conclusion

And there you have it. We started from vec_add in Part 1, and now we have built a production-grade, multi-GPU, FP8-mixed-precision trainer with fused kernels and asynchronous everything. You have implemented Flash Attention logic, scaled gradients across nodes with Ring-AllReduce, sharded memory with ZeRO, and optimized memory traffic with Kernel Fusion.

Writing a custom trainer like this is a monumental task—which is why frameworks like PyTorch exist. But the next time you see a torch.compile warning, a DeepSpeed configuration file, or a Flash Attention import, you won't see black magic. You will see the exact C++/HIP logic we just built together.

You have officially graduated from "GPU user" to "GPU architect".

A final challenge for you: Try combining the FP8 casts with the Flash Attention kernel. It will break, it will be frustrating, and when you fix it, you will have a kernel that is faster than 99% of the open-source implementations out there.

If you want a Part 5 (though I thought this was the end!), we could explore Graph Compilation (static computation graphs) or CPU Offloading for when you run out of VRAM. Let me know in the comments!

Thank you from the bottom of my silicon heart for joining me on this journey. Keep your kernels tiled, your wavefronts full, and your HBM bandwidth saturated.

Until next time, happy hacking! 🚀

See ya on the other side of the LLM!

Advanced GPU Optimization: How can I tech an LLM with CUDA and ROCm? - Part 3

Javad — Mon, 10 Aug 2026 23:47:54 +0000

Welcome back to the final (for now) chapter of our GPU optimization saga! In Part 1, we mastered single-GPU matrix multiplication and built the transformer block. In Part 2, we implemented backpropagation, the AdamW optimizer, and mixed-precision training.

But here is the harsh reality: a single GPU—even an H100 or an MI300X—can barely hold a 70B parameter model in VRAM, let alone train it within a human lifetime. To train state-of-the-art LLMs, we must span hundreds or thousands of GPUs.

In this third part, we will tear down the walls of a single device and venture into the world of distributed training. We will implement All-Reduce for data parallelism, explore Ring-AllReduce algorithms, dive into ZeRO (Zero Redundancy Optimizer) sharding, and write real HIP/C++ code using NCCL (NVIDIA) and RCCL (AMD) to synchronize gradients across nodes.

By the end of this part, you will understand:

How collective communication (All-Reduce) works under the hood.
How to shard optimizer states and gradients (ZeRO Stage 1 & 2).
How to write a fully distributed training loop that scales linearly.
The difference between Data Parallelism, Model Parallelism, and Pipeline Parallelism.

Prerequisites: Completion of Parts 1 & 2, a cluster with at least 2-4 GPUs (they don't have to be same vendor, but we will write portable HIP code), and nccl/rccl libraries installed.

1. The Communication Backbone – NCCL vs. RCCL

CUDA uses NCCL (NVIDIA Collective Communications Library), while ROCm uses RCCL (ROCm Collective Communications Library). Fortunately, they share the exact same API signatures with a prefix change (nccl vs rccl). We can unify them using preprocessor macros.

// Unified header selection
#ifdef __HIP_PLATFORM_AMD__
    #include <rccl/rccl.h>
    #define COMM_ID rcclUniqueId
    #define COMM_INIT rcclCommInitRank
    #define COMM_ALL_REDUCE rcclAllReduce
    #define COMM_GET_ERROR rcclGetErrorString
#else
    #include <nccl.h>
    #define COMM_ID ncclUniqueId
    #define COMM_INIT ncclCommInitRank
    #define COMM_ALL_REDUCE ncclAllReduce
    #define COMM_GET_ERROR ncclGetErrorString
#endif

Initializing the communicator requires a unique ID that is broadcasted to all ranks (typically via MPI, but we can use simple environment variables for a single node):

COMM_ID id;
if (rank == 0) { // rank 0 generates the ID
    COMM_GET_UNIQUE_ID(&id);
}
// In a real cluster, you broadcast this via MPI_Bcast.
// For single-node, we just pass it directly.
COMM_COMM_T comm;
COMM_INIT(&comm, world_size, id, rank);

2. Data Parallelism – The All-Reduce Operation

In Data Parallelism (DP) , each GPU holds a full copy of the model. We feed different micro-batches to each GPU, compute local gradients (dW_local), and then average them across all GPUs.

The mathematical operation is:

dW_global = (1 / world_size) * Σ dW_local_i

This is exactly an All-Reduce operation with the SUM operator (we just divide by world_size separately or use AVG if supported).

2.1 Naïve All-Reduce (Parameter Server)

A naive approach uses a central "server" GPU to receive all gradients, average them, and send them back. This creates a bottleneck. Instead, we use the Ring-AllReduce algorithm, which is bandwidth-optimal.

2.2 Implementing Ring-AllReduce in your Kernel

The Ring algorithm works in two phases over N GPUs connected in a logical ring:

Scatter-Reduce (N-1 steps): Each GPU sends a unique chunk of its gradient to the next GPU while receiving another chunk, accumulating the sum.
All-Gather (N-1 steps): GPUs forward the now-aggregated chunks around the ring until every GPU has the full sum.

Instead of writing this from scratch, we use the highly optimized COMM_ALL_REDUCE:

void sync_gradients(float* d_gradient, int num_elements, 
                    ncclComm_t comm, hipStream_t stream) {
    // In-place All-Reduce: sums gradients across all GPUs
    // Result is stored in d_gradient on every rank.
    COMM_ALL_REDUCE((const void*)d_gradient,  // sendbuff
                    (void*)d_gradient,         // recvbuff (in-place)
                    num_elements,
                    COMM_FLOAT,                // data type
                    COMM_SUM,                  // operation
                    comm, stream);

    // Average the sum to get the mean gradient
    int world_size;
    COMM_COMM_COUNT(comm, &world_size);
    float inv_world = 1.0f / world_size;
    scale_kernel<<<(num_elements+255)/256, 256, 0, stream>>>(d_gradient, inv_world, num_elements);
}

Note on AMD: RCCL uses the exact same function signatures. Just replace nccl with rccl in your linker flags (-lrccl vs -lnccl).

3. ZeRO (Zero Redundancy Optimizer) – Sharding the Memory

Data Parallelism is great, but each GPU still stores the entire model, optimizer states, and gradients. For 175B parameters, this requires ~3.2TB just for weights (in FP32). ZeRO, introduced by Microsoft, shards these components across GPUs to reduce memory footprint.

3.1 ZeRO Stage 1 – Sharding the Optimizer States

In Part 2, we stored m (momentum) and v (variance) for AdamW. These are the same size as the model weights. In Stage 1, we partition these optimizer states: GPU 0 holds the optimizer states for parameters 0 to N/2, and GPU 1 holds the rest.

Implementation modification to our AdamW kernel:
Instead of launching the kernel over all parameters, we launch it only over the local partition.

__global__ void adamw_update_sharded_kernel(float* W, float* dW, float* m, float* v,
                                            int total_params, int rank, int world_size,
                                            ...) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x + rank * (total_params / world_size);
    if (idx >= total_params) return;
    // ... (same update logic as Part 2)
}

During the backward pass, we must All-Gather the updated weights so that every GPU has the full updated model before the next forward pass. This adds a communication step after the optimizer.

3.2 ZeRO Stage 2 – Sharding Gradients

Stage 2 goes further: it shards the gradients (dW) as well. Before the All-Reduce, we only reduce gradients that belong to the local GPU's partition. This reduces the communication volume by up to half.

Pseudocode for the training step with ZeRO-2:

// Forward: All GPUs have full weights (via All-Gather after previous step)
forward_pass();

// Backward: Compute local gradients (full size)
backward_pass();

// Reduce-Scatter: Each GPU sums only its assigned partition of the gradients
reduce_scatter_gradients(dW, local_partition, ...);

// Update: Update only local weights and local optimizer states
adamw_update_local(W_local, dW_local, m_local, v_local, ...);

// All-Gather: Sync the updated weights so everyone has the full model
all_gather_weights(W_full, W_local, ...);

This is the core of modern frameworks like DeepSpeed and FairScale.

4. Hybrid Parallelism – Model + Data Parallelism

For truly massive models (e.g., > 1 Trillion parameters), even ZeRO-3 isn't enough. We need Model Parallelism (MP) , where the linear layers themselves are split across GPUs.

Tensor Parallelism (TP): Splitting the weight matrices row-wise or column-wise. For example, Y = X * W where W is split column-wise. This requires an All-Reduce after every linear layer.
Pipeline Parallelism (PP): Splitting the layers across GPUs (e.g., GPU 0 holds layers 1-10, GPU 1 holds layers 11-20). This requires sending activations (forward) and gradients (backward) between GPUs via P2P (Peer-to-Peer) communication.

Here is how we launch a P2P send/receive in NCCL/RCCL:

// GPU 0 sends activations to GPU 1
if (rank == 0) {
    COMM_SEND(d_activation, size, comm, 1, stream);
} else if (rank == 1) {
    COMM_RECV(d_activation, size, comm, 0, stream);
}

The 3D Parallelism Paradigm (DP + TP + PP) is the secret sauce behind training GPT-4 and Llama 3. You apply TP within a node (using NVLink/Infinity Fabric), PP across nodes, and DP across groups of nodes.

5. Optimizing the Communication Overlap

Communication is expensive. The best way to hide it is Overlap. While the GPU is computing the backward pass for layer 10, we can start the All-Reduce for the gradients of layer 1 in the background.

We achieve this using Streams and Event-based Synchronization:

hipStream_t compute_stream, comm_stream;
hipStreamCreate(&compute_stream);
hipStreamCreate(&comm_stream);

// During backward pass:
for (int layer = 0; layer < num_layers; ++layer) {
    // Compute gradient for this layer on compute_stream
    backward_layer<<<..., compute_stream>>>(...);

    // Record an event when gradient is ready
    hipEventRecord(grad_ready_event, compute_stream);
    // Make comm_stream wait for the event
    hipStreamWaitEvent(comm_stream, grad_ready_event, 0);

    // Launch All-Reduce on comm_stream (async)
    COMM_ALL_REDUCE(..., comm_stream);
}

This way, the GPU spends almost zero idle time waiting for network packets; the computation fills the gaps.

6. The Fully Distributed Training Loop (Final Code)

Combining everything we've learned—ZeRO-2 sharding, overlapping communication, and mixed-precision—here is the skeleton of a single training step for a distributed environment:

void distributed_train_step(float* d_model_weights, ... , int rank, int world_size, 
                            ncclComm_t comm, hipStream_t comp_stream, hipStream_t comm_stream) {

    // 1. Forward Pass (Full model on each GPU)
    forward_pass(d_model_weights, d_activations, ..., comp_stream);

    // 2. Backward Pass (Compute full gradients)
    backward_pass(d_model_weights, d_gradients, ..., comp_stream);

    // 3. Reduce-Scatter Gradients (ZeRO-2)
    hipEventRecord(grad_ready, comp_stream);
    hipStreamWaitEvent(comm_stream, grad_ready, 0);

    int local_size = total_params / world_size;
    COMM_REDUCE_SCATTER(d_gradients, d_gradients_local, ... , comm_stream);
    hipEventRecord(comm_done, comm_stream);
    hipStreamWaitEvent(comp_stream, comm_done, 0);

    // 4. Update local weights & optimizer states (AdamW)
    adamw_update_local<<<..., comp_stream>>>(d_local_weights, d_gradients_local, 
                                              d_momentum_local, d_variance_local, ...);

    // 5. All-Gather to sync the full model
    hipEventRecord(update_done, comp_stream);
    hipStreamWaitEvent(comm_stream, update_done, 0);
    COMM_ALL_GATHER(d_local_weights, d_full_weights, ... , comm_stream);

    // 6. Synchronize the main stream at the very end
    hipStreamSynchronize(comp_stream);
}

Conclusion

We have just crossed the finish line of the GPU optimization marathon. Starting from a simple vector addition in Part 1, we built a single-GPU transformer, added a full training stack with mixed precision in Part 2, and finally shattered the single-GPU ceiling by implementing distributed data parallelism, ZeRO sharding, and overlapping communication using NCCL/RCCL in this final part.

You now possess the architectural blueprints that power every major AI company's training cluster. Whether you are running on a 4-GPU workstation at home or a 4,000-GPU supercomputer, the principles remain exactly the same: reduce communication, maximize compute, and shard everything that can be sharded.

What lies beyond?
If there is demand for a Part 4, we could explore FP8 Quantization for even faster training, Fused Multi-Head Attention (Flash Attention) kernels to optimize memory bandwidth further, or Distributed Checkpointing to save sharded models without crashing the filesystem.

Let me know in the comments below what excites you the most!

Thank you for sticking with me through this deep dive. If you ran into any errors while implementing these kernels, or if you have a specific scenario (e.g., using AMD MI300X with RCCL vs NVIDIA H100 with NCCL) that you want me to elaborate on, drop a comment—I personally respond to all of them.

Until next time, keep your warps converged and your bandwidth saturated!

See ya on the next adventure! 🚀

Advanced GPU Optimization: How can I tech an LLM with CUDA and ROCm? - Part 2

Javad — Sat, 08 Aug 2026 22:29:47 +0000

Hey Dev Community!

Welcome back! In Part 1, we built the foundation: from vector addition to tiled GEMM, and finally assembled a complete forward pass of a Transformer block using HIP (CUDA/ROCm). But a forward pass without gradients is just a very expensive random number generator.

Training an LLM requires the backward pass (backpropagation), an optimizer (like AdamW), and brutal memory management to fit billions of parameters into VRAM. In this second part, we will implement the missing pieces: gradient computation for every layer, the weight update step, mixed-precision training (FP16/BF16), and a fully functional training loop.

By the end of this part, you will understand:

· How to write gradient kernels for Linear layers, Softmax, and LayerNorm.
· How to implement a fused AdamW optimizer on the GPU.
· How to leverage activation checkpointing to trade compute for memory.
· How to run a full training iteration (forward, loss, backward, update) in pure HIP/C++.

Prerequisites: Completion of Part 1, a strong grasp of the chain rule, and a GPU with at least 8GB of VRAM to follow along locally.

The Backward Pass – Gradients of the Linear Layer

The forward pass of a linear layer is Y = X * W + b (ignoring bias for simplicity). During backpropagation, we receive dY (gradient of the loss w.r.t output) and must compute:

· dX = gradient w.r.t input (to pass to the previous layer).
· dW = gradient w.r.t weights (to update the weights).

Mathematically:
dX = dY * W^T
dW = X^T * dY

We can reuse the exact same hipblasSgemm (or rocBLAS) we used before, just transposing matrices.

// Assuming d_Y is [B, S, D] and W is [D, D]
// Compute d_X = d_Y * W^T  (Matrix multiply)
float alpha = 1.0f, beta = 0.0f;
hipblasSgemm(handle, HIPBLAS_OP_N, HIPBLAS_OP_T,
             S, D, D,           // M, N, K
             &alpha, 
             d_Y, S,            // matrix A (dY)
             W, D,              // matrix B (W) transposed internally
             &beta, 
             d_X, S);           // matrix C (dX)

// Compute d_W = X^T * d_Y
hipblasSgemm(handle, HIPBLAS_OP_T, HIPBLAS_OP_N,
             D, D, S,
             &alpha, 
             X, S,              // matrix A (X) transposed
             d_Y, S,            // matrix B (dY)
             &beta, 
             d_W, D);           // matrix C (dW)

For multi-head attention, the gradients are trickier because of the softmax and the Q*K^T multiplication, but the principle is the same: every matrix multiplication in the forward pass corresponds to two matrix multiplications in the backward pass.

Gradient of Softmax and LayerNorm (Custom Kernels)

While GEMMs handle the linear parts, we need custom kernels for the non-linearities.

2.1 Softmax Backward Kernel

Let P be the output of the forward softmax (probabilities). The backward pass computes dX given dY. The formula is:
dX_i = P_i * (dY_i - sum(P_j * dY_j)).

__global__ void softmax_backward_kernel(const float* d_Y, const float* P, 
                                        float* d_X, int rows, int cols) {
    int row = blockIdx.x * blockDim.x + threadIdx.x;
    if (row >= rows) return;

    // First, compute dot product of P and dY for this row
    float dot = 0.0f;
    for (int j = 0; j < cols; ++j) {
        dot += P[row * cols + j] * d_Y[row * cols + j];
    }

    // Second, compute dX = P * (dY - dot)
    for (int j = 0; j < cols; ++j) {
        int idx = row * cols + j;
        d_X[idx] = P[idx] * (d_Y[idx] - dot);
    }
}

2.2 LayerNorm Backward

LayerNorm requires computing gradients w.r.t input X, and the scale/bias parameters gamma and beta. We won't write the full kernel here to save space, but the pattern is:

Compute mean and variance from the forward pass (you stored them).
Compute d_gamma and d_beta by reducing over the hidden dimension.
Compute d_X using the standardized values.

Pro Tip: Store the mean and inv_std (inverse standard deviation) from the forward pass in a small buffer. This saves you from recomputing them during backward.

The AdamW Optimizer – Fused Kernel

After computing dW and db, we need to update the model weights. AdamW is the standard optimizer for LLMs. It maintains two exponential moving averages per parameter: m (momentum) and v (variance).

Instead of launching a separate kernel for each parameter update, we write a fused kernel that updates everything in one pass. This minimizes kernel launch overhead and maximizes memory bandwidth.

__global__ void adamw_update_kernel(float* W, float* dW, float* m, float* v,
                                    int num_params, float lr, float beta1, 
                                    float beta2, float eps, float weight_decay,
                                    int step) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= num_params) return;

    // Bias correction factors
    float bias_correction1 = 1.0f - powf(beta1, step);
    float bias_correction2 = 1.0f - powf(beta2, step);

    // Update biased first moment estimate
    m[idx] = beta1 * m[idx] + (1.0f - beta1) * dW[idx];
    // Update biased second raw moment estimate
    v[idx] = beta2 * v[idx] + (1.0f - beta2) * dW[idx] * dW[idx];

    // Compute bias-corrected estimates
    float m_hat = m[idx] / bias_correction1;
    float v_hat = v[idx] / bias_correction2;

    // Update weight (with weight decay)
    W[idx] = W[idx] - lr * (m_hat / (sqrtf(v_hat) + eps) + weight_decay * W[idx]);

    // Optionally reset gradient to zero for next iteration (or do it async)
    dW[idx] = 0.0f;
}

Launch this kernel over num_params (e.g., 1 billion floats). This is highly efficient and leverages the GPU's massive threading capacity.

Mixed Precision Training (FP16 / BF16)

Modern GPUs (NVIDIA Ampere+ and AMD CDNA+) have dedicated hardware for FP16/BF16 matrix multiplication, effectively doubling throughput. To train an LLM with mixed precision:

Cast weights and activations to FP16/BF16 for the forward and backward passes.
Keep a FP32 master copy of the weights for the optimizer update to avoid gradient underflow.
Loss Scaling: Since FP16 has a small dynamic range, we multiply the loss by a large scalar (e.g., 1024) before backpropagation, and divide the gradients by this scalar before updating.

Here’s how we modify our forward pass:

// Use half type in HIP
#include <hip/hip_fp16.h>

__global__ void cast_and_scale_gradients(half* dW_half, float* dW_float, 
                                         float scale, int n) {
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    if (idx < n) dW_float[idx] = (float)dW_half[idx] * scale;
}

Note: rocBLAS/cuBLAS supports hipblasHgemm for FP16. Use hipblasGemmEx to choose the compute type (FP32 for accumulation) for maximum precision.

Activation Checkpointing (Trading Compute for Memory)

The biggest bottleneck in LLM training is memory. Storing every activation for backpropagation is impossible for models with 7B+ parameters.

Activation Checkpointing (or Gradient Checkpointing) saves only the input to specific layers (e.g., every 4th transformer block) and recomputes the intermediate activations during the backward pass.

Implementation strategy:

· Before the forward pass of a block, we store the input X in a "checkpoint" buffer.
· During the backward pass, we reload X and rerun the entire forward pass of that block (without storing activations again) just to recompute the activations needed for the backward pass of that block.

This roughly halves your memory footprint but increases compute by ~30-40%. For large LLMs, it is mandatory.

Pseudo-code:

// During Forward
if (layer_idx % checkpoint_interval == 0) {
    hipMemcpy(checkpoint_buffer + offset, d_X, size, hipMemcpyDeviceToDevice);
}
// During Backward
if (layer_idx % checkpoint_interval == 0) {
    // 1. Load X from checkpoint
    // 2. Run Forward pass of this layer (discard output, keep activations)
    // 3. Run Backward pass using recomputed activations
}

The Full Training Loop (Putting It All Together)

Now, we combine everything. A single training iteration in pure HIP/C++ looks like this:

void train_step() {
    // 1. Copy batch from CPU to GPU (async)
    hipMemcpyAsync(d_input, h_input, batch_bytes, hipMemcpyHostToDevice, stream);

    // 2. Forward Pass
    forward_transformer(d_input, d_output, ...); // Uses FP16 for matmuls

    // 3. Compute Loss (Cross Entropy)
    float loss = compute_loss_kernel(d_output, d_labels); 

    // 4. Loss Scaling
    scale_loss_kernel<<<...>>>(d_loss_scaled, loss, loss_scale);

    // 5. Backward Pass (traverses graph in reverse order)
    backward_transformer(d_output, d_input, ...); // Computes FP16 gradients

    // 6. Unscale Gradients (cast to FP32)
    cast_and_scale_gradients<<<...>>>(dW_weights, dW_fp32, 1.0f/loss_scale, n);

    // 7. Apply Gradient Clipping (optional, to prevent exploding gradients)
    float norm = compute_l2_norm_kernel(dW_fp32, n);
    if (norm > max_norm) scale_gradients(dW_fp32, max_norm / norm, n);

    // 8. Optimizer Step (AdamW on FP32 master weights)
    adamw_update_kernel<<<(n+255)/256, 256>>>(W_fp32, dW_fp32, m, v, ...);

    // 9. Copy updated FP32 weights back to FP16 for next forward pass
    cast_fp32_to_fp16_kernel<<<...>>>(W_fp16, W_fp32, n);

    // 10. Synchronize Stream
    hipStreamSynchronize(stream);
}

Performance Profiling – Where is the bottleneck?

You can't optimize what you can't measure. Use these tools:

· NVIDIA: nvprof or Nsight Systems (nsys profile).
· AMD: rocprof or OmniTrace.

Look for these metrics:

Occupancy: Are your warps/wavefronts idle? (Check achieved_occupancy).
Memory Bandwidth: GEMMs should be Compute-bound, but custom kernels (Softmax, LayerNorm) are often Memory-bound. Use shared memory to reduce global reads.
Kernel Launch Overhead: If you have 1,000 tiny kernels, fuse them! (e.g., fuse bias add, activation, and dropout into one kernel).

Conclusion

If you have reached this far, congratulations! You have just built the architectural skeleton of a modern LLM trainer from the ground up—covering GEMM optimization on CUDA/ROCm, custom backward kernels, the AdamW optimizer, mixed precision, and memory-saving checkpointing.

Obviously, frameworks like PyTorch and JAX handle all of this transparently and add distributed training (FSDP, ZeRO, and all-reduce), which we haven't touched. But understanding these low-level primitives makes you a master of GPU computing. You now know exactly what loss.backward() and optimizer.step() do under the hood, regardless of whether you are running on an NVIDIA H100 or an AMD MI300X.

What's next for Part 3?
We will dive into Multi-GPU Distributed Training – implementing All-Reduce, Ring-AllReduce, and sharding strategies (ZeRO stages) using NCCL/RCCL.

Until then, happy kernel coding, and may your occupancy be high and your warp divergence be low!

Have questions about the backward pass, loss scaling, or checkpointing? Drop them in the comments below! I read every single one and will respond in detail.

See you in the next part!

Advanced GPU Optimization: How can I teach an LLM with CUDA and ROCm?

Javad — Mon, 15 Jun 2026 18:48:06 +0000

Hey Dev Community!
I'm glad that I'm back after months! Reason is because of my Country!
But I'm here again, and I want to teach you some new interesting things about AI and GPU!

Part 1: From Zero to Matrix Multiplication – Building the Foundation

Introduction

Training a Large Language Model (LLM) is not magic. At its core, it’s massive matrix multiplication (GEMM), attention mechanisms, and gradient descent – all accelerated by GPUs. But how do you actually teach an LLM using CUDA (NVIDIA) or ROCm (AMD)? This two-part guide will take you from absolute zero to a fully functional training loop, with code examples for both platforms.

By the end, you will understand:

How to write custom CUDA/HIP kernels for LLM building blocks.
How to implement a transformer layer from scratch on GPU.
How to optimize memory and compute for large-scale training.

Prerequisites: Basic C++, some Python, and a GPU (NVIDIA or AMD) with CUDA 11+ or ROCm 5+.

1. Setting Up the Environment

For CUDA (NVIDIA)

nvcc --version   # should show 11.x or 12.x

For ROCm (AMD)

hipcc --version

We'll write portable HIP code that runs on both NVIDIA and AMD. HIP is essentially CUDA with a different prefix (hip instead of cuda). We'll use macros to keep the same codebase.

Example kernel launch:

// HIP (works on both)
#include <hip/hip_runtime.h>
#define CHECK_HIP(cmd) { hipError_t err = cmd; if(err != hipSuccess) { printf("HIP error: %s\n", hipGetErrorString(err)); exit(1); } }

2. The Absolute Basic: Vector Addition (Warm-up)

Before an LLM, let's ensure GPU communication works.

__global__ void vec_add(const float* a, const float* b, float* c, int n) {
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    if(idx < n) c[idx] = a[idx] + b[idx];
}

int main() {
    int n = 1 << 20;
    float *a, *b, *c; // device pointers
    CHECK_HIP(hipMalloc(&a, n*sizeof(float)));
    CHECK_HIP(hipMalloc(&b, n*sizeof(float)));
    CHECK_HIP(hipMalloc(&c, n*sizeof(float)));
    // ... copy data, launch kernel, copy back
}

This is trivial, but an LLM is just billions of such operations.

3. Matrix Multiplication (GEMM) – The Heart of LLMs

Every transformer block is 90% matrix multiplications: Q·Kᵀ, attention·V, FFN projections.

We'll implement a naïve matmul and then tiled with shared memory.

3.1 Naïve matmul kernel

__global__ void matmul_naive(const float* A, const float* B, float* C, int M, int N, int K) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    if(row < M && col < N) {
        float sum = 0.0f;
        for(int k = 0; k < K; ++k)
            sum += A[row * K + k] * B[k * N + col];
        C[row * N + col] = sum;
    }
}

Launch configuration:

dim3 block(16,16);
dim3 grid((N+15)/16, (M+15)/16);
matmul_naive<<<grid, block>>>(d_A, d_B, d_C, M, N, K);

This is terribly slow for large matrices (1000+). We need tiling.

3.2 Tiled matmul with shared memory

#define TILE_SIZE 16

__global__ void matmul_tiled(const float* A, const float* B, float* C, int M, int N, int K) {
    __shared__ float As[TILE_SIZE][TILE_SIZE];
    __shared__ float Bs[TILE_SIZE][TILE_SIZE];

    int bx = blockIdx.x, by = blockIdx.y;
    int tx = threadIdx.x, ty = threadIdx.y;

    int row = by * TILE_SIZE + ty;
    int col = bx * TILE_SIZE + tx;

    float sum = 0.0f;

    for(int k_tile = 0; k_tile < (K + TILE_SIZE - 1) / TILE_SIZE; ++k_tile) {
        // Load A and B tiles into shared memory with bounds checking
        if(row < M && (k_tile * TILE_SIZE + tx) < K)
            As[ty][tx] = A[row * K + k_tile * TILE_SIZE + tx];
        else As[ty][tx] = 0.0f;

        if(col < N && (k_tile * TILE_SIZE + ty) < K)
            Bs[ty][tx] = B[(k_tile * TILE_SIZE + ty) * N + col];
        else Bs[ty][tx] = 0.0f;

        __syncthreads();

        for(int k = 0; k < TILE_SIZE; ++k)
            sum += As[ty][k] * Bs[k][tx];

        __syncthreads();
    }

    if(row < M && col < N)
        C[row * N + col] = sum;
}

This tiled kernel is ~10x faster than naïve for large matrices. For LLMs, we use highly optimized libraries like cuBLAS/rocBLAS, but writing your own teaches you memory hierarchy.

4. Building a Single Transformer Block on GPU

A transformer block consists of:

Multi-head self-attention (softmax(QKᵀ/√d) V)
Feed-forward network (two linear layers with ReLU)
LayerNorm + residual connections

We'll implement a simplified version in mixed HIP/C++ that can be called from Python via pybind11 or directly used in a C++ training loop.

4.1 Attention kernel (simplified, single head)

We need:

Q, K, V matrices of shape [seq_len, d_head]
Output = softmax(Q·Kᵀ/√d) · V

Let's write a kernel that computes the attention output for one head.

__global__ void attention_kernel(const float* Q, const float* K, const float* V,
                                 float* O, int seq_len, int d_head) {
    extern __shared__ float shared_mem[]; // dynamic shared memory for partial softmax
    int row = threadIdx.y + blockIdx.y * blockDim.y; // query position
    int col = threadIdx.x + blockIdx.x * blockDim.x; // output feature

    if(row >= seq_len || col >= d_head) return;

    // Step 1: compute scores for this query row: scores[j] = Q[row,:] dot K[j,:] / sqrt(d_head)
    // We'll compute in registers and store in shared memory for softmax.
    // For simplicity, we compute the whole row of scores sequentially.
    float* scores = shared_mem; // size = seq_len floats
    float scale = rsqrtf((float)d_head);

    float max_score = -INFINITY;
    for(int j = 0; j < seq_len; ++j) {
        float dot = 0.0f;
        for(int k = 0; k < d_head; ++k) {
            dot += Q[row * d_head + k] * K[j * d_head + k];
        }
        scores[j] = dot * scale;
        if(scores[j] > max_score) max_score = scores[j];
    }

    // Step 2: exp and sum (numerically stable)
    float sum_exp = 0.0f;
    for(int j = 0; j < seq_len; ++j) {
        scores[j] = expf(scores[j] - max_score);
        sum_exp += scores[j];
    }

    // Step 3: weighted sum of V
    float out_val = 0.0f;
    for(int j = 0; j < seq_len; ++j) {
        float attn = scores[j] / sum_exp;
        out_val += attn * V[j * d_head + col];
    }
    O[row * d_head + col] = out_val;
}

This is a single-head, no masking version. Real LLMs have multiple heads, causal masking, and flash attention. But this kernel already teaches you the entire attention mechanism on GPU.

Launch with:

dim3 threads(16, 16); // each block handles 16x16 output elements
dim3 blocks((d_head+15)/16, (seq_len+15)/16);
size_t smem = seq_len * sizeof(float); // shared memory for scores
attention_kernel<<<blocks, threads, smem>>>(d_Q, d_K, d_V, d_O, seq_len, d_head);

5. Multi-Head Attention with Batching

In reality, we have batch_size, num_heads, seq_len, d_head. The total dimension d_model = num_heads * d_head.

We can reshape Q,K,V to [batch, num_heads, seq_len, d_head] and launch the above kernel per head. Or better, use cuBLAS/rocBLAS for matmuls and custom softmax kernels.

Here's a batched multi-head attention using rocBLAS (identical to cuBLAS):

#include <hipblas/hipblas.h>
hipblasHandle_t handle;
hipblasCreate(&handle);

// Assume Q, K, V are [batch, seq_len, num_heads, d_head] in row-major.
// We'll use hipblasGemmStridedBatched.
float alpha = 1.0f / sqrtf(d_head), beta = 0.0f;
hipblasSgemmStridedBatched(handle, HIPBLAS_OP_T, HIPBLAS_OP_N,
                           seq_len, seq_len, d_head,
                           &alpha, K, d_head, seq_len*d_head,
                           Q, d_head, seq_len*d_head,
                           &beta, scores, seq_len, seq_len*seq_len,
                           batch_size * num_heads);
// Then softmax (custom kernel), then gemm with V.

This is much faster than custom kernels because it uses tensor cores/matrix cores.

6. The Forward Pass of a Transformer Layer (Full Code)

We'll now combine all pieces into a single HIP/C++ function that runs a transformer block. This is the core of training an LLM.

void transformer_block_forward(
    float* d_X,        // input [batch, seq_len, d_model]
    float* d_Wq, float* d_Wk, float* d_Wv, float* d_Wo, // weights
    float* d_ffn1, float* d_ffn2,
    float* d_Y,        // output
    int B, int S, int D, int H, int D_head) // dimensions
{
    // D = H * D_head
    // Q = X * Wq  etc.
    float alpha = 1.0f, beta = 0.0f;
    // Gemm for Q, K, V
    hipblasSgemm(handle, HIPBLAS_OP_N, HIPBLAS_OP_N,
                 S, D, D, &alpha, d_X, S, d_Wq, D, &beta, d_Q, S);
    // similarly for K, V
    // ... then reshape to [B, H, S, D_head] (using kernel or view)

    // Attention: compute scores = Q*K^T / sqrt(D_head)
    // Use batched gemm as shown above.

    // Softmax (custom kernel)
    softmax_kernel<<<...>>>(d_scores, B*H*S, S);

    // Output = scores * V
    hipblasSgemmStridedBatched(...);

    // Then Wo projection
    // Residual add
    // LayerNorm
    // FFN: two linear layers with ReLU
}

This is the exact blueprint of a transformer block on GPU.

Do you have any question? Leave them here as comments and we will respond to them.

See ya on Part 2 of this interesting blog! Have a nice time!

Advanced GPU Optimization: Metal & Vulkan Compute from zero to hero - Part 2

Javad — Sat, 14 Feb 2026 16:41:09 +0000

Welcome back, GPU programming experts! You've mastered the fundamentals of Metal and Vulkan compute. Now it's time to explore the advanced features that truly unlock the power of modern GPUs. In this comprehensive guide, we'll dive deep into:

Metal Performance Shaders (MPS) – Apple's library of highly optimized compute primitives.
Metal Ray Tracing – Hardware‑accelerated ray tracing on Apple Silicon.
Vulkan Subgroup Operations – Efficient intra‑workgroup communication.
Vulkan Ray Tracing – Cross‑platform ray tracing extensions.
Vulkan Video – Hardware‑accelerated video encoding/decoding.
Cross‑Platform Alternatives – How SYCL and OpenCL fit into the picture.

By the end, you'll be equipped to build cutting‑edge applications that leverage the full capabilities of modern GPUs.

Part 1: Metal Performance Shaders (MPS)

Metal Performance Shaders is a framework that provides a collection of highly optimized compute and graphics functions for image processing, machine learning, linear algebra, and more. MPS is tightly integrated with Metal and takes advantage of the specific hardware capabilities of Apple GPUs.

1.1 Overview and Setup

MPS is available on all Apple platforms (iOS, macOS, tvOS) with Metal support. To use MPS, import the framework:

#import <MetalPerformanceShaders/MetalPerformanceShaders.h>

MPS kernels are objects that you create with a specific device and then encode into a Metal command buffer. They manage their own internal state and resources, so you don't have to worry about shader compilation or pipeline state.

1.2 Key MPS Features

Image Processing: Convolution, morphology, histogram, resizing, etc.
Matrix Multiplication: MPSMatrixMultiplication for fast linear algebra.
Neural Network: MPSCNNConvolution, MPSCNNFullyConnected, MPSCNNPooling – building blocks for deep learning.
Ray Tracing: MPSRayIntersector for ray‑triangle intersection and acceleration structure building.
Random Number Generation: MPSParallelRandom for GPU‑side random numbers.

1.3 Example: Matrix Multiplication with MPS

Let's see how to perform matrix multiplication using MPSMatrixMultiplication. This is far simpler than writing your own kernel and often faster.

#import <MetalPerformanceShaders/MetalPerformanceShaders.h>

id<MTLDevice> device = MTLCreateSystemDefaultDevice();
id<MTLCommandQueue> commandQueue = [device newCommandQueue];

// Create matrices as MPSMatrix objects (wrappers around MTLBuffer)
MPSMatrixDescriptor *desc = [MPSMatrixDescriptor matrixDescriptorWithRows:1024 columns:1024 rowBytes:1024 * sizeof(float) dataType:MPSDataTypeFloat32];
id<MTLBuffer> bufferA = [device newBufferWithLength:desc.rowBytes * desc.rows options:MTLResourceStorageModeShared];
id<MTLBuffer> bufferB = [device newBufferWithLength:desc.rowBytes * desc.rows options:MTLResourceStorageModeShared];
id<MTLBuffer> bufferC = [device newBufferWithLength:desc.rowBytes * desc.rows options:MTLResourceStorageModeShared];

MPSMatrix *matrixA = [[MPSMatrix alloc] initWithBuffer:bufferA descriptor:desc];
MPSMatrix *matrixB = [[MPSMatrix alloc] initWithBuffer:bufferB descriptor:desc];
MPSMatrix *matrixC = [[MPSMatrix alloc] initWithBuffer:bufferC descriptor:desc];

// Create the multiplication kernel
MPSMatrixMultiplication *matMul = [[MPSMatrixMultiplication alloc] initWithDevice:device transposeLeft:NO transposeRight:NO resultRows:1024 resultColumns:1024 interiorColumns:1024 alpha:1.0 beta:0.0];

// Encode the kernel
id<MTLCommandBuffer> commandBuffer = [commandQueue commandBuffer];
[matMul encodeToCommandBuffer:commandBuffer leftMatrix:matrixA rightMatrix:matrixB resultMatrix:matrixC];
[commandBuffer commit];
[commandBuffer waitUntilCompleted];

// Result in bufferC

That's it! The kernel automatically chooses the optimal algorithm for the device.

1.4 Example: Convolution for Deep Learning

MPS provides a full suite of neural network layers. Here's a simple convolution layer:

MPSCNNConvolutionDescriptor *convDesc = [MPSCNNConvolutionDescriptor cnnConvolutionDescriptorWithKernelWidth:3 kernelHeight:3 inputFeatureChannels:3 outputFeatureChannels:64];
convDesc.strideInPixelsX = 1; convDesc.strideInPixelsY = 1;

MPSCNNConvolution *conv = [[MPSCNNConvolution alloc] initWithDevice:device convolutionDescriptor:convDesc kernelWeights:weightsData biasTerms:biasData flags:MPSCNNConvolutionFlagsNone];

// Assume you have an MPSImage *srcImage
[conv encodeToCommandBuffer:commandBuffer sourceImage:srcImage destinationImage:dstImage];

MPS automatically handles weight layout and uses hardware‑accelerated matrix multiplication or Winograd algorithms.

1.5 MPSRayIntersector for Ray Tracing

MPSRayIntersector provides a fast, high‑level interface for ray tracing on Apple GPUs. It builds acceleration structures (BVHs) and performs ray‑primitive intersection tests.

Basic usage:

Create MPSRayIntersector with a device.
Build an acceleration structure from vertex and index data.
Encode intersection tests with rays.

MPSRayIntersector *intersector = [[MPSRayIntersector alloc] initWithDevice:device];

// Create a vertex buffer and index buffer
id<MTLBuffer> vertexBuffer = [device newBufferWithBytes:vertices length:vertexDataSize options:MTLResourceStorageModeManaged];
id<MTLBuffer> indexBuffer = [device newBufferWithBytes:indices length:indexDataSize options:MTLResourceStorageModeManaged];

// Build acceleration structure
MPSPolygonBuffer *polygonBuffer = [MPSPolygonBuffer polygonBufferWithDevice:device];
polygonBuffer.vertexBuffer = vertexBuffer;
polygonBuffer.indexBuffer = indexBuffer;
NSArray<MPSPolygonBuffer *> *buffers = @[polygonBuffer];
id<MTLBuffer> vertexCountBuffer = ...; // optional per‑polygon buffer

MPSTriangleAccelerationStructure *accel = [[MPSTriangleAccelerationStructure alloc] initWithDevice:device];
accel.polygonBuffers = buffers;
[accel rebuild];

// Prepare ray buffer (with origin, direction, etc.)
id<MTLBuffer> rayBuffer = [device newBufferWithLength:rayCount * sizeof(MPSRayOriginMaskDirection) options:MTLResourceStorageModeManaged];
MPSRayOriginMaskDirection *rays = (MPSRayOriginMaskDirection*)[rayBuffer contents];
// fill rays...

// Encode intersection
MPSIntersectionDataType intersectionDataType = MPSIntersectionDataTypeDistance;
id<MTLBuffer> intersectionBuffer = [device newBufferWithLength:rayCount * sizeof(MPSIntersectionDistance) options:MTLResourceStorageModeManaged];
[intersector encodeIntersectionToCommandBuffer:commandBuffer intersectionType:MPSIntersectionTypeNearest rayBuffer:rayBuffer rayBufferOffset:0 intersectionBuffer:intersectionBuffer intersectionBufferOffset:0 rayCount:rayCount accelerationStructure:accel];

You can then read the intersection distances from intersectionBuffer. MPSRayIntersector supports both triangle and sphere primitives, and can handle large scenes efficiently.

1.6 Performance Tips for MPS

Reuse MPS kernels and acceleration structures across frames.
Use MPSKernelOptions to tune for performance (e.g., MPSKernelOptionsAllowReducedPrecision for FP16).
For neural networks, consider using MPSNNGraph to fuse multiple layers into a single optimized graph.
Profile with Xcode's Metal System Trace to see how MPS kernels are scheduled.

Part 2: Metal Ray Tracing

Beyond MPSRayIntersector, Metal provides a more explicit ray tracing API that integrates with the Metal shading pipeline. This allows you to write custom ray generation, any‑hit, closest‑hit, and miss shaders.

2.1 Metal Ray Tracing Concepts

Acceleration Structure: A BVH built from geometry (primitives).
Intersection Function: A shader that determines if a ray hits a primitive. Metal provides a default triangle intersection, but you can write custom ones.
Ray Tracing Pipeline: Similar to a compute pipeline, but specialized for ray tracing.
Shader Stages: raygen, intersection, anyhit, closesthit, miss.

2.2 Setting Up a Basic Ray Tracer

Step 1: Create an acceleration structure
Similar to MPS, but using Metal's MTLAccelerationStructure APIs. You build a primitive acceleration structure descriptor, then encode the build into a command buffer.

Step 2: Create a ray tracing pipeline
You need a MTLRayTracingPipelineDescriptor that specifies the shader functions.

Step 3: Allocate resources

Ray buffer (origin, direction, etc.)
Intersection result buffer
Texture for output

Step 4: Encode ray tracing work
Use MTLRayTracingCommandEncoder to set the pipeline, bind resources, and dispatch rays.

2.3 Example: Simple Ray Generation Shader (MSL)

#include <metal_stdlib>
using namespace metal;

struct Ray {
    float3 origin;
    float3 direction;
};

struct Intersection {
    float distance;
    uint primitiveIndex;
    float2 barycentricCoord;
};

kernel void raygen(device Ray* rays [[buffer(0)]],
                   device Intersection* intersections [[buffer(1)]],
                   texture2d<float, access::write> output [[texture(0)]],
                   uint2 tid [[thread_position_in_grid]]) {
    // Generate ray for pixel (tid.x, tid.y)
    float2 uv = float2(tid) / float2(output.get_width(), output.get_height());
    rays[tid.y * output.get_width() + tid.x] = Ray{ float3(uv * 2.0f - 1.0f, 0.0f), float3(0,0,1) };

    // Wait for intersection results (simplified)
    Intersection hit = intersections[tid.y * output.get_width() + tid.x];
    float3 color = (hit.distance < 1e30) ? float3(1,0,0) : float3(0,0,0);
    output.write(float4(color, 1), tid);
}

You'll also need a closest‑hit shader that computes color based on material.

2.4 Building the Pipeline

On the host side, you compile the shaders and create a pipeline:

id<MTLLibrary> library = [device newDefaultLibrary];
id<MTLFunction> raygenFunc = [library newFunctionWithName:@"raygen"];
id<MTLFunction> closestHitFunc = [library newFunctionWithName:@"closesthit"];
id<MTLFunction> missFunc = [library newFunctionWithName:@"miss"];

MTLRayTracingPipelineDescriptor *desc = [MTLRayTracingPipelineDescriptor new];
desc.rayGeneratorFunction = raygenFunc;
desc.closestHitFunction = closestHitFunc;
desc.missFunction = missFunc;
desc.maxCallStackDepth = 1; // for recursion
desc.maxPayloadSize = sizeof(float3); // per‑ray data
desc.maxAttributeSize = sizeof(float2); // barycentrics

id<MTLRayTracingPipelineState> pipeline = [device newRayTracingPipelineStateWithDescriptor:desc error:nil];

2.5 Encoding and Dispatch

id<MTLRayTracingCommandEncoder> rtEncoder = [commandBuffer rayTracingCommandEncoder];
[rtEncoder setRayTracingPipelineState:pipeline];
[rtEncoder setRayBuffer:rayBuffer offset:0 atIndex:0];
[rtEncoder setIntersectionBuffer:intersectionBuffer offset:0 atIndex:1];
[rtEncoder setTexture:outputTexture atIndex:0];

// Dispatch rays (one per pixel)
[rtEncoder dispatchRays:MTLSizeMake(width, height, 1)];
[rtEncoder endEncoding];

This is a minimal example. Full ray tracing requires handling acceleration structures, intersection queries, and multiple hit groups.

2.6 Performance Considerations

Use MTLAccelerationStructure for static geometry; rebuild only when needed.
Batch rays in coherent groups for better coherence.
Use the smallest possible payload and attribute sizes.
Profile with Xcode to identify bottlenecks.

Part 3: Vulkan Subgroup Operations

Subgroup operations are a powerful feature in Vulkan (and other APIs) that allow efficient communication and computation among a set of invocations that execute together (a subgroup, roughly corresponding to a CUDA warp or AMD wavefront). They can dramatically improve performance for reductions, scans, and other data‑parallel patterns.

3.1 What are Subgroups?

A subgroup is a collection of invocations within a workgroup that can communicate and synchronize using specialized operations. The size of a subgroup is implementation‑defined but can be queried. On NVIDIA, it's 32; on AMD, 64; on Intel, maybe 8,16,32.

Subgroup operations include:

Vote: subgroupAll(), subgroupAny(), subgroupBallot().
Broadcast: subgroupBroadcast().
Shuffle: subgroupShuffle(), subgroupShuffleXor().
Arithmetic: subgroupAdd(), subgroupMin(), etc. (reductions).
Scan: subgroupInclusiveAdd(), subgroupExclusiveAdd().

3.2 Enabling Subgroup Support in Vulkan

You need to check for the VK_EXT_subgroup_size_control extension or at least the core subgroup support (since Vulkan 1.1, subgroups are core but optional). In your shader, declare:

#version 450
#extension GL_KHR_shader_subgroup_arithmetic : enable
#extension GL_KHR_shader_subgroup_ballot : enable
#extension GL_KHR_shader_subgroup_shuffle : enable

3.3 Example: Parallel Reduction with Subgroups

Suppose you want to sum an array within a workgroup. Without subgroups, you'd use shared memory and barriers. With subgroups, you can do a hierarchical reduction:

#version 450
#extension GL_KHR_shader_subgroup_arithmetic : enable

layout(local_size_x = 256) in;
layout(binding = 0) buffer In { float data[]; } inBuf;
layout(binding = 1) buffer Out { float data[]; } outBuf;

void main() {
    uint gid = gl_GlobalInvocationID.x;
    uint lid = gl_LocalInvocationID.x;
    float val = inBuf.data[gid];

    // First reduce within subgroup
    float subgroupSum = subgroupAdd(val);

    // The first invocation in each subgroup writes to shared memory
    shared float sharedSums[gl_WorkGroupSize.x / gl_SubgroupSize];
    if (subgroupElect()) {
        sharedSums[lid / gl_SubgroupSize] = subgroupSum;
    }
    barrier();

    // One subgroup reduces the shared sums
    if (gl_SubgroupID == 0) {
        float total = (lid < gl_NumSubgroups) ? sharedSums[lid] : 0;
        total = subgroupAdd(total);
        if (subgroupElect() && lid == 0) {
            outBuf.data[gl_WorkGroupID.x] = total;
        }
    }
}

This is much faster than a naive loop because it leverages hardware acceleration.

3.4 Vulkan Subgroup Query on Host

You can query subgroup properties:

VkPhysicalDeviceSubgroupProperties subgroupProps = {};
subgroupProps.sType = VK_STRUCTURE_TYPE_PHYSICAL_DEVICE_SUBGROUP_PROPERTIES;
VkPhysicalDeviceProperties2 props2 = {};
props2.sType = VK_STRUCTURE_TYPE_PHYSICAL_DEVICE_PROPERTIES_2;
props2.pNext = &subgroupProps;
vkGetPhysicalDeviceProperties2(physicalDevice, &props2);

uint32_t subgroupSize = subgroupProps.subgroupSize;
// Check supported operations in subgroupProps.supportedOperations

3.5 Best Practices

Use subgroups when you need intra‑workgroup communication with low latency.
Avoid branching that diverges within a subgroup (all invocations in a subgroup should take the same path).
Prefer subgroup operations over shared memory for small reductions/scans.
Fall back to shared memory if subgroup operations are not supported (query via gl_SubgroupSize == 0 or use #ifdef).

Part 4: Vulkan Ray Tracing

Vulkan ray tracing extensions (VK_KHR_ray_tracing_pipeline, VK_KHR_acceleration_structure, etc.) provide a cross‑platform API for hardware‑accelerated ray tracing. They are supported on NVIDIA RTX, AMD RDNA2+, and Intel Arc GPUs.

4.1 Core Concepts

Acceleration Structure: A BVH built from bottom‑level (geometry) and top‑level (instances) structures.
Shader Binding Table (SBT): A table of shader handles for ray generation, hit groups, and miss shaders.
Ray Tracing Pipeline: Combines shader stages (raygen, miss, closesthit, anyhit, intersection).
Ray Tracing Commands: vkCmdTraceRaysKHR to launch rays.

4.2 Building Acceleration Structures

First, create buffers with geometry data. Then use vkCmdBuildAccelerationStructuresKHR to build bottom‑level and top‑level AS.

Example: Build a bottom‑level AS from a triangle mesh:

VkAccelerationStructureGeometryKHR geometry = {};
geometry.sType = VK_STRUCTURE_TYPE_ACCELERATION_STRUCTURE_GEOMETRY_KHR;
geometry.geometryType = VK_GEOMETRY_TYPE_TRIANGLES_KHR;
geometry.flags = VK_GEOMETRY_OPAQUE_BIT_KHR;
geometry.geometry.triangles.sType = VK_STRUCTURE_TYPE_ACCELERATION_STRUCTURE_GEOMETRY_TRIANGLES_DATA_KHR;
geometry.geometry.triangles.vertexFormat = VK_FORMAT_R32G32B32_SFLOAT;
geometry.geometry.triangles.vertexData.deviceAddress = vertexBufferAddress;
geometry.geometry.triangles.vertexStride = sizeof(Vertex);
geometry.geometry.triangles.maxVertex = vertexCount - 1;
geometry.geometry.triangles.indexType = VK_INDEX_TYPE_UINT32;
geometry.geometry.triangles.indexData.deviceAddress = indexBufferAddress;

VkAccelerationStructureBuildRangeInfoKHR range = {};
range.primitiveCount = indexCount / 3;
range.primitiveOffset = 0;
range.firstVertex = 0;
range.transformOffset = 0;

// Create and build bottom AS...

You'll need to query the required memory size, allocate, and then build.

4.3 Shader Example (GLSL)

Ray generation shader:

#version 460
#extension GL_EXT_ray_tracing : require

layout(binding = 0) uniform accelerationStructureEXT topLevelAS;
layout(binding = 1, rgba8) uniform image2D image;
layout(location = 0) rayPayloadEXT vec3 hitColor;

void main() {
    ivec2 pixel = ivec2(gl_LaunchIDEXT.xy);
    vec2 dim = vec2(gl_LaunchSizeEXT.xy);
    vec2 d = (vec2(pixel) + vec2(0.5)) / dim - 0.5;
    float aspect = dim.x / dim.y;
    vec3 origin = vec3(0,0,-2);
    vec3 direction = normalize(vec3(d.x * aspect, -d.y, 1));

    traceRayEXT(topLevelAS, gl_RayFlagsOpaqueEXT, 0xff, 0, 0, 0, origin, 0.001, direction, 1000.0, 0);
    imageStore(image, pixel, vec4(hitColor, 1));
}

Closest hit shader:

#version 460
#extension GL_EXT_ray_tracing : require

layout(location = 0) rayPayloadInEXT vec3 hitColor;
hitAttributeEXT vec2 attribs;

void main() {
    hitColor = vec3(1,0,0); // red
}

Miss shader:

#version 460
#extension GL_EXT_ray_tracing : require

layout(location = 0) rayPayloadInEXT vec3 hitColor;

void main() {
    hitColor = vec3(0,0,0); // black
}

4.4 Creating a Ray Tracing Pipeline

You need to create shader modules for each stage, then a ray tracing pipeline with VkRayTracingPipelineCreateInfoKHR. This includes the shader groups (raygen, miss, hit groups).

4.5 Shader Binding Table (SBT)

The SBT is a buffer that holds shader handles and any local data for each shader group. You must fill it according to the pipeline layout. It's complex but crucial.

4.6 Dispatching Rays

After recording the SBT buffer, use:

vkCmdTraceRaysKHR(commandBuffer, &raygenSBT, &missSBT, &hitSBT, &callableSBT, width, height, 1);

4.7 Performance Tips

Use opaque geometry and avoid any‑hit shaders when possible.
Keep ray payload small.
Group rays coherently (e.g., sort by direction).
Use inline ray tracing (VK_KHR_ray_query) for simpler tasks without full pipeline overhead.

Part 5: Vulkan Video

Vulkan Video extensions (VK_KHR_video_queue, VK_KHR_video_decode_queue, VK_KHR_video_encode_queue) provide hardware‑accelerated video encoding and decoding. This is a relatively new area but powerful for applications needing real‑time video processing.

5.1 Overview

Vulkan Video exposes video codec capabilities as Vulkan queues and operations. You can decode H.264, H.265, etc., directly to Vulkan images, and encode from images. It integrates seamlessly with the rest of Vulkan for post‑processing and rendering.

5.2 Setting Up Video Queues

When creating a logical device, you need to request a video decode/encode queue:

float priority = 1.0f;
VkDeviceQueueCreateInfo queueCreateInfo = {};
queueCreateInfo.sType = VK_STRUCTURE_TYPE_DEVICE_QUEUE_CREATE_INFO;
queueCreateInfo.queueFamilyIndex = videoQueueFamilyIndex;
queueCreateInfo.queueCount = 1;
queueCreateInfo.pQueuePriorities = &priority;

You must first query queue families that support VK_QUEUE_VIDEO_DECODE_BIT_KHR or VK_QUEUE_VIDEO_ENCODE_BIT_KHR.

5.3 Decoding a Video Frame

Create a video session (VkVideoSessionKHR) with codec parameters.
Create video picture resources (reference and output images).
Submit decode operations via vkCmdDecodeVideoKHR.

Example decode command:

VkVideoDecodeInfoKHR decodeInfo = {};
decodeInfo.sType = VK_STRUCTURE_TYPE_VIDEO_DECODE_INFO_KHR;
decodeInfo.srcBuffer = bitstreamBuffer;
decodeInfo.srcBufferOffset = bitstreamOffset;
decodeInfo.srcBufferRange = bitstreamSize;
decodeInfo.dstPictureResource = pictureResource; // reference to image view
decodeInfo.pSetupReferenceSlots = ...; // for DPB

vkCmdDecodeVideoKHR(commandBuffer, &decodeInfo);

You'll need to manage the decoded picture buffer (DPB) and reference frames.

5.4 Encoding

Similarly, you can encode frames by providing raw images and receiving a bitstream.

5.5 Limitations and Considerations

Vulkan Video is still evolving; not all vendors support all codecs.
You need to handle bitstream parsing yourself (the API expects raw NAL units).
It's a low‑level interface; for simpler use, consider libraries like FFmpeg with Vulkan acceleration.

Part 6: Cross‑Platform Considerations

While Metal and Vulkan are powerful, they are platform‑specific. For applications that need to run on both Apple and non‑Apple hardware, you have several options:

6.1 SYCL

SYCL (as covered in previous tutorials) is a high‑level C++ abstraction that can target multiple backends including OpenCL, Level Zero, and CUDA. With SYCL 2020, you can write a single source that compiles to various devices. However, SYCL does not yet directly target Metal or Vulkan compute, but it can target OpenCL which runs on many platforms (including macOS via Intel GPU or with limitations). For Apple Silicon, SYCL is not a direct fit; you'd need to use OpenCL (which is deprecated on macOS) or a SYCL implementation that supports Metal (none currently).

6.2 OpenCL

OpenCL is widely supported (including on Intel GPUs, AMD, NVIDIA), but Apple has deprecated OpenCL on macOS in favor of Metal. For cross‑platform, you might write OpenCL kernels and use them on non‑Apple platforms, and rewrite for Metal on Apple. This is double the work.

6.3 Vulkan Everywhere?

Vulkan runs on Windows, Linux, Android, and macOS via MoltenVK. So you could write your compute kernels in Vulkan and use MoltenVK to run on macOS. MoltenVK translates Vulkan to Metal, and for compute, it works quite well. This gives you a single codebase for all platforms. However, you lose some Metal‑specific optimizations, but the portability is valuable.

For ray tracing, Vulkan ray tracing extensions are not yet supported through MoltenVK (Metal ray tracing is different). So if you need ray tracing on Apple, you'd need separate Metal code.

6.4 Recommendations

If you target only Apple: Use Metal and MPS for maximum performance.
If you target multiple platforms and need compute only: Use Vulkan with MoltenVK for macOS, or SYCL with appropriate backends.
If you need ray tracing across platforms: Write separate backends for Metal and Vulkan, as the APIs differ significantly.
If you need video encode/decode: Vulkan Video is cross‑vendor but not yet on Apple; use platform‑specific APIs (VideoToolbox on macOS, MediaFoundation on Windows, VAAPI on Linux).

Conclusion

You've now explored the advanced corners of Metal and Vulkan: the high‑level performance primitives of MPS, the excitement of ray tracing on both platforms, the efficiency of subgroup operations, and the emerging field of Vulkan Video. Armed with this knowledge, you can build applications that push the boundaries of what's possible on modern GPUs.

Remember, the key to mastery is practice. Implement a small path tracer in Metal, try a subgroup‑accelerated reduction in Vulkan, decode a video stream with Vulkan Video, and see how MPS can accelerate your machine learning models. The possibilities are endless.

If you have questions or want to share your own projects, drop a comment below. Happy coding!

Advanced GPU Optimization: Metal & Vulkan Compute from zero to hero

Javad — Sat, 14 Feb 2026 16:38:13 +0000

Welcome to the exciting world of GPU compute programming! In this comprehensive guide, we'll explore two powerful, low‑level APIs: Metal for Apple devices (especially Apple Silicon) and Vulkan Compute for cross‑platform high‑performance computing. By the end, you'll be able to write efficient compute kernels, manage GPU resources, and optimize for real‑world tasks like matrix multiplication and image processing.

Whether you're a macOS/iOS developer wanting to leverage the Neural Engine and GPU on Apple Silicon, or a cross‑platform developer aiming for maximum performance on Windows, Linux, or Android, this tutorial has you covered.

Introduction: Why Metal and Vulkan Compute?

Metal

Metal is Apple's graphics and compute API, available on all modern Apple devices (iPhone, iPad, Mac). With the advent of Apple Silicon (M1, M2, M3 families), the GPU is tightly integrated with the CPU, offering unified memory and exceptional performance. Metal provides a low‑overhead, explicit API that gives you fine‑grained control over the GPU.

Vulkan Compute

Vulkan is a cross‑platform graphics and compute API from the Khronos Group. Its compute capabilities are exposed through compute pipelines and shaders written in GLSL or HLSL (compiled to SPIR‑V). Vulkan is designed for high‑performance, low‑latency applications and runs on Windows, Linux, Android, macOS (via MoltenVK), and more.

Both Metal and Vulkan follow a similar explicit, command‑buffer‑based model. Understanding one makes learning the other easier. Let's dive in!

Part 1: Metal on Apple Silicon

1.1 Setting Up Metal

To develop with Metal, you'll need:

A Mac with Xcode installed (Xcode includes the Metal compiler, framework, and debugging tools).
An Apple Silicon Mac (M1 or later) for native testing; Intel Macs also work but performance characteristics differ.
For iOS/tvOS, you'll need an Apple Developer account (free for device testing).

Create a new Xcode project (e.g., a command‑line tool for macOS) and add the Metal framework. In your source file, import Metal:

#import <Metal/Metal.h>
// or in Swift: import Metal

1.2 Metal Shading Language (MSL)

Metal uses its own C++‑based shading language, MSL. For compute, you write kernels (functions with kernel keyword). Basic MSL features:

Data types: float, int, half, vectors (float4, int2), matrices (float4x4).
Address spaces: device (global memory), threadgroup (shared memory within a workgroup), constant (read‑only), thread (private).
Built‑in variables:
- uint3 [[thread_position_in_grid]] – global thread ID.
- uint3 [[thread_position_in_threadgroup]] – local thread ID within threadgroup.
- uint3 [[threadgroup_position_in_grid]] – threadgroup ID.
- uint3 [[threads_per_threadgroup]] – size of threadgroup.

Example kernel that adds two vectors:

#include <metal_stdlib>
using namespace metal;

kernel void vecAdd(device const float* a [[buffer(0)]],
                   device const float* b [[buffer(1)]],
                   device float* c [[buffer(2)]],
                   uint id [[thread_position_in_grid]]) {
    c[id] = a[id] + b[id];
}

Note the use of [[buffer(n)]] to bind resources.

1.3 Metal Compute Pipeline

The host side (Objective‑C or Swift) sets up the pipeline:

MTLDevice: The abstract representation of the GPU.
MTLCommandQueue: Queue for submitting command buffers.
MTLComputePipelineState: Compiled compute pipeline (kernel + device).
MTLCommandBuffer: Container for commands.
MTLComputeCommandEncoder: Encodes compute commands.
Dispatch: Specify threadgroup dimensions.

1.4 Your First Metal Compute Program: Vector Addition

Let's walk through the host code in Objective‑C. We'll assume you have a Metal device and have compiled the kernel source.

Step 1: Get the default device

id<MTLDevice> device = MTLCreateSystemDefaultDevice();

Step 2: Create a command queue

id<MTLCommandQueue> commandQueue = [device newCommandQueue];

Step 3: Load the kernel

// Assume you have a .metal file in your project named "kernel.metal"
id<MTLLibrary> library = [device newDefaultLibrary];
id<MTLFunction> kernelFunction = [library newFunctionWithName:@"vecAdd"];

Step 4: Create the pipeline state

NSError *error = nil;
id<MTLComputePipelineState> pipelineState = [device newComputePipelineStateWithFunction:kernelFunction error:&error];

Step 5: Create buffers

NSUInteger dataSize = VEC_SIZE * sizeof(float);

id<MTLBuffer> bufferA = [device newBufferWithBytes:h_a length:dataSize options:MTLResourceStorageModeShared];
id<MTLBuffer> bufferB = [device newBufferWithBytes:h_b length:dataSize options:MTLResourceStorageModeShared];
id<MTLBuffer> bufferC = [device newBufferWithLength:dataSize options:MTLResourceStorageModeShared];

MTLResourceStorageModeShared means the buffer is accessible by both CPU and GPU (unified memory on Apple Silicon).

Step 6: Create command buffer and encoder

id<MTLCommandBuffer> commandBuffer = [commandQueue commandBuffer];
id<MTLComputeCommandEncoder> encoder = [commandBuffer computeCommandEncoder];
[encoder setComputePipelineState:pipelineState];
[encoder setBuffer:bufferA offset:0 atIndex:0];
[encoder setBuffer:bufferB offset:0 atIndex:1];
[encoder setBuffer:bufferC offset:0 atIndex:2];

Step 7: Dispatch threads

MTLSize threadgroupSize = MTLSizeMake(256, 1, 1); // threads per threadgroup
MTLSize threadgroupCount = MTLSizeMake((VEC_SIZE + 255) / 256, 1, 1);
[encoder dispatchThreadgroups:threadgroupCount threadsPerThreadgroup:threadgroupSize];

Step 8: End encoding and commit

[encoder endEncoding];
[commandBuffer commit];
[commandBuffer waitUntilCompleted]; // or use completion handler

Step 9: Read results

float *c_ptr = (float*)[bufferC contents];
// c_ptr now contains the result

Step 10: Verify and clean up (ARC will handle objects).

That's it! You've just run your first Metal compute kernel.

1.5 Memory Management in Metal

Metal buffers are created with a storage mode:

MTLResourceStorageModeShared: CPU and GPU access same memory (ideal for Apple Silicon).
MTLResourceStorageModePrivate: Only GPU accessible; used for intermediate buffers.
MTLResourceStorageModeManaged: For discrete GPUs (Intel Macs); requires explicit synchronization.

On Apple Silicon, use Shared for simplicity and performance.

You can also create textures (id<MTLTexture>) for image data. Textures are useful for compute kernels that need spatial locality.

1.6 Performance Considerations

Threadgroup memory (shared memory)

Metal provides threadgroup memory (like CUDA shared memory) for fast data exchange within a workgroup. Declare it in the kernel:

kernel void myKernel(device float* in [[buffer(0)]],
                     device float* out [[buffer(1)]],
                     threadgroup float* shared [[threadgroup(0)]],
                     uint tid [[thread_position_in_threadgroup]],
                     uint gid [[thread_position_in_grid]]) {
    // use shared[tid] ...
}

On the host, you specify the size of threadgroup memory when dispatching.

Coalescing and bank conflicts

Threads in a threadgroup access memory. Ensure that within a warp (32 threads on Apple GPUs), accesses are coalesced. Avoid bank conflicts in threadgroup memory by accessing different banks.

Occupancy and threadgroup size

Use pipelineState.threadExecutionWidth to get the optimal threadgroup size multiple. Typical sizes are multiples of 32.

1.7 Advanced Metal Compute

Using textures

Compute kernels can read/write textures. Useful for image processing.

kernel void sobel(texture2d<float, access::read> inTexture [[texture(0)]],
                  texture2d<float, access::write> outTexture [[texture(1)]],
                  uint2 gid [[thread_position_in_grid]]) {
    float4 color = inTexture.read(gid);
    // compute...
    outTexture.write(result, gid);
}

Multiple dispatches and synchronization

You can encode multiple dispatches in one command buffer, and even use MTLComputeCommandEncoder's setThreadgroupMemoryLength to adjust per‑dispatch.

Profiling

Use Xcode's GPU Debugger and Metal System Trace to analyze performance, view occupancy, and find bottlenecks.

1.8 Example: Tiled Matrix Multiplication in Metal

Let's implement a tiled matrix multiplication kernel using threadgroup memory.

Kernel (MSL):

kernel void matMulTiled(device const float* A [[buffer(0)]],
                        device const float* B [[buffer(1)]],
                        device float* C [[buffer(2)]],
                        constant uint& N [[buffer(3)]],
                        threadgroup float* As [[threadgroup(0)]],
                        threadgroup float* Bs [[threadgroup(1)]],
                        uint2 tgid [[threadgroup_position_in_grid]],
                        uint2 lid [[thread_position_in_threadgroup]],
                        uint2 dims [[threads_per_threadgroup]]) {
    uint tileSize = dims.x; // assuming square threadgroup
    uint row = tgid.y * tileSize + lid.y;
    uint col = tgid.x * tileSize + lid.x;

    float sum = 0.0f;
    for (uint tile = 0; tile < N / tileSize; ++tile) {
        // Load tile of A and B
        As[lid.y * tileSize + lid.x] = A[row * N + (tile * tileSize + lid.x)];
        Bs[lid.y * tileSize + lid.x] = B[(tile * tileSize + lid.y) * N + col];
        threadgroup_barrier(mem_flags::mem_threadgroup);

        // Compute partial product
        for (uint k = 0; k < tileSize; ++k) {
            sum += As[lid.y * tileSize + k] * Bs[k * tileSize + lid.x];
        }
        threadgroup_barrier(mem_flags::mem_threadgroup);
    }
    C[row * N + col] = sum;
}

Host setup:

Buffer A, B, C.
Set N as a constant buffer.
Dispatch with threadgroup size (16,16) and grid size (N/16, N/16).

This kernel leverages threadgroup memory for data reuse, significantly improving performance.

Part 2: Vulkan Compute

2.1 Setting Up Vulkan

To develop Vulkan compute applications, you'll need:

Vulkan SDK from LunarG. Install and set VULKAN_SDK environment variable.
A Vulkan‑capable GPU and drivers (most modern GPUs from NVIDIA, AMD, Intel).
For macOS, you can use MoltenVK (included in the SDK) which translates Vulkan to Metal. Performance is good for compute.

Create a C++ project and include vulkan/vulkan.h. Link against the Vulkan loader (vulkan-1 on Windows, vulkan on Linux).

2.2 Vulkan Compute Pipeline

The Vulkan compute pipeline involves several steps:

Create a VkInstance (with validation layers for debugging).
Pick a physical device (VkPhysicalDevice) that supports compute.
Create a logical device (VkDevice) with a compute queue.
Create a command pool and allocate command buffers.
Create shader module from SPIR‑V.
Create descriptor set layout and pipeline layout.
Create compute pipeline.
Create buffers and allocate device memory.
Update descriptor sets to bind buffers.
Record command buffer: bind pipeline, bind descriptor sets, dispatch.
Submit to queue and synchronize.

Let's go through each step with code.

2.3 GLSL for Compute Shaders

We'll write compute shaders in GLSL (or HLSL compiled to SPIR‑V). A simple vector addition shader:

#version 450
layout(local_size_x = 256, local_size_y = 1, local_size_z = 1) in;

layout(binding = 0) buffer InputA { float data[]; } a;
layout(binding = 1) buffer InputB { float data[]; } b;
layout(binding = 2) buffer OutputC { float data[]; } c;

void main() {
    uint id = gl_GlobalInvocationID.x;
    c.data[id] = a.data[id] + b.data[id];
}

Compile this to SPIR‑V using glslangValidator or glslc. For example:

glslc -c -o vecAdd.spv vecAdd.comp

2.4 Memory and Resources

Vulkan uses VkBuffer objects backed by VkDeviceMemory. You must allocate memory, bind it, and map it for host access.

Descriptor sets are used to bind resources (buffers, images) to shaders. You need:

VkDescriptorSetLayout describing the bindings.
VkDescriptorPool to allocate descriptor sets.
VkDescriptorSet containing actual buffer references.

Push constants are small amounts of data that can be set directly in command buffers (fast).

2.5 Your First Vulkan Compute Program: Vector Addition

Let's outline the host code in C++. Error checking omitted for brevity.

Step 1: Create instance and pick device (standard Vulkan initialization). We'll assume you have a VkDevice and a compute queue.

Step 2: Create shader module

VkShaderModuleCreateInfo createInfo = {};
createInfo.sType = VK_STRUCTURE_TYPE_SHADER_MODULE_CREATE_INFO;
createInfo.codeSize = codeSize;
createInfo.pCode = (uint32_t*)spvCode;
vkCreateShaderModule(device, &createInfo, nullptr, &shaderModule);

Step 3: Create descriptor set layout (bindings for three buffers)

VkDescriptorSetLayoutBinding bindings[3] = {};
bindings[0].binding = 0; bindings[0].descriptorType = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER; bindings[0].descriptorCount = 1; bindings[0].stageFlags = VK_SHADER_STAGE_COMPUTE_BIT;
// similarly for 1 and 2
VkDescriptorSetLayoutCreateInfo layoutInfo = {};
layoutInfo.sType = VK_STRUCTURE_TYPE_DESCRIPTOR_SET_LAYOUT_CREATE_INFO;
layoutInfo.bindingCount = 3; layoutInfo.pBindings = bindings;
vkCreateDescriptorSetLayout(device, &layoutInfo, nullptr, &descriptorSetLayout);

Step 4: Create pipeline layout

VkPipelineLayoutCreateInfo pipelineLayoutInfo = {};
pipelineLayoutInfo.sType = VK_STRUCTURE_TYPE_PIPELINE_LAYOUT_CREATE_INFO;
pipelineLayoutInfo.setLayoutCount = 1; pipelineLayoutInfo.pSetLayouts = &descriptorSetLayout;
vkCreatePipelineLayout(device, &pipelineLayoutInfo, nullptr, &pipelineLayout);

Step 5: Create compute pipeline

VkComputePipelineCreateInfo pipelineInfo = {};
pipelineInfo.sType = VK_STRUCTURE_TYPE_COMPUTE_PIPELINE_CREATE_INFO;
pipelineInfo.stage.sType = VK_STRUCTURE_TYPE_PIPELINE_SHADER_STAGE_CREATE_INFO;
pipelineInfo.stage.stage = VK_SHADER_STAGE_COMPUTE_BIT;
pipelineInfo.stage.module = shaderModule;
pipelineInfo.stage.pName = "main";
pipelineInfo.layout = pipelineLayout;
vkCreateComputePipelines(device, VK_NULL_HANDLE, 1, &pipelineInfo, nullptr, &pipeline);

Step 6: Create buffers and allocate memory

VkBufferCreateInfo bufferInfo = {};
bufferInfo.sType = VK_STRUCTURE_TYPE_BUFFER_CREATE_INFO;
bufferInfo.size = dataSize;
bufferInfo.usage = VK_BUFFER_USAGE_STORAGE_BUFFER_BIT | VK_BUFFER_USAGE_TRANSFER_DST_BIT; // for copying from host
bufferInfo.sharingMode = VK_SHARING_MODE_EXCLUSIVE;
vkCreateBuffer(device, &bufferInfo, nullptr, &bufferA);
// Similarly for B, C

// Get memory requirements, allocate memory, bind buffers
VkMemoryRequirements memReqs;
vkGetBufferMemoryRequirements(device, bufferA, &memReqs);
VkMemoryAllocateInfo allocInfo = {};
allocInfo.sType = VK_STRUCTURE_TYPE_MEMORY_ALLOCATE_INFO;
allocInfo.allocationSize = memReqs.size;
allocInfo.memoryTypeIndex = findMemoryType(memReqs.memoryTypeBits, VK_MEMORY_PROPERTY_HOST_VISIBLE_BIT | VK_MEMORY_PROPERTY_HOST_COHERENT_BIT);
vkAllocateMemory(device, &allocInfo, nullptr, &memoryA);
vkBindBufferMemory(device, bufferA, memoryA, 0);

Step 7: Map memory and copy initial data

void* mapped;
vkMapMemory(device, memoryA, 0, dataSize, 0, &mapped);
memcpy(mapped, h_a, dataSize);
vkUnmapMemory(device, memoryA);
// Similarly for B

Step 8: Create descriptor pool and descriptor set

VkDescriptorPoolSize poolSize = {};
poolSize.type = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER;
poolSize.descriptorCount = 3;
VkDescriptorPoolCreateInfo poolInfo = {};
poolInfo.sType = VK_STRUCTURE_TYPE_DESCRIPTOR_POOL_CREATE_INFO;
poolInfo.maxSets = 1;
poolInfo.poolSizeCount = 1; poolInfo.pPoolSizes = &poolSize;
vkCreateDescriptorPool(device, &poolInfo, nullptr, &descriptorPool);

VkDescriptorSetAllocateInfo allocInfo = {};
allocInfo.sType = VK_STRUCTURE_TYPE_DESCRIPTOR_SET_ALLOCATE_INFO;
allocInfo.descriptorPool = descriptorPool;
allocInfo.descriptorSetCount = 1; allocInfo.pSetLayouts = &descriptorSetLayout;
vkAllocateDescriptorSets(device, &allocInfo, &descriptorSet);

Step 9: Update descriptor set with buffer info

VkDescriptorBufferInfo bufferInfoA = { bufferA, 0, dataSize };
VkWriteDescriptorSet writeA = {};
writeA.sType = VK_STRUCTURE_TYPE_WRITE_DESCRIPTOR_SET;
writeA.dstSet = descriptorSet;
writeA.dstBinding = 0;
writeA.descriptorCount = 1;
writeA.descriptorType = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER;
writeA.pBufferInfo = &bufferInfoA;
// Similarly for B, C (bindings 1,2)
vkUpdateDescriptorSets(device, 3, writes, 0, nullptr);

Step 10: Create command pool and command buffer

VkCommandPoolCreateInfo poolCmdInfo = {};
poolCmdInfo.sType = VK_STRUCTURE_TYPE_COMMAND_POOL_CREATE_INFO;
poolCmdInfo.queueFamilyIndex = computeQueueFamilyIndex;
vkCreateCommandPool(device, &poolCmdInfo, nullptr, &commandPool);

VkCommandBufferAllocateInfo cmdAllocInfo = {};
cmdAllocInfo.sType = VK_STRUCTURE_TYPE_COMMAND_BUFFER_ALLOCATE_INFO;
cmdAllocInfo.commandPool = commandPool;
cmdAllocInfo.level = VK_COMMAND_BUFFER_LEVEL_PRIMARY;
cmdAllocInfo.commandBufferCount = 1;
vkAllocateCommandBuffers(device, &cmdAllocInfo, &commandBuffer);

Step 11: Record commands

VkCommandBufferBeginInfo beginInfo = {};
beginInfo.sType = VK_STRUCTURE_TYPE_COMMAND_BUFFER_BEGIN_INFO;
vkBeginCommandBuffer(commandBuffer, &beginInfo);

vkCmdBindPipeline(commandBuffer, VK_PIPELINE_BIND_POINT_COMPUTE, pipeline);
vkCmdBindDescriptorSets(commandBuffer, VK_PIPELINE_BIND_POINT_COMPUTE, pipelineLayout, 0, 1, &descriptorSet, 0, nullptr);

uint32_t groupCountX = (VEC_SIZE + 255) / 256;
vkCmdDispatch(commandBuffer, groupCountX, 1, 1);

vkEndCommandBuffer(commandBuffer);

Step 12: Submit and wait

VkSubmitInfo submitInfo = {};
submitInfo.sType = VK_STRUCTURE_TYPE_SUBMIT_INFO;
submitInfo.commandBufferCount = 1; submitInfo.pCommandBuffers = &commandBuffer;
vkQueueSubmit(computeQueue, 1, &submitInfo, VK_NULL_HANDLE);
vkQueueWaitIdle(computeQueue);

Step 13: Read back results (map memory of buffer C and copy).

That's the complete flow. It's verbose, but gives you full control.

2.6 Synchronization

Vulkan uses fences (GPU‑CPU sync) and semaphores (GPU‑GPU sync). In compute, you often just wait for the queue idle. For finer control, use vkWaitForFences after submitting with a fence.

Memory dependencies are handled by pipeline barriers. For example, if you write to a buffer then read it in a later dispatch, you need a barrier to ensure visibility.

2.7 Performance Optimization

Workgroup size

Choose a workgroup size that is a multiple of the subgroup size (typically 32 or 64) and that maximizes occupancy. Use VkPhysicalDeviceProperties::maxComputeWorkGroupSize.

Shared memory

Declare shared memory in GLSL with shared float tile[16][16]; and use barrier() for synchronization. Similar to Metal threadgroup memory.

Subgroup operations

Vulkan exposes subgroup operations via the VK_EXT_subgroup extension. These allow efficient reductions, scans, and shuffles within a subgroup (like warp intrinsics).

Avoiding bank conflicts – same principle as Metal.

2.8 Example: Tiled Matrix Multiplication in Vulkan

Shader (GLSL):

#version 450
layout(local_size_x = 16, local_size_y = 16) in;

layout(binding = 0) buffer InputA { float data[]; } A;
layout(binding = 1) buffer InputB { float data[]; } B;
layout(binding = 2) buffer OutputC { float data[]; } C;
layout(push_constant) uniform PushConsts { uint N; } pc;

shared float As[16][16];
shared float Bs[16][16];

void main() {
    uint row = gl_GlobalInvocationID.y;
    uint col = gl_GlobalInvocationID.x;
    uint localRow = gl_LocalInvocationID.y;
    uint localCol = gl_LocalInvocationID.x;

    float sum = 0.0;
    for (uint tile = 0; tile < pc.N / 16; ++tile) {
        As[localRow][localCol] = A.data[row * pc.N + (tile * 16 + localCol)];
        Bs[localRow][localCol] = B.data[(tile * 16 + localRow) * pc.N + col];
        barrier();

        for (uint k = 0; k < 16; ++k) {
            sum += As[localRow][k] * Bs[k][localCol];
        }
        barrier();
    }
    C.data[row * pc.N + col] = sum;
}

Host changes:

Use push constant to pass N.
Dispatch with grid size (N/16, N/16) and local size (16,16).
Ensure descriptor set binds buffers A, B, C.

Part 3: Comparative Analysis and Cross‑Platform Strategies

3.1 Metal vs Vulkan Compute: Similarities and Differences

Aspect	Metal	Vulkan Compute
API Style	Objective‑C / Swift (objects, references)	C (handles, structures)
Shader Language	MSL (C++‑based)	GLSL / HLSL compiled to SPIR‑V
Resource Binding	[[buffer(n)]] and [[texture(n)]]	Descriptor sets with bindings
Memory Model	Unified on Apple Silicon; discrete otherwise	Explicit memory management
Workgroup/Threadgroup	`threadgroup` memory, `threadgroup_barrier`	`shared` memory, `barrier()`
Error Handling	NSError, exceptions (Swift)	Return codes (VkResult)
Debugging/Profiling	Xcode GPU Debugger, Metal System Trace	RenderDoc, Vulkan Validation Layers, GPU vendors' tools

Both are explicit, low‑level, and provide fine‑grained control. Metal is slightly simpler to use, especially on Apple platforms, while Vulkan offers cross‑platform portability.

3.2 Writing Portable Compute Code

If you need to support both Metal and Vulkan, consider:

Writing compute kernels in a common language (e.g., using a library like OpenCL C or SYCL that compiles to both, but that's beyond this guide).
Maintaining two separate codebases with shared algorithms.
Using conditional compilation and abstraction layers (like MoltenVK to run Vulkan on Metal, though compute performance is good).

For compute kernels, the core logic is often identical; only the syntax for resource binding and built‑in variables differs. You can use macros to abstract differences.

Example kernel snippet (pseudo‑code):

#ifdef METAL
#define GLOBAL_ID [[thread_position_in_grid]]
#define SHARED threadgroup
#define BARRIER threadgroup_barrier(mem_flags::mem_threadgroup)
#else
#define GLOBAL_ID gl_GlobalInvocationID
#define SHARED shared
#define BARRIER barrier()
#endif

Part 4: Advanced Topics

4.1 Indirect Dispatch

Both APIs support indirect dispatch where the workgroup counts are read from a buffer, allowing GPU‑generated workloads.

Metal: [encoder dispatchThreadgroupsWithIndirectBuffer:indirectBuffer indirectBufferOffset:0 threadsPerThreadgroup:threadgroupSize];

Vulkan: vkCmdDispatchIndirect with a buffer containing the dispatch dimensions.

4.2 GPU Timestamps and Queries

Metal: Use MTLCounterSampleBuffer (on supported devices) to sample GPU timestamps.

Vulkan: Use query pools (VK_QUERY_TYPE_TIMESTAMP) to record timestamps in command buffers.

4.3 Interop with Graphics Pipelines

Both APIs allow mixing compute and graphics in the same command buffer. For example, you can use a compute shader to generate data for a vertex buffer, then render it.

4.4 Multi‑GPU (Vulkan) and Multi‑Device (Metal)

Vulkan supports multi‑GPU via logical devices with multiple physical devices, but it's complex. Metal has MTLDevice per GPU; you can create separate queues and manage workloads manually.

4.5 Real‑world Example: Sobel Edge Detection

Let's outline a Sobel kernel in both APIs.

Metal Sobel (texture version):

kernel void sobel(texture2d<float, access::read> inTex [[texture(0)]],
                  texture2d<float, access::write> outTex [[texture(1)]],
                  uint2 gid [[thread_position_in_grid]]) {
    if (gid.x < outTex.get_width() && gid.y < outTex.get_height()) {
        float sumX = 0, sumY = 0;
        for (int dy = -1; dy <= 1; dy++) {
            for (int dx = -1; dx <= 1; dx++) {
                uint2 coord = uint2(gid.x + dx, gid.y + dy);
                float val = inTex.read(coord).r;
                sumX += val * sobelX[dy+1][dx+1];
                sumY += val * sobelY[dy+1][dx+1];
            }
        }
        float mag = min(1.0f, hypot(sumX, sumY));
        outTex.write(float4(mag, mag, mag, 1), gid);
    }
}

Vulkan Sobel (buffer version, assuming image stored as float array):

#version 450
layout(local_size_x = 16, local_size_y = 16) in;
layout(binding = 0) buffer Input { float pixels[]; } inBuf;
layout(binding = 1) buffer Output { float pixels[]; } outBuf;
layout(push_constant) uniform Dims { uint width; uint height; } dims;

void main() {
    uint x = gl_GlobalInvocationID.x;
    uint y = gl_GlobalInvocationID.y;
    if (x >= dims.width || y >= dims.height) return;
    float sumX = 0, sumY = 0;
    for (int dy = -1; dy <= 1; dy++) {
        for (int dx = -1; dx <= 1; dx++) {
            int nx = int(x) + dx, ny = int(y) + dy;
            if (nx < 0 || nx >= dims.width || ny < 0 || ny >= dims.height) continue;
            float val = inBuf.pixels[ny * dims.width + nx];
            sumX += val * sobelX[dy+1][dx+1];
            sumY += val * sobelY[dy+1][dx+1];
        }
    }
    float mag = min(1.0, length(vec2(sumX, sumY)));
    outBuf.pixels[y * dims.width + x] = mag;
}

Conclusion

Congratulations! You've journeyed through the depths of Metal compute on Apple Silicon and Vulkan compute. You now know how to set up the environment, write kernels, manage memory, optimize performance, and implement complex patterns like matrix multiplication and edge detection.

Next Steps

Metal: Explore Metal Performance Shaders (MPS) for highly optimized primitives (convolution, matrix multiplication, etc.). Look into Metal Ray Tracing for advanced graphics.
Vulkan: Dive into the Vulkan specification, experiment with subgroup operations, and try the Vulkan Ray Tracing extensions. Explore Vulkan Video for encode/decode.
Cross‑platform: Consider using SYCL or OpenCL for broader portability, but now you have the low‑level knowledge to understand what's happening under the hood.

The world of GPU compute is vast and exciting. Keep experimenting, profiling, and optimizing. And remember – the best way to learn is to write code. So fire up Xcode or your Vulkan SDK and start computing!

If you found this guide helpful, share it with your fellow developers. Have questions or want to share your own experiences? Leave a comment below. Happy coding!

Advanced GPU Optimization: OpenCL & SYCL from zero to hero - Part 2

Javad — Sat, 14 Feb 2026 16:31:54 +0000

Welcome back, heterogeneous computing enthusiasts! In our previous guides, we covered the fundamentals of OpenCL and SYCL, from setting up your environment to writing basic kernels and optimizing them with local memory. Now it's time to go deeper. In this advanced tutorial, we'll explore:

Key features of the OpenCL specification (2.x/3.0) that enable sophisticated programming.
What SYCL 2020 brings to the table, including unified shared memory and group algorithms.
Implementation of complex parallel patterns: reductions, scans, and stencils.
How to target CPUs, GPUs, and FPGAs with SYCL and compare performance.
A tour of oneAPI's domain‑specific libraries (oneMKL, oneDNN, oneDPL, etc.) that supercharge your development.

Whether you're looking to squeeze every drop of performance from your hardware or write portable code that scales across architectures, this guide has you covered.

Part 1: Deep Dive into the OpenCL Specification

The OpenCL specification has evolved significantly. While OpenCL 1.2 is the most widely supported baseline, modern devices support newer features that can greatly simplify and accelerate development. Let's look at the highlights.

1.1 OpenCL 2.x Features

Shared Virtual Memory (SVM)

OpenCL 2.0 introduced SVM, allowing the host and devices to share pointers. This eliminates explicit buffer mappings and can simplify programming. Three types exist:

Coarse-grained SVM: Sharing at buffer granularity; synchronization needed.
Fine-grained SVM: Sharing at any memory location; optional atomics.
Fine-grained buffer SVM: Like fine-grained but limited to a buffer region.

With SVM, you can pass complex data structures (like linked lists) to kernels without deep copies.

Device-side Enqueue

Kernels can enqueue other kernels, creating nested parallelism. This is great for recursive algorithms or dynamic workloads. The parent kernel creates a child kernel and enqueues it to a device queue.

Pipes

Pipes are FIFO data structures connecting kernels. They facilitate streaming patterns without global memory round trips. One kernel writes to a pipe, another reads from it, with synchronization handled by the hardware.

Work-group Collective Functions

OpenCL 2.0 added built‑in work‑group functions like work_group_reduce_add, work_group_scan_exclusive_add, etc., enabling efficient intra‑group reductions and scans.

C11 Atomics

Kernels can use atomic operations conforming to the C11 memory model, essential for lock‑free algorithms.

1.2 OpenCL 3.0: The New Baseline

OpenCL 3.0, released in 2020, takes a different approach: it makes OpenCL 1.2 the mandatory baseline, with all 2.x features optional. This allows vendors to implement only what they support while maintaining compatibility. Key additions:

A new query mechanism to check for optional features.
Support for newer C++ features in kernels (like auto, lambda expressions) via the cl_khr_c++11 extension.
Better integration with Vulkan and other APIs.

1.3 Using Advanced Features in Practice

To use these features, you must check for extensions and enable them at program creation. For example, to use SVM:

cl_device_svm_capabilities caps;
clGetDeviceInfo(device, CL_DEVICE_SVM_CAPABILITIES, sizeof(caps), &caps, NULL);
if (caps & CL_DEVICE_SVM_COARSE_GRAIN_BUFFER) {
    // Allocate coarse‑grain SVM buffer
    clSVMAlloc(context, CL_MEM_READ_WRITE, size, 0);
}

Device‑side enqueue requires building the program with the -cl-std=CL2.0 flag and using the cl_khr_device_enqueue extension.

While these features are powerful, they are not universally supported across all devices. Always fall back to portable code when necessary.

Part 2: Deep Dive into the SYCL 2020 Specification

SYCL 2020 is a major update that aligns with modern C++ and adds many features inspired by OpenCL 2.x and beyond. Let's explore the most impactful additions.

2.1 Unified Shared Memory (USM)

SYCL 2020 standardizes USM, providing pointer‑based memory management. Three allocation types:

malloc_device: Device‑accessible memory (not directly on host).
malloc_host: Host‑accessible memory (may be cached on device).
malloc_shared: Automatically migrated memory, accessible from both sides.

Example:

queue q;
int *data = malloc_shared<int>(N, q);
q.parallel_for(range<1>(N), [=](id<1> i) { data[i] = i; }).wait();
// data is updated on host after synchronization
free(data, q);

USM simplifies porting CUDA/HIP code and reduces verbosity compared to buffers.

2.2 Reductions

SYCL 2020 introduces built‑in reduction support. You can use reduction objects with parallel_for to perform reductions efficiently.

queue q;
std::vector<float> vec(N, 1.0f);
buffer<float> buf(vec);
float init = 0.0f;
{
    buffer<float> sum_buf(&init, 1);
    q.submit([&](handler& h) {
        auto in = buf.get_access<access::mode::read>(h);
        auto out = sum_buf.get_access<access::mode::write>(h);
        h.parallel_for(range<1>(N), reduction(out, 0.0f, std::plus<>()), [=](id<1> i, auto& sum) {
            sum += in[i];
        });
    });
}
// sum_buf now contains the total

The implementation chooses the most efficient strategy (tree reduction, atomics, etc.) based on the device.

2.3 Group Algorithms and Sub‑Groups

SYCL 2020 provides collective algorithms for work‑groups, similar to OpenCL 2.0:

group_barrier()
group_broadcast()
group_reduce(), group_scan() (exclusive/inclusive)
sub_group class for warp/wavefront operations.

Example of group reduction:

h.parallel_for(nd_range<1>(N, 256), [=](nd_item<1> it) {
    auto g = it.get_group();
    int local_sum = ...;
    int group_sum = reduce_over_group(g, local_sum, std::plus<>());
    if (it.get_local_id(0) == 0)
        atomic_ref<int, memory_order::relaxed, memory_scope::device>(global_sum) += group_sum;
});

Sub‑groups allow portable SIMD programming. You can query sub‑group size and use shuffles.

2.4 In‑Order Queues and Properties

SYCL 2020 adds queue properties to specify in‑order execution, which can simplify dependencies:

queue q{ property::queue::in_order() };

2.5 Interoperability with OpenCL

SYCL 2020 improves interop: you can create a SYCL platform, device, context, or queue from native OpenCL handles, and vice versa. This allows incremental migration or mixing low‑level control with high‑level productivity.

Part 3: Complex Parallel Patterns

Now that we've covered the specifications, let's implement three fundamental parallel patterns: reduction, scan, and stencil. We'll show both OpenCL and SYCL versions, with optimizations.

3.1 Reduction (Sum of an Array)

Reduction combines elements using an associative operator (e.g., sum). A common approach is hierarchical: each work‑group reduces its chunk in local memory, then one thread per group adds to global atomic.

OpenCL kernel:

__kernel void reduce(__global const float *in, __global float *out, __local float *localMem, int n) {
    int gid = get_global_id(0);
    int lid = get_local_id(0);
    int groupSize = get_local_size(0);
    int groupId = get_group_id(0);

    // Load data into local memory
    float sum = (gid < n) ? in[gid] : 0.0f;
    localMem[lid] = sum;
    barrier(CLK_LOCAL_MEM_FENCE);

    // Reduce within group
    for (int stride = groupSize / 2; stride > 0; stride >>= 1) {
        if (lid < stride) {
            localMem[lid] += localMem[lid + stride];
        }
        barrier(CLK_LOCAL_MEM_FENCE);
    }

    // Write partial sum
    if (lid == 0) {
        out[groupId] = localMem[0];
    }
}

After this kernel, you have partial sums. Run a second kernel with one work‑group to combine them, or copy to host and finalize.

SYCL with reductions:
SYCL's built‑in reduction handles this automatically, but for educational purposes, here's a manual version:

q.submit([&](handler& h) {
    auto in = buf_in.get_access<access::mode::read>(h);
    auto out = buf_out.get_access<access::mode::write>(h);
    accessor<float, 1, access::mode::read_write, access::target::local> localMem(range<1>(256), h);

    h.parallel_for(nd_range<1>(n, 256), [=](nd_item<1> it) {
        int gid = it.get_global_id(0);
        int lid = it.get_local_id(0);
        float sum = (gid < n) ? in[gid] : 0.0f;
        localMem[lid] = sum;
        it.barrier();

        for (int stride = 128; stride > 0; stride >>= 1) {
            if (lid < stride) {
                localMem[lid] += localMem[lid + stride];
            }
            it.barrier();
        }

        if (lid == 0) {
            out[it.get_group(0)] = localMem[0];
        }
    });
});

But the SYCL 2020 way is simpler:

q.submit([&](handler& h) {
    auto in = buf_in.get_access<access::mode::read>(h);
    h.parallel_for(range<1>(n), reduction(sum_buf, h, 0.0f, std::plus<>()), [=](id<1> i, auto& sum) {
        sum += in[i];
    });
});

3.2 Scan (Prefix Sum)

Scan computes cumulative sums. An efficient work‑group scan uses a double‑buffered local memory approach.

OpenCL kernel for inclusive scan (one block):

__kernel void scan(__global const float *in, __global float *out, __local float *temp, int n) {
    int gid = get_global_id(0);
    int lid = get_local_id(0);
    int groupSize = get_local_size(0);
    int groupId = get_group_id(0);

    // Load data
    temp[2 * lid] = (gid < n) ? in[gid] : 0;
    temp[2 * lid + 1] = (gid + groupSize < n) ? in[gid + groupSize] : 0;
    barrier(CLK_LOCAL_MEM_FENCE);

    // Upsweep phase (parallel reduction)
    for (int stride = 1; stride < 2 * groupSize; stride <<= 1) {
        int index = (lid + 1) * stride * 2 - 1;
        if (index < 2 * groupSize && lid < groupSize) {
            temp[index] += temp[index - stride];
        }
        barrier(CLK_LOCAL_MEM_FENCE);
    }

    // Downsweep phase
    for (int stride = groupSize; stride > 0; stride >>= 1) {
        int index = (lid + 1) * stride * 2 - 1;
        if (index + stride < 2 * groupSize && lid < groupSize) {
            temp[index + stride] += temp[index];
        }
        barrier(CLK_LOCAL_MEM_FENCE);
    }

    // Write results
    if (gid < n) {
        out[gid] = temp[2 * lid];
        if (gid + groupSize < n)
            out[gid + groupSize] = temp[2 * lid + 1];
    }
}

After this, you need a second pass to add the total of previous groups. For large arrays, you'd combine with a top‑level scan of group sums.

SYCL using group_scan:
SYCL 2020 provides group_scan_inclusive and group_scan_exclusive:

q.submit([&](handler& h) {
    auto in = buf_in.get_access<access::mode::read>(h);
    auto out = buf_out.get_access<access::mode::write>(h);
    accessor<float, 1, access::mode::read_write, access::target::local> localMem(range<1>(256), h);

    h.parallel_for(nd_range<1>(n, 256), [=](nd_item<1> it) {
        int gid = it.get_global_id(0);
        int lid = it.get_local_id(0);
        float val = (gid < n) ? in[gid] : 0.0f;
        auto grp = it.get_group();
        float scan_val = group_scan_inclusive(grp, val, std::plus<>());
        if (gid < n)
            out[gid] = scan_val;
    });
});

This performs an inclusive scan within each work‑group. To get a global scan, you'd again need to handle the group offsets.

3.3 Stencil (2D Heat Conduction)

Stencils update each element based on its neighbors. Optimizations: use local memory to cache a tile plus halo.

OpenCL kernel for 2D heat stencil:

__kernel void heat(__global const float *in, __global float *out, int width, int height) {
    int x = get_global_id(0);
    int y = get_global_id(1);
    int localX = get_local_id(0);
    int localY = get_local_id(1);
    int localSizeX = get_local_size(0);
    int localSizeY = get_local_size(1);

    // Local memory with halo of 1
    __local float tile[LOCAL_SIZE_Y + 2][LOCAL_SIZE_X + 2];

    // Load interior
    if (x < width && y < height)
        tile[localY + 1][localX + 1] = in[y * width + x];
    else
        tile[localY + 1][localX + 1] = 0.0f;

    // Load halo regions (simplified: assumes boundary conditions)
    if (localY == 0 && y > 0)
        tile[0][localX + 1] = in[(y-1) * width + x];
    if (localY == localSizeY - 1 && y < height - 1)
        tile[localSizeY + 1][localX + 1] = in[(y+1) * width + x];
    // similarly for left/right
    barrier(CLK_LOCAL_MEM_FENCE);

    if (x < width && y < height) {
        float center = tile[localY + 1][localX + 1];
        float left   = tile[localY + 1][localX];
        float right  = tile[localY + 1][localX + 2];
        float up     = tile[localY][localX + 1];
        float down   = tile[localY + 2][localX + 1];
        out[y * width + x] = center + 0.125f * (left + right + up + down - 4.0f * center);
    }
}

SYCL version (with local accessor and nd_range):

q.submit([&](handler& h) {
    auto in = buf_in.get_access<access::mode::read>(h);
    auto out = buf_out.get_access<access::mode::write>(h);
    accessor<float, 2, access::mode::read_write, access::target::local> tile(range<2>(LOCAL_Y+2, LOCAL_X+2), h);

    h.parallel_for(nd_range<2>(range<2>(height, width), range<2>(LOCAL_Y, LOCAL_X)), [=](nd_item<2> it) {
        int x = it.get_global_id(1);
        int y = it.get_global_id(0);
        int lx = it.get_local_id(1);
        int ly = it.get_local_id(0);

        // Load interior
        if (x < width && y < height)
            tile[ly+1][lx+1] = in[y][x];

        // Load halos (simplified)
        // ... similar to OpenCL
        it.barrier();

        if (x < width && y < height) {
            float center = tile[ly+1][lx+1];
            float left   = tile[ly+1][lx];
            float right  = tile[ly+1][lx+2];
            float up     = tile[ly][lx+1];
            float down   = tile[ly+2][lx+1];
            out[y][x] = center + 0.125f * (left + right + up + down - 4.0f * center);
        }
    });
});

These patterns are building blocks for many scientific and engineering applications.

Part 4: SYCL on Different Hardware and Performance Comparison

One of SYCL's greatest promises is "write once, run anywhere." But how does performance compare across CPUs, GPUs, and FPGAs? Let's explore the practical considerations.

4.1 Targeting Devices with SYCL

Selecting a device is straightforward:

// Default selector (picks a device, often GPU if available)
queue q;

// Specific selectors
cpu_selector cpu;
gpu_selector gpu;
accelerator_selector fpga; // for FPGA

queue q_cpu(cpu);
queue q_gpu(gpu);
queue q_fpga(fpga);

You can also write custom selectors based on device info (e.g., preferring a GPU with larger local memory).

4.2 Performance Considerations per Device

CPU:

Optimize for vectorization: use wide SIMD (work‑group size = vector width).
Minimize data transfers; CPU shares memory with host, so USM can be efficient.
Avoid excessive local memory; CPU caches are good.

GPU:

Maximize occupancy: choose work‑group size that hides latency (e.g., 128‑256).
Use local memory for data reuse.
Ensure coalesced global memory accesses.

FPGA:

Kernels are synthesized into hardware pipelines.
Loop unrolling, pipelining, and on‑chip memory are key.
SYCL for FPGA supports [[intel::max_work_group_size(1)]] to create single‑work‑item kernels that are highly pipelined.
Use #pragma unroll and [[intel::initiation_interval(1)]] directives.

4.3 Benchmarking and Comparison

To compare performance across devices, you need a common benchmark. Let's take matrix multiplication (tiled) and run it on different hardware.

Methodology:

Compile with appropriate flags: for CPU, use -O3 -march=native; for GPU, -O3 -fsycl; for FPGA, you need the Intel FPGA SDK and use -Xshardware for actual hardware runs.
Measure kernel execution time using SYCL events:

q.submit([&](handler& h) {
    // ...
}).wait();
auto start = std::chrono::steady_clock::now();
q.submit([&](handler& h) {
    // kernel
}).wait();
auto end = std::chrono::steady_clock::now();

Better: use SYCL profiling info via property::queue::enable_profiling().

Expected Observations:

CPU: Good for small matrices, latency may be higher due to OS scheduling.
GPU: Excels at large, regular computations.
FPGA: Often slower for compute‑bound tasks but can be very power‑efficient and excel at pipeline‑friendly streaming.

Case Study: Matrix Multiplication
We can run the tiled SYCL kernel from Part 2.7 on a CPU, GPU, and FPGA (if available). The performance will vary dramatically. For example, on an Intel Core i9, the CPU might achieve ~50 GFLOPs; on an Intel UHD Graphics, maybe ~200 GFLOPs; on an Arria 10 FPGA, perhaps ~50 GFLOPs but with lower power. However, the FPGA could achieve lower latency for small batches and better predictability.

4.4 Portable Optimization Tips

Use #ifdef SYCL_DEVICE_CPU etc. to specialize for each architecture.
Tune work‑group size dynamically using device.get_info<info::device::max_work_group_size>().
For FPGA, design kernels to be highly parallel with pipeline parallelism rather than data parallelism.

Part 5: oneAPI Domain-Specific Libraries

oneAPI includes a rich set of libraries that provide optimized building blocks for common tasks. Using them can drastically accelerate development while ensuring high performance across devices.

5.1 oneMKL (Math Kernel Library)

oneMKL provides BLAS, LAPACK, FFT, RNG, and more. It has SYCL interfaces that accept SYCL buffers or USM pointers.

Example: Matrix multiplication using oneMKL:

#include <oneapi/mkl/blas.hpp>
queue q;
float alpha = 1.0f, beta = 0.0f;
oneapi::mkl::blas::gemm(q, oneapi::mkl::transpose::nontrans, oneapi::mkl::transpose::nontrans,
                        N, N, N, alpha, d_A, N, d_B, N, beta, d_C, N);
q.wait();

This will call an optimized GEMM for the target device (CPU, GPU, or FPGA). No need to write your own kernel!

FFT example:

#include <oneapi/mkl/dft.hpp>
using namespace oneapi::mkl::dft;
descriptor<precision::single, domain::real> transform(N);
transform.commit(q);
// in and out buffers...
transform.compute_forward(in, out);

5.2 oneDNN (Deep Neural Networks)

oneDNN (formerly DNNL) provides primitives for deep learning: convolution, pooling, normalization, etc. It has SYCL and OpenCL backends.

Example: Convolution:

#include <oneapi/dnnl/dnnl.hpp>
using namespace dnnl;

engine eng(engine::kind::gpu, 0);
stream s(eng);
memory::dims conv_src_tz = {N, C, H, W};
memory::dims conv_weights_tz = {K, C, R, S};
memory::dims conv_dst_tz = {N, K, H', W'};
// Create memory objects
auto conv_src_mem = memory({{conv_src_tz}, memory::data_type::f32, memory::format_tag::any}, eng);
// Initialize convolution descriptor and primitive
convolution_forward::desc conv_desc(prop_kind::forward, algorithm::convolution_direct,
                                    conv_src_md, conv_weights_md, conv_dst_md,
                                    strides, padding_l, padding_r);
convolution_forward::primitive_desc conv_pd(conv_desc, eng);
auto conv = convolution_forward(conv_pd);
conv.execute(s, {{DNNL_ARG_SRC, conv_src_mem}, {DNNL_ARG_WEIGHTS, conv_weights_mem}, {DNNL_ARG_DST, conv_dst_mem}});
s.wait();

oneDNN automatically selects the best implementation for the hardware.

5.3 oneDPL (DPC++ Library)

oneDPL provides parallel algorithms similar to C++17's Parallel STL, implemented with SYCL.

Example: parallel sort:

#include <oneapi/dpl/algorithm>
#include <oneapi/dpl/execution>

std::vector<float> data(N);
// fill data...
auto policy = oneapi::dpl::execution::make_device_policy(q);
std::sort(policy, data.begin(), data.end());

Other algorithms: reduce, transform, scan, etc.

5.4 oneTBB (Threading Building Blocks)

While not exclusively for devices, oneTBB integrates with oneAPI to provide task parallelism on the host, which can coordinate with device work.

5.5 Other oneAPI Libraries

oneVPL – Video Processing Library for encode/decode.
oneDAL – Data Analytics Library.
oneCCL – Collective Communications Library for distributed computing.

5.6 Combining Libraries in a Project

A typical SYCL + oneAPI project might use:

oneDPL for data preparation on CPU.
oneMKL for FFT on GPU.
oneDNN for inference on FPGA.
oneTBB to manage host-side tasks.

All with the same SYCL queue and memory model.

Conclusion

We've journeyed through the advanced features of OpenCL and SYCL, implemented essential parallel patterns, explored performance across diverse hardware, and surveyed the powerful oneAPI libraries. Armed with this knowledge, you can now tackle complex computational problems with confidence and write code that runs efficiently on CPUs, GPUs, and FPGAs alike.

Remember, the key to mastering heterogeneous computing is practice. Experiment with these patterns on your own hardware, profile them, and see where optimizations take you. And when you need a productivity boost, don't reinvent the wheel – leverage the oneAPI libraries.

Further Resources:

Thank you for reading! If you have questions or want to share your experiences with OpenCL, SYCL, or oneAPI, leave a comment below. Happy coding!

Advanced GPU Optimization: OpenCL & SYCL from zero to hero

Javad — Sat, 14 Feb 2026 16:29:50 +0000

Welcome to the world of heterogeneous programming! In this comprehensive guide, we’ll explore two powerful frameworks that let you harness the computational power of CPUs, GPUs, FPGAs, and other accelerators: OpenCL and SYCL (the heart of Intel’s oneAPI initiative). Whether you’re a student, researcher, or industry developer, by the end of this tutorial you’ll be able to write portable, efficient code that runs on a variety of hardware.

We’ll start from the absolute basics (no prior experience required) and progress to advanced topics like local memory optimization, profiling, and multi-device programming. Along the way, you’ll see how OpenCL provides a low‑level, C‑based interface for fine‑grained control, while SYCL offers a modern C++ single‑source approach that builds on OpenCL’s concepts.

Let’s dive in!

Introduction: Why OpenCL and SYCL?

OpenCL (Open Computing Language) is the industry standard for heterogeneous computing. It allows you to write kernels (programs) that execute on any OpenCL‑compliant device—CPUs, GPUs, DSPs, FPGAs—using a C‑based language. It gives you explicit control over memory and execution, but at the cost of some verbosity.

SYCL (pronounced “sickle”) is a higher‑level C++ abstraction built on top of OpenCL. It enables single‑source programming: you write host and device code in the same C++ file, using modern C++ features. SYCL is the foundation of Intel’s oneAPI initiative, which aims to provide a unified programming model across CPUs, GPUs, and FPGAs. SYCL code can be compiled for different backends (OpenCL, Level Zero, CUDA) without modification.

Both are essential tools in the heterogeneous computing landscape. Learning them gives you the flexibility to target a wide range of hardware while keeping your code portable.

Part 1: OpenCL Fundamentals

1.1 Setting Up OpenCL

Before we write code, we need an OpenCL development environment. OpenCL consists of a host API (to manage devices, memory, and execution) and a kernel language (based on C99 with extensions).

Installation:

On Linux: Install the OpenCL headers and a runtime (e.g., from your GPU vendor). For Intel CPUs/GPUs, install the Intel Compute Runtime. For NVIDIA, install the CUDA toolkit which includes OpenCL. For AMD, install ROCm or the AMD APP SDK.
On Windows: Download and install the appropriate SDK from your hardware vendor (Intel, NVIDIA, AMD). You can also use the open‑source PoCL (Portable Computing Language) for CPU development.

Basic includes and linking:

#include <CL/cl.h>

Link against the OpenCL library (e.g., -lOpenCL on Linux, OpenCL.lib on Windows).

1.2 The OpenCL Execution Model

OpenCL organizes work into an NDRange (N‑dimensional range). The key concepts:

Work‑item: The smallest unit of execution. Each work‑item executes the same kernel but works on different data.
Work‑group: A collection of work‑items that execute together on a compute unit. Work‑items in the same group can synchronize via barriers and share memory through local memory.
NDRange: The entire grid of work‑items, partitioned into work‑groups. You specify the global size (total work‑items) and local size (work‑group size).

When you launch a kernel, you provide the number of dimensions (1, 2, or 3), the global size in each dimension, and optionally the local size. The runtime then schedules work‑groups onto compute units.

1.3 OpenCL Memory Model

OpenCL defines several memory regions:

Global memory: Accessible by all work‑items; large but high latency.
Constant memory: Read‑only global memory, cached.
Local memory: Shared within a work‑group; fast, programmer‑managed.
Private memory: Per work‑item, typically registers.

Host code manages device memory via buffer objects. Data is explicitly moved between host and device.

1.4 Your First OpenCL Program: Vector Addition

Let’s write a simple vector addition kernel and the host code to run it.

Kernel (stored in a string or separate .cl file):

__kernel void vecAdd(__global const float *a,
                     __global const float *b,
                     __global float *c) {
    int i = get_global_id(0);
    c[i] = a[i] + b[i];
}

Host code (C):

#include <stdio.h>
#include <stdlib.h>
#include <CL/cl.h>

#define VEC_SIZE 1048576  // 1M elements

int main() {
    // Step 1: Get platform and device
    cl_platform_id platform;
    cl_device_id device;
    clGetPlatformIDs(1, &platform, NULL);
    clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL);

    // Step 2: Create context and command queue
    cl_context context = clCreateContext(NULL, 1, &device, NULL, NULL, NULL);
    cl_command_queue queue = clCreateCommandQueue(context, device, 0, NULL);

    // Step 3: Allocate host memory and initialize data
    size_t bytes = VEC_SIZE * sizeof(float);
    float *h_a = (float*)malloc(bytes);
    float *h_b = (float*)malloc(bytes);
    float *h_c = (float*)malloc(bytes);
    for (int i = 0; i < VEC_SIZE; i++) {
        h_a[i] = i * 1.0f;
        h_b[i] = (i + 1) * 1.0f;
    }

    // Step 4: Create device buffers
    cl_mem d_a = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, bytes, h_a, NULL);
    cl_mem d_b = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, bytes, h_b, NULL);
    cl_mem d_c = clCreateBuffer(context, CL_MEM_WRITE_ONLY, bytes, NULL, NULL);

    // Step 5: Build the program
    const char *kernelSource = "__kernel void vecAdd(__global const float *a, __global const float *b, __global float *c) { int i = get_global_id(0); c[i] = a[i] + b[i]; }";
    cl_program program = clCreateProgramWithSource(context, 1, &kernelSource, NULL, NULL);
    clBuildProgram(program, 1, &device, NULL, NULL, NULL);
    cl_kernel kernel = clCreateKernel(program, "vecAdd", NULL);

    // Step 6: Set kernel arguments
    clSetKernelArg(kernel, 0, sizeof(cl_mem), &d_a);
    clSetKernelArg(kernel, 1, sizeof(cl_mem), &d_b);
    clSetKernelArg(kernel, 2, sizeof(cl_mem), &d_c);

    // Step 7: Execute kernel
    size_t global_size = VEC_SIZE;
    size_t local_size = 256;  // must divide global_size
    clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &global_size, &local_size, 0, NULL, NULL);

    // Step 8: Read result back
    clEnqueueReadBuffer(queue, d_c, CL_TRUE, 0, bytes, h_c, 0, NULL, NULL);

    // Step 9: Verify (simplified)
    for (int i = 0; i < 10; i++)
        printf("%f + %f = %f\n", h_a[i], h_b[i], h_c[i]);

    // Step 10: Clean up
    clReleaseMemObject(d_a); clReleaseMemObject(d_b); clReleaseMemObject(d_c);
    clReleaseKernel(kernel); clReleaseProgram(program);
    clReleaseCommandQueue(queue); clReleaseContext(context);
    free(h_a); free(h_b); free(h_c);

    return 0;
}

Compile and run (on Linux):

gcc -o vecAdd vecAdd.c -lOpenCL
./vecAdd

Explanation:

We query a platform and device (here, we take the first GPU).
Create a context and command queue (the queue is where we enqueue commands).
Allocate host buffers and device buffers. clCreateBuffer with CL_MEM_COPY_HOST_PTR copies data immediately.
Build the program from source string; in real applications, you’d load from a file and check build errors.
Set kernel arguments, enqueue the kernel with global and local sizes.
Read the result back to host with clEnqueueReadBuffer.
Clean up OpenCL objects.

This is the basic pattern for any OpenCL program.

1.5 Building and Running OpenCL Programs

Always check error codes. Each OpenCL function returns a cl_int error code. Use macros like:

#define CL_CHECK(err) \
    if (err != CL_SUCCESS) { \
        fprintf(stderr, "OpenCL error %d at %s:%d\n", err, __FILE__, __LINE__); \
        exit(EXIT_FAILURE); \
    }

// Usage:
cl_int err;
cl_context context = clCreateContext(NULL, 1, &device, NULL, NULL, &err);
CL_CHECK(err);

When building a program, you should check build logs:

clBuildProgram(program, 1, &device, NULL, NULL, NULL);
if (err != CL_SUCCESS) {
    size_t log_size;
    clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, 0, NULL, &log_size);
    char *log = (char*)malloc(log_size);
    clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, log_size, log, NULL);
    fprintf(stderr, "Build log:\n%s\n", log);
    free(log);
    exit(EXIT_FAILURE);
}

1.6 Profiling OpenCL with Events

OpenCL events can measure execution time. Enable profiling on the command queue:

cl_command_queue queue = clCreateCommandQueue(context, device, CL_QUEUE_PROFILING_ENABLE, &err);

Then attach an event to kernel launch and read timestamps:

cl_event event;
clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &global_size, &local_size, 0, NULL, &event);
clWaitForEvents(1, &event);

cl_ulong start, end;
clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &start, NULL);
clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_END, sizeof(cl_ulong), &end, NULL);
double time_ms = (end - start) / 1.0e6;
printf("Kernel execution time: %f ms\n", time_ms);

clReleaseEvent(event);

1.7 Example: Matrix Multiplication (Naive and Tiled)

Let's implement matrix multiplication to demonstrate local memory usage.

Naive kernel (each work-item computes one element of C):

__kernel void matMulNaive(__global const float *A,
                          __global const float *B,
                          __global float *C,
                          int N) {
    int row = get_global_id(1);
    int col = get_global_id(0);
    if (row < N && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < N; k++) {
            sum += A[row * N + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

Host side (simplified):

size_t global_size[2] = {N, N};
clEnqueueNDRangeKernel(queue, kernel, 2, NULL, global_size, NULL, 0, NULL, NULL);

Tiled kernel using local memory:

#define TILE_SIZE 16

__kernel void matMulTiled(__global const float *A,
                          __global const float *B,
                          __global float *C,
                          int N) {
    __local float As[TILE_SIZE][TILE_SIZE];
    __local float Bs[TILE_SIZE][TILE_SIZE];

    int row = get_global_id(1);
    int col = get_global_id(0);
    int localRow = get_local_id(1);
    int localCol = get_local_id(0);

    float sum = 0.0f;
    for (int tile = 0; tile < N / TILE_SIZE; tile++) {
        // Load tile of A and B into local memory
        As[localRow][localCol] = A[row * N + (tile * TILE_SIZE + localCol)];
        Bs[localRow][localCol] = B[(tile * TILE_SIZE + localRow) * N + col];
        barrier(CLK_LOCAL_MEM_FENCE);

        // Compute partial sum
        for (int k = 0; k < TILE_SIZE; k++) {
            sum += As[localRow][k] * Bs[k][localCol];
        }
        barrier(CLK_LOCAL_MEM_FENCE);
    }

    if (row < N && col < N) {
        C[row * N + col] = sum;
    }
}

Launch configuration:

size_t local_size[2] = {TILE_SIZE, TILE_SIZE};
size_t global_size[2] = {N, N};
clEnqueueNDRangeKernel(queue, kernel, 2, NULL, global_size, local_size, 0, NULL, NULL);

This tiled version reduces global memory traffic by reusing data loaded into local memory. It's a key optimization pattern.

1.8 Advanced OpenCL Topics

Multiple devices: You can create a context with multiple devices and use separate command queues.
Image objects: OpenCL supports 2D/3D images with samplers for graphics‑oriented workloads.
SVM (Shared Virtual Memory): Allows sharing pointers between host and device (OpenCL 2.0).
Device-side enqueue: Kernels can enqueue other kernels (OpenCL 2.0).

Part 2: SYCL (oneAPI) Fundamentals

2.1 What is SYCL?

SYCL is a C++ abstraction layer that sits on top of OpenCL (and other backends). It allows you to write host and device code in the same file using standard C++17 (or later). Key features:

Single-source: Kernel code is written as C++ functors or lambdas, compiled by a SYCL compiler.
High-level abstractions: Buffers and accessors manage memory dependencies; command groups encapsulate kernel launches.
Portability: Same code can target CPUs, GPUs, FPGAs, etc., via different SYCL implementations (Intel oneAPI DPC++, ComputeCpp, hipSYCL).
Interoperability: Can mix with OpenCL code if needed.

Intel’s oneAPI is a complete suite of tools and libraries built around SYCL (DPC++ – Data Parallel C++). We'll focus on DPC++ as the most widely adopted SYCL implementation.

2.2 Setting Up oneAPI / SYCL

Option 1: Intel oneAPI base toolkit (free):

Download from Intel oneAPI.
Install and source the environment: source /opt/intel/oneapi/setvars.sh.

Option 2: Open-source DPC++:

Follow instructions on intel/llvm to build the DPC++ compiler.

Option 3: hipSYCL (supports NVIDIA and AMD backends):

See hipSYCL documentation.

We'll assume you have dpcpp (DPC++ compiler) available.

2.3 SYCL Execution Model

SYCL programs revolve around:

queue: Submits work to a device (similar to OpenCL command queue).
buffer: Manages data across host and device.
accessor: Requests access to a buffer inside a command group; specifies access mode (read, write, read_write) and target (host, device).
handler: Used inside command groups to set kernel arguments and launch parallel work.
nd_range: Defines global and local size (like OpenCL NDRange).
parallel_for: Launches a kernel over a range.

A typical SYCL program:

#include <CL/sycl.hpp>
using namespace sycl;

int main() {
    queue q;  // selects default device

    buffer<float, 1> a_buf{ h_a, range<1>{N} };
    buffer<float, 1> b_buf{ h_b, range<1>{N} };
    buffer<float, 1> c_buf{ range<1>{N} };

    q.submit([&](handler& h) {
        auto a = a_buf.get_access<access::mode::read>(h);
        auto b = b_buf.get_access<access::mode::read>(h);
        auto c = c_buf.get_access<access::mode::write>(h);

        h.parallel_for(range<1>{N}, [=](id<1> i) {
            c[i] = a[i] + b[i];
        });
    });

    // Accessor on host to read result
    auto c_host = c_buf.get_access<access::mode::read>();
    // ... use c_host ...
}

2.4 Memory Model: Buffers and Accessors

Buffers abstract data movement. Accessors create dependencies: when you request an accessor in a command group, the runtime ensures data is available on the device before the kernel runs, and updates the host after.

USM (Unified Shared Memory) is an alternative that uses pointers, similar to CUDA/HIP. SYCL supports device USM, shared USM, and host USM.

2.5 Your First SYCL Program: Vector Addition

Let's rewrite the vector addition using SYCL.

#include <CL/sycl.hpp>
#include <iostream>
#include <vector>

using namespace sycl;

int main() {
    const int N = 1 << 20;
    std::vector<float> a(N), b(N), c(N);
    for (int i = 0; i < N; i++) {
        a[i] = i * 1.0f;
        b[i] = (i + 1) * 1.0f;
    }

    try {
        queue q;  // default device

        // Create buffers from host data (moves data to device automatically)
        buffer<float, 1> a_buf(a.data(), range<1>(N));
        buffer<float, 1> b_buf(b.data(), range<1>(N));
        buffer<float, 1> c_buf(c.data(), range<1>(N));

        // Submit command group
        q.submit([&](handler& h) {
            // Request accessors
            auto a_acc = a_buf.get_access<access::mode::read>(h);
            auto b_acc = b_buf.get_access<access::mode::read>(h);
            auto c_acc = c_buf.get_access<access::mode::write>(h);

            // Launch kernel
            h.parallel_for(range<1>(N), [=](id<1> i) {
                c_acc[i] = a_acc[i] + b_acc[i];
            });
        });

        // Buffers go out of scope, so data is copied back to host automatically
        // (destructor of c_buf waits for completion and updates host)
    } catch (sycl::exception& e) {
        std::cerr << "SYCL exception: " << e.what() << std::endl;
        return 1;
    }

    // Verify
    for (int i = 0; i < 10; i++)
        std::cout << a[i] << " + " << b[i] << " = " << c[i] << std::endl;

    return 0;
}

Compile with:

dpcpp -o vecAdd vecAdd.cpp
./vecAdd

Explanation:

queue q creates a queue to the default device. You can select a specific device using device_selector.
Buffers wrap host data; they don't copy immediately but manage data lazily.
Inside the command group, accessors specify how the kernel uses the data. The runtime ensures data is on the device before execution.
The kernel is a lambda capturing accessors by value. It executes in parallel over a 1D range.
When buffers are destroyed (at end of scope), they wait for any dependent operations and copy data back to host if needed.

2.6 Error Handling and Asynchronous Exceptions

SYCL uses exceptions for synchronous errors. Asynchronous errors (during kernel execution) can be caught via an async handler:

void async_handler(sycl::exception_list e) {
    for (auto& ep : e) {
        try { std::rethrow_exception(ep); }
        catch (sycl::exception& e) { std::cerr << "Async SYCL exception: " << e.what() << std::endl; }
    }
}

queue q(device_selector{}, async_handler);

2.7 Example: Matrix Multiplication with SYCL

Now let's implement matrix multiplication using SYCL's nd_range for work‑groups and local memory.

#include <CL/sycl.hpp>
#include <iostream>
#include <vector>

using namespace sycl;

#define TILE_SIZE 16

int main() {
    const int N = 1024;
    std::vector<float> A(N * N), B(N * N), C(N * N, 0.0f);
    // Initialize A, B...

    try {
        queue q;

        buffer<float, 2> a_buf(A.data(), range<2>(N, N));
        buffer<float, 2> b_buf(B.data(), range<2>(N, N));
        buffer<float, 2> c_buf(C.data(), range<2>(N, N));

        q.submit([&](handler& h) {
            auto a = a_buf.get_access<access::mode::read>(h);
            auto b = b_buf.get_access<access::mode::read>(h);
            auto c = c_buf.get_access<access::mode::write>(h);

            // Local memory accessor
            accessor<float, 2, access::mode::read_write, access::target::local> As(TILE_SIZE, TILE_SIZE, h);
            accessor<float, 2, access::mode::read_write, access::target::local> Bs(TILE_SIZE, TILE_SIZE, h);

            h.parallel_for(nd_range<2>({N, N}, {TILE_SIZE, TILE_SIZE}), [=](nd_item<2> it) {
                int row = it.get_global_id(1);
                int col = it.get_global_id(0);
                int localRow = it.get_local_id(1);
                int localCol = it.get_local_id(0);

                float sum = 0.0f;
                for (int tile = 0; tile < N / TILE_SIZE; ++tile) {
                    // Load tiles
                    As[localRow][localCol] = a[row][tile * TILE_SIZE + localCol];
                    Bs[localRow][localCol] = b[tile * TILE_SIZE + localRow][col];
                    it.barrier(access::fence_space::local_space);

                    // Compute
                    for (int k = 0; k < TILE_SIZE; ++k) {
                        sum += As[localRow][k] * Bs[k][localCol];
                    }
                    it.barrier(access::fence_space::local_space);
                }

                if (row < N && col < N) {
                    c[row][col] = sum;
                }
            });
        });
    } catch (sycl::exception& e) {
        std::cerr << e.what() << std::endl;
        return 1;
    }

    // Use C...
    return 0;
}

Notes:

We use 2D buffers with range<2>.
Local memory is allocated via accessor with target::local.
nd_range combines global and local sizes.
Barriers are called on the nd_item to synchronize work‑group.

2.8 Unified Shared Memory (USM)

USM simplifies memory management by using pointers. There are three types:

device: Allocated on device, accessible only from device.
host: Allocated on host, accessible only from host.
shared: Managed pointer that can be accessed from both host and device (migrates data automatically).

Example with shared USM:

queue q;
float *a = malloc_shared<float>(N, q);
float *b = malloc_shared<float>(N, q);
float *c = malloc_shared<float>(N, q);

// Initialize on host
for (int i = 0; i < N; i++) {
    a[i] = i * 1.0f;
    b[i] = (i + 1) * 1.0f;
}

q.parallel_for(range<1>(N), [=](id<1> i) {
    c[i] = a[i] + b[i];
}).wait();

// c is ready on host
// ...

free(a, q); free(b, q); free(c, q);

USM is more convenient for porting CUDA/HIP code but requires explicit waits or dependencies.

2.9 Advanced SYCL Topics

Device selection: Use default_selector, gpu_selector, cpu_selector, or custom.
Kernel programming model: SYCL supports basic parallel_for, nd_range, and single_task.
Reductions: SYCL 2020 adds built‑in reductions.
Graph dependencies: Use depends_on and events.
Interoperability with OpenCL: SYCL can create buffers from cl_mem, etc.

Part 3: Advanced Topics and Optimization

3.1 Work-Group Size and Occupancy

Both OpenCL and SYCL rely on work‑groups (or blocks) to map to hardware compute units. Choosing the right local size is crucial.

OpenCL: Local size must divide global size evenly (or use clEnqueueNDRangeKernel with NULL local size to let runtime choose). Typically, use multiples of the preferred work‑group size from device info.
SYCL: Local size specified in nd_range. You can query device info: device.get_info<info::device::max_work_group_size>().

A good rule of thumb: start with a work‑group size of 64–256 and tune based on profiling.

3.2 Local Memory and Barriers

OpenCL: Use __local pointers or __local arrays. Synchronize with barrier(CLK_LOCAL_MEM_FENCE).
SYCL: Use local accessors. Synchronize with item.barrier(access::fence_space::local_space).

Local memory is key for data reuse (e.g., tiled matrix multiplication) and reduction.

3.3 Memory Access Patterns

Coalescing: Work‑items in a work‑group should access consecutive global memory addresses to maximize bandwidth. In OpenCL, this means the global id (or get_global_id(0)) should be the fastest‑changing dimension.
Avoiding bank conflicts: In local memory, ensure that within a work‑group, different work‑items access different banks (or all access the same address for broadcast). This may require padding.

3.4 Using SIMD and Vectorization

OpenCL kernels can use vector types (float4, int4, etc.) to explicitly use SIMD. The compiler may also auto‑vectorize. In SYCL, you can use the vec class.

3.5 Profiling and Debugging

OpenCL:

Use events as shown earlier.
Vendor tools: Intel VTune, AMD CodeXL, NVIDIA Visual Profiler (with OpenCL support).

SYCL:

oneAPI provides Intel VTune Profiler and Intel Advisor.
The queue can be constructed with property::queue::enable_profiling() to get kernel times via events.

Debugging:

OpenCL: printf inside kernels works (if supported). Use clBuildProgram logs.
SYCL: Use std::cout in kernels (SYCL 2020 supports device iostream emulation, but it's slow). For serious debugging, use Intel's Debugger for OpenCL or gdb‑oneapi.

3.6 Interoperability: OpenCL and SYCL

SYCL can interoperate with OpenCL objects:

Create SYCL platform/device from cl_platform_id/cl_device_id.
Create SYCL queue from cl_command_queue.
Create SYCL buffer from cl_mem.

This allows incremental porting or mixing low‑level control with high‑level productivity.

3.7 Targeting Different Devices

OpenCL: Code can run on any OpenCL device, but performance tuning may be device‑specific.
SYCL with oneAPI: You can target CPU, GPU, FPGA with the same source. For FPGA, SYCL supports pipeline‑oriented optimizations via [[intel::kernel_args_restrict]] and other attributes.

Example: selecting a GPU device explicitly:

std::vector<device> devices = device::get_devices(info::device_type::gpu);
queue q(devices[0]);

Part 4: Portability and Choosing Between OpenCL and SYCL

4.1 OpenCL Strengths and Weaknesses

Strengths:

Mature, widely supported.
Fine‑grained control over devices and memory.
C‑based kernel language (easy to learn for C programmers).

Weaknesses:

Verbose host code.
No single‑source, so host and device code are separate.
Kernel language lacks modern C++ features.

4.2 SYCL Strengths

Strengths:

Single‑source C++ with lambdas and templates.
RAII and automatic dependency management via buffers.
Portability across backends (OpenCL, Level Zero, CUDA).
Backed by oneAPI ecosystem with optimized libraries.

Weaknesses:

Newer, still evolving (SYCL 2020 is a major step).
Compiler support not as universal as OpenCL.
Learning curve for the abstraction model.

4.3 When to Use Which

Use OpenCL if you need to support very old hardware, have strict control requirements, or are maintaining an existing OpenCL codebase.
Use SYCL/oneAPI for new projects where productivity and portability are priorities, especially if you target Intel GPUs/FPGAs or want to keep options open for AMD/NVIDIA via hipSYCL.

Conclusion

You've now journeyed through the fundamentals of both OpenCL and SYCL, from setting up your environment to writing optimized kernels for vector addition and matrix multiplication. You've seen how OpenCL gives you explicit control, while SYCL provides a modern C++ interface that simplifies many aspects of heterogeneous programming.

Both frameworks are powerful tools in the developer's toolbox. As you continue, explore real‑world applications, contribute to open‑source projects, and leverage the rich ecosystem of libraries (oneMKL, oneDNN, etc.). The world of heterogeneous computing is vast, and you're now equipped to navigate it.

Next Steps:

Dive into the OpenCL specification and SYCL 2020 specification.
Experiment with more complex patterns: reductions, scans, stencils.
Try SYCL on different hardware (CPU, GPU, FPGA) and compare performance.
Explore oneAPI's domain‑specific libraries for deep learning, ray tracing, etc.

If you found this guide valuable, share it with your peers. Have questions or want to share your experiences? Leave a comment below. Happy coding!

Advanced GPU Optimization: CUDA & HIP from zero to hero

Javad — Sat, 14 Feb 2026 16:23:19 +0000

Welcome, fellow developer! If you’ve ever wanted to tap into the immense power of graphics processors for general-purpose computing, you’re in the right place. GPU programming used to be an NVIDIA-only game with CUDA, but the rise of AMD GPUs in workstations, servers, and supercomputers has made vendor lock‑in a thing of the past. Enter HIP – a C++ runtime API that lets you write code that runs on both NVIDIA and AMD GPUs with minimal changes.

In this comprehensive guide, we’ll start from the very basics (no prior GPU experience assumed) and work our way up to advanced topics like streams, shared memory optimizations, and multi‑GPU programming. Along the way, you’ll learn CUDA inside out and then discover how HIP gives you portability without sacrificing performance. By the end, you’ll be comfortable writing efficient, portable GPU kernels and ready to tackle real‑world parallel problems.

Let’s dive in!

Part 1: CUDA Fundamentals

1.1 Setting Up CUDA

Before we write any code, we need a working CUDA development environment.

Check your hardware: You need an NVIDIA GPU with Compute Capability 3.0 or higher. Run lspci | grep -i nvidia on Linux or check Device Manager on Windows.
Install the CUDA Toolkit:
- Download the appropriate installer from NVIDIA CUDA Toolkit.
- Follow the installation guide for your OS.
Verify installation:

   nvcc --version

You should see the CUDA version. Also run nvidia-smi to check driver and GPU status.

Set environment variables (if not done automatically):

   export PATH=/usr/local/cuda/bin:$PATH
   export LD_LIBRARY_PATH=/usr/local/cuda/lib64:$LD_LIBRARY_PATH

Now you’re ready to compile CUDA code with nvcc.

1.2 The CUDA Programming Model

A GPU contains hundreds or thousands of cores. To use them effectively, we organize work into a hierarchy:

Thread: The smallest unit of execution. Each thread runs the same kernel function but works on different data.
Block: A group of threads that can cooperate via shared memory and synchronize with __syncthreads(). Threads in the same block run on the same streaming multiprocessor (SM).
Grid: A collection of blocks that together cover the entire problem.

When you launch a kernel, you specify the grid and block dimensions (in 1D, 2D, or 3D). Each thread gets unique indices:

threadIdx.x, threadIdx.y, threadIdx.z – position inside its block.
blockIdx.x, blockIdx.y, blockIdx.z – position of the block in the grid.
blockDim.x, blockDim.y, blockDim.z – size of each block.
gridDim.x, gridDim.y, gridDim.z – size of the grid.

The typical formula to map a thread to a data element is:

int i = blockIdx.x * blockDim.x + threadIdx.x;

1.3 Memory Hierarchy

Global memory: Large (several GB), high latency, accessible by all threads. It persists across kernel launches.
Shared memory: Small (tens of KB per block), low latency, shared among threads in the same block. Programmer‑managed cache.
Registers: Private to each thread, fastest.
Constant memory: Read‑only, cached, good for broadcast data.
Texture memory: Specialized for 2D spatial locality.

Understanding this hierarchy is key to writing fast kernels.

1.4 Your First CUDA Kernel: Vector Addition

Let’s implement a simple vector addition. This will introduce all the essential steps.

#include <iostream>
#include <cuda_runtime.h>

// Kernel definition
__global__ void vecAdd(const float *a, const float *b, float *c, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {  // Guard against extra threads
        c[i] = a[i] + b[i];
    }
}

int main() {
    int n = 1 << 20;  // 1,048,576 elements
    size_t bytes = n * sizeof(float);

    // Allocate host memory
    float *h_a = new float[n];
    float *h_b = new float[n];
    float *h_c = new float[n];

    // Initialize input arrays
    for (int i = 0; i < n; i++) {
        h_a[i] = i * 1.0f;
        h_b[i] = (i + 1) * 1.0f;
    }

    // Allocate device memory
    float *d_a, *d_b, *d_c;
    cudaMalloc(&d_a, bytes);
    cudaMalloc(&d_b, bytes);
    cudaMalloc(&d_c, bytes);

    // Copy data from host to device
    cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice);

    // Launch kernel
    int threadsPerBlock = 256;
    int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;
    vecAdd<<<blocksPerGrid, threadsPerBlock>>>(d_a, d_b, d_c, n);

    // Copy result back to host
    cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);

    // Verify (simple check)
    for (int i = 0; i < n; i++) {
        if (h_c[i] != h_a[i] + h_b[i]) {
            std::cout << "Mismatch at index " << i << "\n";
            break;
        }
    }

    // Cleanup
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);
    delete[] h_a;
    delete[] h_b;
    delete[] h_c;

    return 0;
}

Explanation:

__global__ marks a function as a kernel – callable from host and executed on device.
We compute a global index i and, if within bounds, perform the addition.
<<<blocksPerGrid, threadsPerBlock>>> sets the execution configuration.
cudaMalloc allocates memory on the GPU.
cudaMemcpy transfers data between host and device (the last argument specifies direction).
Always free device memory with cudaFree to avoid leaks.

Compile with:

nvcc -o vecAdd vecAdd.cu
./vecAdd

If everything works, you should see no output (or a success message if you add one).

1.5 Understanding Kernel Launch Parameters

The launch configuration <<<grid, block, sharedMem, stream>>> has two required parameters and two optional:

grid (dim3): number of blocks per grid. Can be 1D, 2D, or 3D. E.g., dim3 grid(16, 16) creates a 16x16 grid (256 blocks).
block (dim3): number of threads per block. Usually 1D, but can be 2D/3D. Total threads per block ≤ 1024 (for modern GPUs).
sharedMem (size_t): dynamic shared memory per block (default 0).
stream (cudaStream_t): specifies which stream to enqueue the kernel in (default 0, the null stream).

In our example we used 1D grids and blocks. For 2D data (like images), you might use 2D blocks:

dim3 block(16, 16);
dim3 grid((width + block.x - 1) / block.x, (height + block.y - 1) / block.y);
kernel<<<grid, block>>>(...);

Inside the kernel, you would then compute:

int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
int idx = y * width + x;

1.6 Memory Management in Depth

cudaMalloc(void **ptr, size_t size) – allocates linear memory on device.
cudaFree(void *ptr) – frees device memory.
cudaMemcpy(void *dst, const void *src, size_t count, cudaMemcpyKind kind) – synchronously copies memory.
cudaMemcpyAsync – asynchronous version, requires pinned host memory.
cudaMallocHost – allocates page‑locked (pinned) host memory, which enables faster transfers and asynchronous copies.
cudaMemset – sets device memory to a value.

Pinned memory example:

float *h_a;
cudaMallocHost(&h_a, bytes);   // pinned
cudaMemcpyAsync(d_a, h_a, bytes, cudaMemcpyHostToDevice, stream);

Use pinned memory when you have many transfers or overlapping with computation.

1.7 Error Handling

CUDA functions return an error code of type cudaError_t. You should always check them:

cudaError_t err = cudaMalloc(&d_a, bytes);
if (err != cudaSuccess) {
    std::cerr << "cudaMalloc failed: " << cudaGetErrorString(err) << std::endl;
    exit(EXIT_FAILURE);
}

For kernel launches, which are asynchronous, you can check for launch errors with cudaGetLastError() and wait for completion with cudaDeviceSynchronize() to catch execution errors:

kernel<<<grid, block>>>(args);
err = cudaGetLastError();   // catch launch errors
if (err != cudaSuccess) { ... }
cudaDeviceSynchronize();    // wait and catch execution errors
err = cudaGetLastError();
if (err != cudaSuccess) { ... }

Better yet, wrap calls in a macro:

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            std::cerr << "CUDA error in " << __FILE__ << ":" << __LINE__ << ": " \
                      << cudaGetErrorString(err) << std::endl; \
            exit(EXIT_FAILURE); \
        } \
    } while(0)

// Usage:
CUDA_CHECK(cudaMalloc(&d_a, bytes));

1.8 Profiling and Debugging Basics

Profiling: NVIDIA provides nvprof (deprecated but still useful) and the newer Nsight Systems / Nsight Compute. To profile a simple app:

  nvprof ./vecAdd

This shows kernel execution time, memory transfer times, etc. For deeper analysis, use Nsight Compute: ncu ./vecAdd.

Debugging: Use cuda-gdb (Linux) or the Visual Studio debugger on Windows. Compile with -g -G to generate debug info:

  nvcc -g -G -o vecAdd vecAdd.cu
  cuda-gdb ./vecAdd

Inside cuda-gdb, you can set breakpoints, inspect threads, and step through device code.

Memory error checking: compute-sanitizer (formerly cuda-memcheck) helps detect out‑of‑bounds accesses, race conditions, etc.

  compute-sanitizer ./vecAdd

1.9 Practice: Matrix Multiplication

Let's solidify our understanding with a classic: matrix multiplication C = A * B (N x N). We'll start with a naive kernel.

Naive implementation (each thread computes one element of C):

__global__ void matMulNaive(const float *A, const float *B, float *C, int N) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < N && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < N; k++) {
            sum += A[row * N + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

Launch configuration (assuming 16x16 blocks):

dim3 block(16, 16);
dim3 grid((N + block.x - 1) / block.x, (N + block.y - 1) / block.y);
matMulNaive<<<grid, block>>>(d_A, d_B, d_C, N);

Performance issues: This naive version has poor memory access patterns (especially for B, where each thread accesses a different column, causing non‑coalesced reads) and no data reuse. We'll improve it later with shared memory.

Part 2: Intermediate CUDA

Now that you're comfortable with basic CUDA, let's explore techniques to write faster kernels.

2.1 Performance Considerations

Occupancy: The ratio of active warps to maximum possible warps on an SM. Higher occupancy can hide latency but isn't always the goal. Use the CUDA Occupancy Calculator or the cudaOccupancyMaxPotentialBlockSize function to tune block size.

Warp divergence: Threads in a warp execute in lockstep. If they take different paths (e.g., an if‑else), the warp serially executes both paths, reducing performance. Try to keep control flow uniform within warps.

Coalesced memory access: Global memory bandwidth is best utilized when consecutive threads access consecutive memory locations. In our vector addition, thread i accesses a[i] – perfectly coalesced. In naive matrix multiplication, row‑major access of A is coalesced, but column‑major access of B is not. We'll fix that with shared memory.

2.2 Shared Memory and Bank Conflicts

Shared memory is divided into banks. If multiple threads access the same bank (with some exceptions), we get a bank conflict and the accesses are serialized. To avoid conflicts, ensure that threads within a warp access different banks (or all access the same word in broadcast mode). This is an advanced topic, but be aware of it when designing shared memory layouts.

2.3 Tiled Matrix Multiplication with Shared Memory

Let's rewrite matrix multiplication to use shared memory. The idea: load tiles of A and B into shared memory, then compute partial sums.

#define TILE_SIZE 16

__global__ void matMulShared(const float *A, const float *B, float *C, int N) {
    __shared__ float As[TILE_SIZE][TILE_SIZE];
    __shared__ float Bs[TILE_SIZE][TILE_SIZE];

    int bx = blockIdx.x, by = blockIdx.y;
    int tx = threadIdx.x, ty = threadIdx.y;

    // Identify the row and column of the C element to compute
    int row = by * TILE_SIZE + ty;
    int col = bx * TILE_SIZE + tx;

    float sum = 0.0f;

    // Loop over all tiles
    for (int tile = 0; tile < N / TILE_SIZE; ++tile) {
        // Load one tile of A and B into shared memory
        As[ty][tx] = A[row * N + (tile * TILE_SIZE + tx)];
        Bs[ty][tx] = B[(tile * TILE_SIZE + ty) * N + col];
        __syncthreads();  // Ensure tile is loaded before computation

        // Compute partial product for this tile
        for (int k = 0; k < TILE_SIZE; ++k) {
            sum += As[ty][k] * Bs[k][tx];
        }
        __syncthreads();  // Avoid overwriting shared memory before next tile
    }

    if (row < N && col < N) {
        C[row * N + col] = sum;
    }
}

What changed?

Each block loads a tile of A and a tile of B into shared memory.
All threads in the block cooperate to load the tiles (coalesced because threads in a warp load consecutive elements from global memory).
Computation uses the fast shared memory.
__syncthreads() ensures that all threads have finished loading before any thread starts using the data, and that no thread starts the next tile before all have finished the current tile.

This kernel can achieve near‑peak performance for large matrices.

2.4 Atomic Operations

When multiple threads need to update the same location (e.g., histogram, reduction), you can use atomic operations to avoid race conditions.

__global__ void histogram(const unsigned char *data, unsigned int *bins, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        atomicAdd(&bins[data[i]], 1);
    }
}

CUDA provides atomicAdd, atomicSub, atomicExch, atomicMin, atomicMax, atomicAnd, atomicOr, atomicXor for various data types.

Performance: Atomics can be slow if many threads contend for the same location. Use shared memory to reduce contention when possible (e.g., each block computes its own partial histogram, then combine with atomics).

2.5 Streams and Concurrency

Streams allow you to overlap kernel execution with data transfers and other operations. A stream is a sequence of operations that execute in order. Operations from different streams may overlap.

Example: overlapping copy and compute:

cudaStream_t stream1, stream2;
cudaStreamCreate(&stream1);
cudaStreamCreate(&stream2);

// Split data into two halves
int half = n / 2;

// Asynchronous memory copies (requires pinned host memory)
cudaMemcpyAsync(d_a1, h_a1, half_bytes, cudaMemcpyHostToDevice, stream1);
cudaMemcpyAsync(d_a2, h_a2, half_bytes, cudaMemcpyHostToDevice, stream2);

// Launch kernels in respective streams
kernel<<<grid, block, 0, stream1>>>(d_a1, d_b1, d_c1, half);
kernel<<<grid, block, 0, stream2>>>(d_a2, d_b2, d_c2, half);

// Copy back
cudaMemcpyAsync(h_c1, d_c1, half_bytes, cudaMemcpyDeviceToHost, stream1);
cudaMemcpyAsync(h_c2, d_c2, half_bytes, cudaMemcpyDeviceToHost, stream2);

cudaStreamSynchronize(stream1);
cudaStreamSynchronize(stream2);

cudaStreamDestroy(stream1);
cudaStreamDestroy(stream2);

Note: Use pinned memory (cudaMallocHost) for asynchronous transfers. Also, kernels from different streams may run concurrently if resources permit.

2.6 Events for Timing

CUDA events provide a way to measure elapsed time on the device accurately.

cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);

cudaEventRecord(start);
kernel<<<grid, block>>>(args);
cudaEventRecord(stop);

cudaEventSynchronize(stop);  // Wait for the stop event to be recorded
float milliseconds = 0;
cudaEventElapsedTime(&milliseconds, start, stop);
std::cout << "Kernel time: " << milliseconds << " ms\n";

cudaEventDestroy(start);
cudaEventDestroy(stop);

2.7 Unified Memory

Unified Memory (UM) creates a single pointer accessible from both CPU and GPU. Data is automatically migrated on demand. It simplifies programming but may introduce overhead.

Allocate UM:

float *x;
cudaMallocManaged(&x, N * sizeof(float));

You can then use x on both host and device without explicit cudaMemcpy. To ensure data is ready on the device, you can prefetch:

cudaMemPrefetchAsync(x, N * sizeof(float), deviceId, 0);

UM is great for prototyping and for algorithms with irregular data access. For performance‑critical code, manual transfers are still recommended.

2.8 Multi‑GPU Programming

Modern systems may have multiple GPUs. You can use them to scale your application.

Use cudaGetDeviceCount to find number of devices.
Set current device with cudaSetDevice(deviceId) before allocating memory or launching kernels.
Each device has its own memory space; you must manage data distribution and communication between devices (e.g., via PCIe or NVLink). Use peer‑to‑peer access if available (cudaDeviceCanAccessPeer, cudaDeviceEnablePeerAccess).

Simple multi‑GPU example:

int numDevices;
cudaGetDeviceCount(&numDevices);

#pragma omp parallel for
for (int dev = 0; dev < numDevices; ++dev) {
    cudaSetDevice(dev);
    // Allocate memory and launch kernel on this device
    // Copy data to/from host (strided or partitioned)
}

2.9 Advanced Example: Parallel Reduction

Reduction (summing an array) is a common pattern. We'll implement an optimized version that uses shared memory and avoids warp divergence.

__global__ void reduce(float *g_idata, float *g_odata, int n) {
    extern __shared__ float sdata[];  // Dynamic shared memory

    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;

    // Load data into shared memory (with bounds check)
    sdata[tid] = (i < n) ? g_idata[i] : 0.0f;
    __syncthreads();

    // Perform reduction in shared memory
    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }

    // Write result for this block to global memory
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}

Launch:

int threads = 256;
int blocks = (n + threads - 1) / threads;
reduce<<<blocks, threads, threads * sizeof(float)>>>(d_in, d_out, n);

After this kernel, you have partial sums per block. You may need a second kernel to combine them, or run a single block on the CPU.

Part 3: Introduction to HIP

3.1 What is HIP? Why Portability?

HIP (Heterogeneous‑compute Interface for Portability) is an AMD initiative to provide a C++ runtime API that can target both NVIDIA and AMD GPUs. It is essentially a thin abstraction layer: HIP code compiles to CUDA (using a CUDA backend) when targeting NVIDIA, and to AMD’s ROCm when targeting AMD. This means you can write your kernel once and run it everywhere.

Major benefits:

Vendor independence: No lock‑in to NVIDIA.
Single source: One codebase for multiple platforms.
Ease of porting: CUDA code can be mechanically converted using hipify tools.

3.2 Setting Up HIP

On AMD systems with ROCm:

Follow the ROCm installation guide. On Ubuntu:

  sudo apt update
  sudo apt install rocm-dev

Verify with hipcc --version.

On NVIDIA systems (HIP on CUDA):

Install the CUDA Toolkit first.
Then install HIP from ROCm’s GitHub or via package manager. Many Linux distributions provide a hip-cuda package.
After installation, hipcc will use the CUDA backend.

Environment:

Ensure hipcc is in your PATH.
For AMD, you may need to set ROCM_PATH if not default.

3.3 HIP API Overview

HIP intentionally mimics CUDA. Most CUDA functions have direct HIP equivalents with the hip prefix:

CUDA	HIP
`cudaMalloc`	`hipMalloc`
`cudaFree`	`hipFree`
`cudaMemcpy`	`hipMemcpy`
`cudaMemcpyAsync`	`hipMemcpyAsync`
`cudaDeviceSynchronize`	`hipDeviceSynchronize`
`cudaStreamCreate`	`hipStreamCreate`
`cudaEventCreate`	`hipEventCreate`
`cudaGetLastError`	`hipGetLastError`
`cudaError_t`	`hipError_t`
`cudaSuccess`	`hipSuccess`

The kernel language is essentially the same: __global__, __device__, __shared__, threadIdx, blockIdx, etc., are all supported.

The only major syntactic difference is the kernel launch. Instead of triple chevrons, HIP uses a function macro hipLaunchKernelGGL:

hipLaunchKernelGGL(kernel_name, dim3(grid), dim3(block), sharedMem, stream, args...);

If you prefer the <<<>>> syntax, you can define a macro (but it's not provided by default).

3.4 Porting a CUDA Program to HIP

Let’s port our vector addition from CUDA to HIP.

Original CUDA kernel (unchanged):

__global__ void vecAdd(const float *a, const float *b, float *c, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        c[i] = a[i] + b[i];
    }
}

Host code converted to HIP:

#include <iostream>
#include <hip/hip_runtime.h>

int main() {
    int n = 1 << 20;
    size_t bytes = n * sizeof(float);

    float *h_a = new float[n];
    float *h_b = new float[n];
    float *h_c = new float[n];

    for (int i = 0; i < n; i++) {
        h_a[i] = i * 1.0f;
        h_b[i] = (i + 1) * 1.0f;
    }

    float *d_a, *d_b, *d_c;
    hipMalloc(&d_a, bytes);
    hipMalloc(&d_b, bytes);
    hipMalloc(&d_c, bytes);

    hipMemcpy(d_a, h_a, bytes, hipMemcpyHostToDevice);
    hipMemcpy(d_b, h_b, bytes, hipMemcpyHostToDevice);

    int threadsPerBlock = 256;
    int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;

    // HIP kernel launch
    hipLaunchKernelGGL(vecAdd, dim3(blocksPerGrid), dim3(threadsPerBlock), 0, 0, d_a, d_b, d_c, n);

    hipMemcpy(h_c, d_c, bytes, hipMemcpyDeviceToHost);

    // Verification (omitted for brevity)

    hipFree(d_a); hipFree(d_b); hipFree(d_c);
    delete[] h_a; delete[] h_b; delete[] h_c;

    return 0;
}

Compile with:

hipcc -o vecAdd vecAdd.cpp   # on AMD
# or on NVIDIA
hipcc --offload-arch=sm_70 -o vecAdd vecAdd.cpp

Notice that the kernel itself is identical to the CUDA version – that’s the beauty of HIP!

3.5 HIP Kernel Launch Syntax Explained

hipLaunchKernelGGL takes the following arguments:

Kernel name (symbol)
Grid dimensions (dim3)
Block dimensions (dim3)
Dynamic shared memory size per block (size_t, default 0)
Stream (hipStream_t, default 0 for null stream)
Kernel arguments (comma‑separated)

If you have many arguments, it can become verbose. Some projects define a macro like:

#define HIP_LAUNCH(kernel, grid, block, shared, stream, ...) \
    hipLaunchKernelGGL(kernel, grid, block, shared, stream, __VA_ARGS__)

3.6 Error Handling in HIP

HIP error handling is identical to CUDA:

hipError_t err = hipMalloc(&d_a, bytes);
if (err != hipSuccess) {
    std::cerr << "HIP error: " << hipGetErrorString(err) << std::endl;
    exit(EXIT_FAILURE);
}

After a kernel launch, you can check for launch errors with hipGetLastError() and for execution errors after synchronization with hipDeviceSynchronize().

3.7 Debugging and Profiling HIP

Debugging:

On AMD, use rocgdb (similar to cuda‑gdb). Compile with -g and run rocgdb ./vecAdd.
On NVIDIA, you can still use cuda-gdb because the HIP code is translated to CUDA under the hood.

Profiling:

On AMD: Use rocprof to collect performance counters and traces.

  rocprof --hip-trace ./vecAdd

This generates a results.json file that can be viewed with tools like perfetto.

On NVIDIA: Use nvprof or nsight-compute as usual.

Memory checking:

On AMD, use compute-sanitizer (if available) or ROCm’s hip-memcheck.

3.8 HIP Examples: Matrix Multiplication

Our tiled matrix multiplication kernel works unchanged in HIP. Simply change the header and launch syntax. Let's see:

#include <hip/hip_runtime.h>
#include <iostream>

#define TILE_SIZE 16

__global__ void matMulShared(const float *A, const float *B, float *C, int N) {
    __shared__ float As[TILE_SIZE][TILE_SIZE];
    __shared__ float Bs[TILE_SIZE][TILE_SIZE];

    int bx = blockIdx.x, by = blockIdx.y;
    int tx = threadIdx.x, ty = threadIdx.y;
    int row = by * TILE_SIZE + ty;
    int col = bx * TILE_SIZE + tx;

    float sum = 0.0f;
    for (int tile = 0; tile < N / TILE_SIZE; ++tile) {
        As[ty][tx] = A[row * N + (tile * TILE_SIZE + tx)];
        Bs[ty][tx] = B[(tile * TILE_SIZE + ty) * N + col];
        __syncthreads();

        for (int k = 0; k < TILE_SIZE; ++k) {
            sum += As[ty][k] * Bs[k][tx];
        }
        __syncthreads();
    }

    if (row < N && col < N) {
        C[row * N + col] = sum;
    }
}

int main() {
    int N = 1024;
    size_t bytes = N * N * sizeof(float);

    // ... host allocations and initialization ...

    float *d_A, *d_B, *d_C;
    hipMalloc(&d_A, bytes);
    hipMalloc(&d_B, bytes);
    hipMalloc(&d_C, bytes);

    hipMemcpy(d_A, h_A, bytes, hipMemcpyHostToDevice);
    hipMemcpy(d_B, h_B, bytes, hipMemcpyHostToDevice);

    dim3 block(TILE_SIZE, TILE_SIZE);
    dim3 grid((N + block.x - 1) / block.x, (N + block.y - 1) / block.y);

    hipLaunchKernelGGL(matMulShared, grid, block, 0, 0, d_A, d_B, d_C, N);

    hipMemcpy(h_C, d_C, bytes, hipMemcpyDeviceToHost);

    // ... verification and cleanup ...

    return 0;
}

Compile with hipcc matmul.cpp -o matmul.

Part 4: Advanced HIP and Porting

4.1 Subtle Differences Between CUDA and HIP

While HIP aims for source compatibility, there are some differences to be aware of:

Warp vs Wavefront: NVIDIA warps are 32 threads; AMD wavefronts are 64 (on current architectures). Code that assumes a warp size of 32 (e.g., warp shuffle intrinsics) may need adjustment. HIP provides __hip_warp_size which returns the actual warp/wavefront size (32 on NVIDIA, 64 on AMD). Use it to write portable code.
Intrinsics: Some CUDA intrinsics (e.g., __popc, __ballot_sync) have HIP equivalents (__popc, __ballot_sync work, but check the HIP documentation for availability). For warp‑level primitives, prefer the portable __sync_warp and __all, __any.
Math functions: Most math functions (sin, exp, etc.) are identical. However, double‑precision performance may vary between vendors.
Architecture‑specific features: Features like tensor cores (NVIDIA) or matrix cores (AMD) are not portable. If you need them, you may have to write vendor‑specific code guarded by macros.

4.2 Using hipify Tools

AMD provides tools to automatically convert CUDA source to HIP:

hipify-perl: A Perl script that performs text‑based replacement.

  hipify-perl vecAdd.cu > vecAdd.cpp

It handles most API renames and some launch syntax, but you'll likely need to manually adjust the kernel launches (the triple chevrons are converted to hipLaunchKernelGGL).

hipify-clang: A more robust tool based on Clang’s AST. It can handle complex code and provides better accuracy.

  hipify-clang vecAdd.cu --cuda-path=/usr/local/cuda

It outputs HIP code to stdout or a file.

After conversion, review the output, especially for any use of CUDA‑specific features that don’t have direct HIP equivalents.

4.3 Writing Portable Code: Best Practices

To ensure your code runs efficiently on both platforms:

Use hip runtime API for memory management, streams, events.
Avoid hard‑coding warp size: Use __hip_warp_size or compute dynamically.
Use portable intrinsics: Prefer __syncthreads, __threadfence, __syncwarp (though __syncwarp may need a mask; HIP provides __syncwarp(mask) with similar semantics).
Conditional compilation for vendor‑specific optimizations:

  #ifdef __HIP_PLATFORM_NVIDIA__
      // NVIDIA-specific code (e.g., using CUDA tensor cores)
  #elif __HIP_PLATFORM_AMD__
      // AMD-specific code
  #endif

Profile on both platforms to identify bottlenecks. What's fast on NVIDIA may not be on AMD due to different cache hierarchies, wavefront size, etc. Tune block sizes and tile sizes accordingly.

4.4 Performance Optimization on AMD GPUs

When targeting AMD, use ROCm tools:

rocprof: For profiling.
rocm-smi: To monitor GPU utilization, clock speeds, temperature.
OmniTrace: For tracing API calls and kernel execution.

Key optimization considerations for AMD:

Wavefront size: 64 threads. Your block size should be a multiple of 64 to avoid idle threads.
Shared memory bank conflicts: AMD’s shared memory has 32 banks (like NVIDIA) but the mapping may differ. Experiment with padding to reduce conflicts.
Memory coalescing: Still crucial; consecutive threads should access consecutive memory.
Occupancy: Similar concepts apply, but the register file and shared memory sizes differ. Use hipOccupancyMaxPotentialBlockSize to get recommendations.

4.5 Real-world Example: Sobel Edge Detection (Portable)

Let's combine everything into a real‑world kernel: Sobel edge detection on an 8‑bit grayscale image. We'll write it portably for HIP/CUDA.

#include <hip/hip_runtime.h>  // or <cuda_runtime.h> for CUDA

#define TILE_SIZE 16
#define MASK_SIZE 3

__global__ void sobel(const unsigned char *in, unsigned char *out, int width, int height) {
    // Shared memory for tile + halo (1 pixel halo on each side)
    __shared__ unsigned char tile[TILE_SIZE + 2][TILE_SIZE + 2];

    int x = blockIdx.x * TILE_SIZE + threadIdx.x;
    int y = blockIdx.y * TILE_SIZE + threadIdx.y;

    // Global coordinates including halo
    int gx = x - 1;
    int gy = y - 1;

    // Load tile with boundary check
    if (gx >= 0 && gx < width && gy >= 0 && gy < height)
        tile[threadIdx.y + 1][threadIdx.x + 1] = in[gy * width + gx];
    else
        tile[threadIdx.y + 1][threadIdx.x + 1] = 0;

    // Load halo (simplified: we rely on neighboring threads)
    // For a full implementation, you'd also load the border regions.
    __syncthreads();

    // Compute Sobel only for inner tile pixels (not halo)
    if (threadIdx.x < TILE_SIZE && threadIdx.y < TILE_SIZE && x < width && y < height) {
        int Gx = -tile[threadIdx.y][threadIdx.x] + tile[threadIdx.y][threadIdx.x + 2]
                 -2*tile[threadIdx.y+1][threadIdx.x] + 2*tile[threadIdx.y+1][threadIdx.x+2]
                 -tile[threadIdx.y+2][threadIdx.x] + tile[threadIdx.y+2][threadIdx.x+2];
        int Gy = -tile[threadIdx.y][threadIdx.x] -2*tile[threadIdx.y][threadIdx.x+1] - tile[threadIdx.y][threadIdx.x+2]
                 +tile[threadIdx.y+2][threadIdx.x] +2*tile[threadIdx.y+2][threadIdx.x+1] + tile[threadIdx.y+2][threadIdx.x+2];
        int mag = min(255, (int)hypot(Gx, Gy));
        out[y * width + x] = mag;
    }
}

Launch:

dim3 block(TILE_SIZE, TILE_SIZE);
dim3 grid((width + block.x - 1) / block.x, (height + block.y - 1) / block.y);
sobel<<<grid, block>>>(d_in, d_out, width, height); // CUDA
// HIP: hipLaunchKernelGGL(sobel, grid, block, 0, 0, d_in, d_out, width, height);

This kernel uses shared memory to reduce global memory traffic. Each block loads a tile plus halo, then computes Sobel for its interior. It's a great example of a real‑world, portable GPU kernel.

Part 5: Conclusion and Next Steps

Congratulations! You've journeyed from the absolute basics of GPU programming to writing advanced, portable kernels that run on both NVIDIA and AMD hardware. You've learned:

The CUDA programming model, memory hierarchy, and essential APIs.
How to write and optimize kernels for maximum performance.
The HIP ecosystem and how to write truly portable code.
Debugging and profiling techniques for both platforms.

Where to Go from Here

Explore official documentation:
- CUDA C++ Programming Guide
- HIP Documentation
Study advanced examples: Reduction, scan, stencils, FFTs.
Contribute to open‑source projects that use CUDA/HIP (e.g., PyTorch, TensorFlow, ROCm libraries).
Experiment with multi‑GPU and distributed computing (e.g., using MPI with GPUs).
Learn about new hardware features like tensor cores (NVIDIA) and matrix cores (AMD) – but remember to keep your code portable or use conditional compilation.

GPU programming is a powerful skill that opens doors to high‑performance computing, machine learning, and graphics. With CUDA and HIP under your belt, you're well‑equipped to tackle the most demanding computational challenges. Now go forth and write some fast, portable code!

If you found this guide helpful, please share it with your fellow developers. Have questions or want to share your own experiences? Drop a comment below – I’d love to hear from you. Happy coding!

OLSRT: The Runtime I Built on Windows 7 When I Couldn't Run MSVC (And Why It Might Be the Last Runtime You'll Ever Need)

Javad — Sat, 14 Feb 2026 16:04:54 +0000

Hey Dev Community!

A while back, I was sitting on my old Windows 7 machine, trying to build a PHP extension. I wanted to bring modern async I/O, WebSocket support, and an event loop to PHP — something like what Node.js has, but for PHP. The problem? phpize requires MSVC, and on Windows 7 with MinGW, I couldn't make it work. No matter what I tried, MSVC just wouldn't install.

So I sat there, frustrated, staring at my screen. And then a thought hit me:

"C can talk to any language. So why am I limiting myself to just one?"

That question changed everything. Instead of building a PHP extension, I started building something much bigger: a language-agnostic runtime that any language could talk to. A runtime that would give every language access to async I/O, WebSocket, event loops, and more.

That was the birth of OLSRT.

The first version (0.1) was a mess. Just a single src folder, no header files, and barely any structure. But it worked — enough to prove the concept.

Then I rewrote everything from scratch. Today, OLSRT is at v1.2 (Atom) , and it's production-ready. Here's what's in it right now:

✅ What OLSRT v1.2 Offers Today

1. Actor Model with Full Process Isolation

Each actor runs in its own isolated process with a private memory arena. If one actor crashes, the rest keep running — just like Erlang/OTP, but in C. The implementation includes:

Lock‑free mailboxes with zero‑copy message passing
Supervisors with restart strategies (permanent, transient, temporary)
Process linking and monitoring (Erlang‑style)

2. Green Threads (Not Just Coroutines)

OLSRT has a full green thread implementation with:

Work stealing across multiple CPU cores
NUMA‑aware memory allocation
SIMD‑accelerated memory operations (memcpy, memset, memcmp)
Guard pages to detect stack overflows
Context switching in assembly for x86_64, ARM64, ARM, and RISC‑V

3. High‑Performance Memory Management

Arena‑based memory allocators (2MB default, expandable)
Zero‑copy IPC between processes
Lock‑free data structures for hot paths
Cache‑line padding to prevent false sharing

4. Async I/O and Event Loop

Cross‑platform event loop with epoll (Linux), kqueue (macOS/BSD), and select (fallback)
Timer support with deadlines
I/O event registration with callbacks

5. Channels and Streams

Thread‑safe FIFO channels with bounded/unbounded capacity
Reactive streams with operators (map, filter, take, merge, debounce)
Backpressure support with demand‑based flow control

6. Promises and Futures

Thread‑safe promise/future implementation
Continuations (.then())
Integration with the event loop for waking

7. Synchronization Primitives

Mutexes, condition variables, read‑write locks, semaphores
All implemented natively on Windows, Linux, macOS, and BSD

8. Networking (TCP/UDP)

Non‑blocking TCP and UDP sockets
Full integration with the event loop
Promise‑based async API (send/recv/accept/connect)

🚧 What's Coming in v1.3 (In Development Now)

Based on our latest commits (as of February 2026), v1.3 will include:

Hot‑Coding References (HCR) – Never pause the main with very lightweight branches.
OLSRT Coroutines (ORoutines) – Even lighter‑weight concurrency
Fully documented runtime
Optimized build systems (CMake and Makefile improvements)

Why OLSRT Exists (The Philosophy)

I didn't set out to build a universal runtime. I just wanted to make PHP async. But along the way, I realized something:

The future of software is polyglot. We use PHP for the web, Python for data, Rust for systems, and JavaScript for the frontend. But today, each language has its own runtime, its own concurrency model, and its own set of limitations. OLSRT is an attempt to unify them — to give every language access to the same powerful primitives: actors, green threads, zero‑copy IPC, and fault‑tolerant supervisors.

Where We Stand Today

OLSRT v1.2 is stable and production‑ready on Linux and BSD (Windows support is coming in v1.3). It has:

118 commits
Full Doxygen documentation (for v1.0, v1.2 docs are in progress)
A growing set of bindings in the works (PHP, JavaScript, Python, Rust, Go, .NET, Dart, Java, Kotlin, D, Elixir, Erlang, Ruby…)
A small but passionate team of 1 core developer ( And it's me :) )

And yes, right now we have exactly 3 GitHub stars. That's it.

Why I'm Writing This

Because OLSRT needs you. Not just as a user, but as a contributor.

This project has the potential to change how we build software — but it can't do that alone. We need:

Developers to star the repo (it sounds small, but it makes a huge difference in visibility)
People to open issues (every bug report, every feature request, every question)
Contributors to send pull requests — even if it's just fixing a typo in the docs
Community members to spread the word (share this post, tweet about OLSRT, talk about it at meetups)

How You Can Help Right Now

Visit the repository: github.com/OverLab-Group/OLSRT
Star it. One click. It takes 2 seconds.
Read the code. It's well‑documented (I promise).
Open an issue. Found a bug? Have an idea? Let me know.
Submit a PR. Even a tiny one counts. Fix a typo, improve a comment, add a test.
Tell a friend. Share this post on Twitter, LinkedIn, or your local dev community.

The Road Ahead

We have a roadmap to v12. But we won't get there alone. OLSRT is open source, and it belongs to the community. With your help, it can become the universal runtime that powers the next generation of software — from real‑time chat servers to telecom infrastructure to embedded IoT devices.

It all started on a Windows 7 machine with a broken MSVC install. Now it's up to you to help it grow.

Join us on GitHub

🚀 OLSRT v1.2: A Powerful Runtime for All Programming Languages!

Javad — Sun, 08 Feb 2026 23:10:18 +0000

Hey Dev Community!

I'm thrilled to be writing blogs for you again after a while!

📖 Introduction

OLSRT (OverLab Streams Runtime) is a powerful, minimal runtime (~90 KB) delivered as shared libraries (.so for Linux/BSD) and dynamic libraries (.dll for NT), written in C, and designed to power all programming languages!

If your language can talk to C (yes, even COBOL!), you can use OLSRT directly via FFI.

✅ Supported Languages: All of them! (JavaScript (with nodeffi), PHP (with C extensions), Python (with cffi), Rust, D, C++, Go, Dart, Kotlin, Java, C#, Julia, and any language you use!)

📌 Note for C, C++, D, Rust, and other languages with C at their core:

You can use OLSRT directly in your programs—even without any .so/.dll files! Just include the header for your required module and use OLSRT with zero bridging.

🔥 What’s New in Version 1.2?

This release brings exciting and advanced features to OLSRT:

🎭 Improved Actors & Supervisors

The first C-based Actors designed for all programming languages! OLSRT now delivers an advanced Actor model for concurrent systems.

⚡ Enhanced Async I/O, Await, Promises & Futures

Full support for modern concurrency patterns across all languages.

🧵 Green Threads & Modern Fibers

Now you can leverage OLSRT's powerful Green Threads in all your projects.

✨ Many Other Improvements

Check out the full release notes for OLSRT Atom v1.2.

📘 Documentation for v1.2 is still in progress, but you can refer to the v1.0 docs (the core concepts remain similar). We’re actively working on examples and comprehensive documentation.

🚧 In Development: OLSRT Version 1.3

We’re already working on v1.3 with revolutionary features:

🧩 Fully Isolated Actors

💥 “Let It Crash!” Philosophy

🔄 ORoutines (OLSRT’s Modern Coroutines)

🧠 HCR (Hot-Coding References) – A Brand New Paradigm!

🤔 What is HCR? A New Paradigm!

OverLab Group divides concurrency paradigms into two categories:

🕐 Waiting Paradigms	⚡ Non-Waiting Paradigms
Blocking:	Branch-based:
– Sync	– HCR (First of its kind!)
– Lock
– Semaphore
Non-Blocking:
– Async
– Await
– Actors

🧪 HCR in Action:

In waiting paradigms, you write code like this:

int x = a + b;

This pauses the main program for milliseconds.

But with HCR:

ol_caller* x = ol_reference_do(a + b);

Here, OLSRT creates a lightweight branch (~1 KB) for x, a, and b, processes them without stopping the main program, and merges the result back seamlessly!

📢 Call to Action

If you like the OLSRT project and want to give us some energy, even doing one of the following helps a ton:

⭐ Star the repository (One energy)
🔀 Fork the repository (One super energy)
💬 Comment, share your ideas, and follow us (One mega energy)
🗣️ Discuss with us on GitHub (Three mega energies)
🛠️ Clone and use the project (One ultra energy)
🐛 Open an issue with your ideas (Three ultra energies)
🤝 Contribute to the repository (One mega ultra super energy)

🔗 GitHub Repository:

https://github.com/OverLab-Group/OLSRT

We welcome all forms of feedback and contribution!

Hope you enjoyed this update.

Have an ORoutine day! 🚀