Tiling Matrix Multiplication in OpenCL with __local Memory

The naive kernel is memory‑bound

A straightforward OpenCL matrix‑multiply kernel that computes C = A × B usually looks like this:

__kernel void matmul_naive(__global const float *A,
                           __global const float *B,
                           __global float *C,
                           const int M, const int N, const int K)
{
    int row = get_global_id(0);
    int col = get_global_id(1);
    if (row >= M || col >= N) return;
    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;
}

Each work‑item reads one element of A and one element of B for every iteration of the inner loop. For an output element we perform N multiply‑adds but 2N global‑memory loads, giving a low arithmetic intensity. On GPUs the compute units spend most of their time waiting for data, so the kernel is bandwidth‑bound.

Work‑groups and __local scratchpad

OpenCL divides the N‑dimensional index space into work‑groups. All work‑items in a group execute on the same compute unit and share a programmer‑managed region called __local memory. This memory lives on‑chip (often tens of kilobytes) and can be accessed much faster than global memory.

To reuse data we cooperatively load a tile of A and a tile of B into __local arrays, synchronize with a barrier, then each work‑item multiplies the loaded tiles without touching global memory again. After the computation we place another barrier before overwriting the tiles for the next iteration.

The two barriers are essential:

• barrier(CLK_LOCAL_MEM_FENCE) ensures that all stores to local memory are visible to every work‑item in the group before any reads.
• A second barrier after the compute step guarantees that no work‑item starts loading the next tile while another is still using the current tile.

If a barrier appears inside divergent control flow (e.g., only some work‑items take a branch), the behavior is undefined.

A 16×16 tiled example

The following kernel assumes square tiles of size TILE = 16. Each work‑item computes one element of C. Host code must launch an NDRange that covers the output matrix and query device limits at runtime (see the checklist later).

#define TILE 16
__kernel void matmul_tiled(__global const float *A,
                           __global const float *B,
                           __global float *C,
                           const int M, const int N, const int K)
{
    int row = get_group_id(0) * TILE + get_local_id(0);
    int col = get_group_id(1) * TILE + get_local_id(1);

    __local float As[TILE][TILE];
    __local float Bs[TILE][TILE];

    float sum = 0.0f;
    for (int t = 0; t < (K + TILE - 1) / TILE; ++t) {
        int tiledK = t * TILE + get_local_id(0);
        int tiledK2 = t * TILE + get_local_id(1);

        // Load tile of A
        if (row < M && tiledK < K)
            As[get_local_id(1)][get_local_id(0)] = A[row * K + tiledK];
        else
            As[get_local_id(1)][get_local_id(0)] = 0.0f;

        // Load tile of B
        if (col < N && tiledK2 < K)
            Bs[get_local_id(1)][get_local_id(0)] = B[tiledK2 * N + col];
        else
            Bs[get_local_id(1)][get_local_id(0)] = 0.0f;

        barrier(CLK_LOCAL_MEM_FENCE);

        for (int k = 0; k < TILE; ++k)
            sum += As[get_local_id(1)][k] * Bs[k][get_local_id(0)];

        barrier(CLK_LOCAL_MEM_FENCE);
    }

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

Explanation:

• Each work‑item calculates its global row and col from the group ID and local ID.
• Two __local arrays hold the current tile of A and B.
• Inside the tiling loop we load elements into local memory, guard against out‑of‑bounds with zero‑fill, then synchronize.
• The inner product accumulates using the locally cached tile.
• A second barrier prevents race conditions when the next iteration overwrites the tile.

Trade‑offs and tuning checklist

Tiling improves data reuse but introduces synchronization overhead, extra register pressure, and boundary‑handling code. Larger tiles increase reuse per global load but reduce the number of work‑groups that can reside simultaneously on a compute unit, potentially lowering occupancy. The optimal tile size depends on:

• CL_DEVICE_LOCAL_MEM_SIZE – total scratchpad available.
• CL_KERNEL_LOCAL_MEM_SIZE – local memory required by the kernel (query with clGetKernelWorkGroupInfo).
• CL_DEVICE_MAX_WORK_GROUP_SIZE – limits how many work‑items can be in a group.

On CPUs the __local region is often mapped to cache; many OpenCL compilers already perform loop blocking, so explicit tiling may give little or no benefit. Always profile first.

Actionable checklist

1. Run clinfo (or query via the host program) to obtain CL_DEVICE_LOCAL_MEM_SIZE, CL_DEVICE_LOCAL_MEM_TYPE, and CL_DEVICE_MAX_WORK_GROUP_SIZE for your target device.
2. Compile both the naive and tiled kernels with -cl-std=CL1.2 (or the version you target).
3. Create a command queue with CL_QUEUE_PROFILING_ENABLE.
4. Enqueue the kernel, record the event, and after completion call clGetEventProfilingInfo for CL_PROFILING_COMMAND_START and CL_PROFILING_COMMAND_END to get elapsed time.
5. Validate correctness by copying C back to the host and comparing to a reference gemm implementation on small matrices (e.g., 64×64).
6. If the tiled version is not faster, experiment with different TILE values that satisfy the local‑memory constraints, or consider using vendor‑provided libraries (e.g., clBLAS).

Remember: never hard‑code tile sizes; query limits at runtime and handle matrices whose dimensions are not multiples of the tile size with the zero‑fill guards shown above.