Programación GPU: CUDA vs ROCm vs Vulkan Compute — Warp Scheduling, Shared Memory y Multiplicación de Matrices

El Modelo de Ejecución: Warps, Wavefronts y Subgroups

A continuación se describe qué ocurre cuando el código se ejecuta en la GPU.

NVIDIA (CUDA): Warps de 32

La unidad de ejecución fundamental en NVIDIA (CUDA) es el warp — 32 threads que ejecutan la misma instrucción en lockstep (SIMT: Single Instruction, Multiple Threads). Un Streaming Multiprocessor (SM) puede schedular múltiples warps, pero dentro de un solo warp, los 32 threads están en el mismo program counter.

// Esto parece inocente pero asesina el rendimiento
__global__ void bad_branch(float* data, int n) {
 int tid = threadIdx.x + blockIdx.x * blockDim.x;
 if (tid < n) {
 if (data[tid] > 0.0f) {
 data[tid] = sqrtf(data[tid]); // Camino A
 } else {
 data[tid] = -sqrtf(-data[tid]); // Camino B
 }
 }
}
// Si threads 0-15 toman Camino A y threads 16-31 toman Camino B,
// el warp ejecuta AMBOS caminos secuencialmente. 50% utilización.

AMD (ROCm/HIP): Wavefronts de 64

El equivalente de AMD es el wavefront — 64 threads en arquitecturas CDNA/RDNA. Esto no es solo una diferencia de nombre; cambia fundamentalmente cómo considerars sobre occupancy y primitise va a nivel warp.

// Código HIP — se ve idéntico a CUDA pero wavefront = 64
__global__ void hip_reduce(float* data, float* result, int n) {
 __shared__ float sdata[256];
 int tid = threadIdx.x;

 sdata[tid] = (tid < n) ? data[tid] : 0.0f;
 __syncthreads();

 // PELIGRO: esta reducción asume warp size de 32
 // En AMD, las últimas iteraciones 'unrolled' necesitan cubrir 64 threads
 for (int s = blockDim.x / 2; s > 32; s >>= 1) {
 if (tid < s) sdata[tid] += sdata[tid + s];
 __syncthreads();
 }

 // Este unrolling a nivel warp está MAL en AMD
 if (tid < 32) { // Debería ser < 64 en AMD!
 volatile float* vsmem = sdata;
 vsmem[tid] += vsmem[tid + 32];
 vsmem[tid] += vsmem[tid + 16];
 vsmem[tid] += vsmem[tid + 8];
 vsmem[tid] += vsmem[tid + 4];
 vsmem[tid] += vsmem[tid + 2];
 vsmem[tid] += vsmem[tid + 1];
 }
}

Este bug exacto aparece al portar codebases CUDA a ROCm para clusters MI250X. hipify-perl traduce la sintaxis perfecto pero no marca la suposición de warp size. Los resultados pueden estar mal por ~3% — lo suficientemente cerca como para pasar tests unitarios con tolerancias flojas pero completamente incorrecto para simulaciones de física.

Vulkan Compute: Subgroups

Vulkan no especifica un tamaño de subgroup. Puede ser 4, 8, 16, 32, 64 o incluso 128 dependiendo del hardware y driver. Lo se consulta en runtime:

#version 450
#extension GL_KHR_shader_subgroup_basic : enable
#extension GL_KHR_shader_subgroup_arithmetic : enable

layout(local_size_x = 256) in;
layout(set = 0, binding = 0) buffer Data { float data[]; };
layout(set = 0, binding = 1) buffer Result { float result; };

void main() {
 uint tid = gl_LocalInvocationID.x;

 // subgroupSize es un built-in — puede ser 32, 64, etc.
 float val = data[tid];
 float subgroup_sum = subgroupAdd(val);

 if (subgroupElect()) {
 atomicAdd(result, subgroup_sum);
 }
}

Shared Memory: El Recurso que Hace o Deshace

La shared memory (CUDA/HIP) es la SRAM rápida on-chip que hace la computación GPU práctica. En GPUs NVIDIA modernas, está organizada en 32 banks, cada uno de 4 bytes de ancho.

Bank Conflicts: El Asesino Silencioso

// Matriz almacenada en shared memory: 32 x 32 floats
__shared__ float tile[32][32];

// Thread (tx, ty) lee por columna:
// Thread 0 lee tile[0][0] -> bank 0
// Thread 1 lee tile[1][0] -> bank 0 (!!!)
// Thread 2 lee tile[2][0] -> bank 0 (!!!)
// ... los 32 threads pegan en bank 0 -> conflicto de 32 vías!

// Fix: paddear cada fila con 1 elemento
__shared__ float tile[32][33]; // <-- el número mágico

// Ahora:
// Thread 0 lee tile[0][0] -> bank 0
// Thread 1 lee tile[1][0] -> bank 1 (offset por 33 % 32 = 1)
// Thread 2 lee tile[2][0] -> bank 2
// ... los 32 threads pegan en banks distintos -> sin conflicto!

Ese [32][33] en vez de [32][32] es el tipo de cosa que convierte un kernel al 40% de pico en uno al 80%. Un elemento extra por fila. Este solo cambio puede dar speedups de 2x en kernels memory-bound.

Multiplicación de Matrices: SGEMM Optimizado

El triple loop naive es insuficiente para uso en producción. El siguiente kernel SGEMM tileado fue optimizado iterativamente. Esta versión consigue ~75% del rendimiento de cuBLAS en una A100:

#define BM 128
#define BN 128
#define BK 8
#define TM 8
#define TN 8

__global__ void sgemm_optimized(
 const float* __restrict__ A,
 const float* __restrict__ B,
 float* __restrict__ C,
 int M, int N, int K
) {
 const int bx = blockIdx.x;
 const int by = blockIdx.y;
 const int tx = threadIdx.x;
 const int ty = threadIdx.y;

 __shared__ float As[BK][BM + 1]; // Transpuesta para coalescing
 __shared__ float Bs[BK][BN + 1];

 float accum[TM][TN] = {0.0f};
 float a_reg[TM];
 float b_reg[TN];

 const int row_base = by * BM;
 const int col_base = bx * BN;

 for (int k = 0; k < K; k += BK) {
 // Carga colaborativa de tile A en shared memory
 #pragma unroll
 for (int i = 0; i < BM; i += blockDim.y) {
 int row = row_base + ty + i;
 int col = k + tx;
 if (row < M && col < K) {
 As[tx][ty + i] = A[row * K + col];
 } else {
 As[tx][ty + i] = 0.0f;
 }
 }

 #pragma unroll
 for (int i = 0; i < BN; i += blockDim.y) {
 int row = k + ty;
 int col = col_base + tx + i;
 if (row < K && col < N) {
 Bs[ty][tx + i] = B[row * N + col];
 } else {
 Bs[ty][tx + i] = 0.0f;
 }
 }

 __syncthreads();

 #pragma unroll
 for (int kk = 0; kk < BK; kk++) {
 #pragma unroll
 for (int m = 0; m < TM; m++) {
 a_reg[m] = As[kk][ty * TM + m];
 }
 #pragma unroll
 for (int n = 0; n < TN; n++) {
 b_reg[n] = Bs[kk][tx * TN + n];
 }
 #pragma unroll
 for (int m = 0; m < TM; m++) {
 #pragma unroll
 for (int n = 0; n < TN; n++) {
 accum[m][n] += a_reg[m] * b_reg[n];
 }
 }
 }

 __syncthreads();
 }

 #pragma unroll
 for (int m = 0; m < TM; m++) {
 #pragma unroll
 for (int n = 0; n < TN; n++) {
 int row = row_base + ty * TM + m;
 int col = col_base + tx * TN + n;
 if (row < M && col < N) {
 C[row * N + col] = accum[m][n];
 }
 }
 }
}

Decisiones clave de optimización:

1. Store transpuesto en shared memory para A — Cuando los threads leen una columna de As, acceden a filas consecutivas que mapean a banks distintos.
2. Register tiling (TM x TN) — Cada thread computa un sub-tile de 8x8, amortizando loads de shared memory sobre 64 operaciones FMA por paso de K.
3. Punteros __restrict__ — Le dice al compilador que A, B, C no se solapan en memoria.
4. #pragma unroll — Crítico para los inner loops.

El Port a ROCm: No es Solo Find-and-Replace

Al portar esto a HIP para MI250X, surgen varios problemas no obvios. La shared memory de AMD es 64 KB por CU vs los 164 KB configurables de NVIDIA en A100. Si tu estrategia de tiling necesita más de 64 KB, se necesita tiles más pequeños en AMD, lo que cambia toda tu estrategia de optimización.

Benchmarks: Números Reales en Hardware Real

Los siguientes benchmarks muestran SGEMM 4096x4096 en tres setups:

| Implementación | Hardware | TFLOPS | % del Pico | |---------------|----------|--------|-----------| | cuBLAS | A100 80GB | 17.8 | 91% | | Kernel CUDA custom | A100 80GB | 14.7 | 75% | | rocBLAS | MI250X | 15.2 | 88% | | Kernel HIP custom | MI250X | 11.9 | 69% | | Vulkan Compute (GLSL) | RTX 4090 | 9.3 | 56% | | Vulkan Compute (GLSL) | RX 7900 XTX | 7.1 | 52% |

Las bibliotecas del vendor (cuBLAS, rocBLAS) hacen trucos a nivel assembly que son impracticables de replicar en código de alto nivel — cosas como software pipelining de loads de shared memory y asignación explícita de registros.

Profiling: Distribución del Tiempo de Ejecución

La habilidad más valiosa en programación GPU es leer un output de profiler. Para el kernel CUDA:

Kernel: sgemm_optimized
Duration: 2.31 ms
Theoretical Occupancy: 75%
Achieved Occupancy: 71.2%
SM Throughput: 87.3%
Memory Throughput: 62.1% <-- No es memory bound
Compute Throughput: 87.3% <-- Es compute bound

L1/TEX Hit Rate: 94.7%
Shared Memory Bank Conflicts: 0 <-- El padding funcionó

Warp Stall Reasons:
 stall_mio_throttle: 23.1% <-- Ancho de banda de shared memory
 stall_math_pipe_throttle: 18.7% <-- Bien! La ALU está saturada
 stall_short_scoreboard: 15.2% <-- Esperando resultados MIO
 stall_not_selected: 12.4% <-- Scheduler no pudo elegir este warp

Ese stall_mio_throttle al 23.1% indica que el ancho de banda de shared memory es el próximo bottleneck — el kernel carga de shared memory más rápido de lo que el hardware puede servir. El fix sería aumentar el register tile size (TM/TN) para hacer más computación por load de shared memory, al costo de más presión de registros y potencialmente menor occupancy.

Este es el juego. Cada optimización crea un nuevo bottleneck. Terminaste cuando el profiler te muestra que se está limitado por algo fundamental — throughput de ALU, ancho de banda de memoria, o límites de occupancy por presión de registros.

Framework de Decisión

Una evaluación honesta de los tres ecosistemas tras uso extensivo:

CUDA tiene el mejor tooling por lejos. Si se está en hardware NVIDIA, no hay razón para usar otra cosa para compute.

ROCm mejoró dramáticamente pero todavía se siente 2 años atrás de CUDA en tooling. El lenguaje HIP en sí es 98% source-compatible con CUDA, lo que hace el porting directo. Vale la pena si teners hardware AMD.

Vulkan Compute es la opción correcta cuando se necesita una codebase corriendo en Intel, NVIDIA, AMD y hasta GPUs mobile. El trade-off es menos control sobre específicos del hardware y una API más verbosa.

Elegir basándote en tus restricciones, no en tribalismo. Y por favor, profileá antes de optimizar.