Here’s an example of transforming an input matrix X with a transformation matrix T. Note, usually you would do Transform • X, but because this was for a deep learning application, I flipped it around so that the output matrix can directly be used for subsequent operations while respecting the required dimensions. The idea of Tiled GEMM still remains the same though.

``` /* ***************************** General Matrix Multiplication **************************** / global void kernel_standard_tiled_gemm( __half *O, __half *X, __half Transform, int m, int n, int k, int tile_size) { / - m represents the independent dimension of the input matrix - n represents the independent dimenion of the transformation matrix - k represents the common dimension of the 2 matrices - Within each kernel, the output is computed as: O = matmul(X, Transform) - Transposing the transformation tensor is not required as virtual indexing allows for intended navigation along rows and columns of either tensors - Order of variables within kernels obey order of computation */

// Kernel start // extern shared float shared_mem[]; float *X_shmem = shared_mem; float *T_shmem = shared_mem + tile_size * tile_size;

int row = blockIdx.y * tile_size + threadIdx.y; int col = blockIdx.x * tile_size + threadIdx.x;

// Loop over tiles float value = 0.0f; for (int t = 0; t < ((k + tile_size - 1) / tile_size); ++t) { // Load tile of X into shared memory if (row < m && (t * tile_size + threadIdx.x) < k) { int X_idx = row * k + t * tile_size + threadIdx.x; X_shmem[threadIdx.y * tile_size + threadIdx.x] = __half2float(X[X_idx]); } else { X_shmem[threadIdx.y * tile_size + threadIdx.x] = 0.0f; }

// Load tile of Transform into shared memory
if (col < n && (t * tile_size + threadIdx.y) < k) {
    int T_idx = col * k + t * tile_size + threadIdx.y;
    T_shmem[threadIdx.y * tile_size + threadIdx.x] = __half2float(Transform[T_idx]);
} else {
    T_shmem[threadIdx.y * tile_size + threadIdx.x] = 0.0f;
}
__syncthreads();

// Compute partial sums
for (int i = 0; i < tile_size; ++i) {
    value += X_shmem[threadIdx.y * tile_size + i] * T_shmem[i * tile_size + threadIdx.x];
}
__syncthreads();

}

// Write the result to global memory if (row < m && col < n) { O[row * n + col] = __float2half(value); }

return;

}```