Coalescencia y patrón de acceso

Dos kernels que leen exactamente el mismo número de valores pueden diferir en un factor de treinta y dos en tiempo. La diferencia no está en cuántos bytes pides sino en qué direcciones los pides, porque la memoria no se sirve por valores sueltos: se sirve por transacciones de una línea entera. Si los treinta y dos carriles de un grupo SIMD piden direcciones dispersas, el bus mueve treinta y dos líneas completas y tira el noventa y siete por ciento de lo que ha movido. Y lo peor es que el código de las dos versiones se parece tanto que nadie sospecha.

La unidad real de lectura es la transacción

El controlador de memoria no sabe leer cuatro bytes. Su unidad de trabajo es la línea de caché, que según la arquitectura y el nivel de la jerarquía mide 32, 64 o 128 bytes; los sectores de DRAM más habituales son de 32 bytes y las líneas de L2 de 64 o 128. Pidas lo que pidas, se mueve una línea entera y se cobra una línea entera.

De ahí sale la aritmética de la coalescencia. Si los 32 carriles de un grupo SIMD leen 32 f32 consecutivos y alineados, eso son 128 bytes contiguos: una transacción, cien por cien de aprovechamiento. Si cada carril lee un f32 separado 256 bytes del anterior, cada uno cae en una línea distinta: 32 transacciones de 128 bytes de las que aprovechas 4. Has movido 4096 bytes para usar 128. Un 3,1 por ciento de eficiencia y treinta y dos veces más tráfico por el mismo resultado.

Fíjate en la segunda fila, que es la trampa silenciosa: un acceso perfectamente contiguo pero que empieza a mitad de una línea toca dos líneas en vez de una y desperdicia la mitad del bus. Por eso minStorageBufferOffsetAlignment y minUniformBufferOffsetAlignment valen 256 en los límites por defecto: la especificación te obliga a que las bases de tus vistas caigan en fronteras razonables, y con eso te ahorra el caso patológico. Lo que no puede hacer por ti es alinear los índices que calculas dentro del shader.

Por filas o por columnas: la localidad que cuenta no es la que crees

El experimento canónico es sumar cada fila de una matriz de N por N en f32 guardada por filas. Hay dos formas obvias y una de ellas es un desastre.

const N: u32 = 4096u;

@group(0) @binding(0) var<storage, read>       m: array<f32>;
@group(0) @binding(1) var<storage, read_write> sumas: array<f32>;

// VERSION A: un hilo por fila. Cada hilo recorre su fila entera.
@compute @workgroup_size(64)
fn porFila(@builtin(global_invocation_id) gid: vec3u) {
  let fila = gid.x;
  if (fila >= N) { return; }
  var s = 0.0;
  for (var k: u32 = 0u; k < N; k = k + 1u) {
    s = s + m[fila * N + k];       // hilos vecinos leen a 16 KB de distancia
  }
  sumas[fila] = s;
}

var<workgroup> parcial: array<f32, 64>;

// VERSION B: un workgroup por fila. Los 64 hilos recorren la fila en paralelo.
@compute @workgroup_size(64)
fn porFilaCooperativo(@builtin(workgroup_id) wid: vec3u,
                      @builtin(local_invocation_index) li: u32) {
  let fila = wid.x;
  var s = 0.0;
  var k = li;
  loop {
    if (k >= N) { break; }
    s = s + m[fila * N + k];       // hilos vecinos leen posiciones vecinas
    k = k + 64u;
  }
  parcial[li] = s;
  workgroupBarrier();

  // Reduccion en arbol dentro del grupo.
  var paso = 32u;
  loop {
    if (paso == 0u) { break; }
    if (li < paso) { parcial[li] = parcial[li] + parcial[li + paso]; }
    workgroupBarrier();
    paso = paso / 2u;
  }
  if (li == 0u) { sumas[fila] = parcial[0]; }
}

Las dos versiones leen exactamente los mismos N² valores y producen el mismo resultado. La versión A es la que escribiría cualquiera que venga de la CPU, porque recorre la memoria por filas y en la CPU eso es lo correcto. En la GPU es la mala: en cada iteración, los hilos 0 y 1 leen direcciones separadas por N por 4 bytes, es decir 16 KB, así que un grupo de 32 carriles genera 32 transacciones para consumir 128 bytes útiles. La versión B lee, en cada iteración, 64 valores consecutivos: dos transacciones perfectas.

La localidad que le importa a una GPU es la que hay entre invocaciones vecinas en el mismo instante, no la que hay dentro de una invocación a lo largo del tiempo. Ese cambio de eje mental es la mitad de la optimización de memoria en compute.

Si montas el experimento con N igual a 4096, que son 64 MiB de matriz, la diferencia medida estará entre ocho y veinte veces según la GPU. Conviene entender por qué no llega a las treinta y dos que predice la teoría, porque el motivo es instructivo: en la versión A, las líneas que se traen en la iteración k contienen también los valores de las iteraciones k+1 a k+31 de esa misma fila. El conjunto de trabajo instantáneo es de unos 512 KB, y si la L2 de tu GPU es mayor que eso, la caché acaba capturando buena parte de la reutilización y convierte el desastre en algo meramente malo. Sube N o ejecútalo en una GPU con L2 pequeña y el factor se acerca al teórico. Es un buen recordatorio de que las cachés perdonan los patrones malos justo hasta el tamaño en que dejan de perdonarlos, que suele ser el tamaño de producción.

Array de estructuras frente a estructura de arrays

El patrón de acceso no se arregla del todo en el kernel: se decide cuando eliges cómo se guardan los datos. Un sistema de partículas es el caso claro.

// AoS: un array de estructuras. 64 bytes por particula.
struct Particula {
  pos:   vec4f,   // 16 B
  vel:   vec4f,   // 16 B
  color: vec4f,   // 16 B
  attrs: vec4f,   // vida, tamano, semilla, masa
};
@group(0) @binding(0) var<storage, read_write> particulas: array<Particula>;

// SoA: una estructura de arrays. Los mismos 64 bytes, repartidos.
@group(0) @binding(0) var<storage, read_write> pos:   array<vec4f>;
@group(0) @binding(1) var<storage, read_write> vel:   array<vec4f>;
@group(0) @binding(2) var<storage, read_write> color: array<vec4f>;
@group(0) @binding(3) var<storage, read_write> attrs: array<vec4f>;

Ahora pasa un kernel que solo necesita las posiciones: el que asigna cada partícula a una celda de una rejilla espacial, o el que calcula la caja envolvente. Con AoS, los hilos vecinos leen pos a 64 bytes de distancia entre sí, así que el bus trae las líneas completas y de cada 64 bytes movidos usas 16: el 25 por ciento. Con SoA, pos es un array denso y los hilos vecinos leen posiciones vecinas: el 100 por cien.

Con un millón de partículas la cuenta es directa. AoS mueve 64 MB para usar 16; SoA mueve 16 MB. A 400 GB/s eso son 160 microsegundos frente a 40. El kernel es idéntico salvo por dónde lee.

La conclusión no es “usa siempre SoA”, y aquí es donde la mayoría de los artículos se quedan cortos. El criterio es qué fracción de la estructura toca cada kernel:

• Si un kernel toca todos los campos de una partícula y ninguno más, AoS es igual de eficiente —trae 64 bytes y usa 64— y encima gasta un solo binding.
• Si tienes kernels que tocan subconjuntos distintos, SoA gana en todos ellos.
• SoA cuesta bindings, y ahí hay un límite duro: maxStorageBuffersPerShaderStage son 8 por defecto y maxBindGroups son 4. Doce arrays separados no caben.

El compromiso que usa la industria es partir por temperatura, no por campo: un buffer con lo que se toca en cada frame —posición y velocidad— y otro con lo que se toca raramente —color, semilla, material—. Dos bindings, la eficiencia casi de SoA, y sitio de sobra dentro de los límites. Y vigila el tamaño: maxStorageBufferBindingSize son 134217728 bytes, 128 MiB, así que con partículas de 64 bytes el techo de una sola vista está en dos millones.

Texturas y alineación: el tráfico que no ves

Las texturas juegan con otras reglas, y por una buena razón. Una imagen no se recorre en una dimensión sino en dos: los píxeles vecinos en pantalla necesitan téxeles vecinos en las dos direcciones. Si una textura se guardara linealmente por filas, una línea de caché de 128 bytes cubriría 32 téxeles RGBA8 en una fila y ninguno de la fila de abajo, y el filtrado bilineal —que necesita un bloque de 2 por 2— tocaría dos líneas siempre.

Por eso el hardware guarda las texturas con un orden en zigzag, entrelazando los bits de las coordenadas para que un bloque cuadrado de téxeles caiga en direcciones contiguas. Con ese orden, una línea cubre algo parecido a un bloque de 8 por 4 téxeles y la localidad bidimensional se convierte en localidad en memoria. Las consecuencias prácticas son tres:

1. Muestrear con coordenadas coherentes es barato. Si los píxeles vecinos muestrean téxeles vecinos, todo el grupo SIMD comparte unas pocas líneas. Con coordenadas aleatorias, cada muestra trae su línea entera para usar 4 bytes: la misma eficiencia del 3 por ciento de la tabla de arriba.
2. Los mipmaps son una optimización de ancho de banda antes que de aliasing. Cuando una textura se ve minificada, los téxeles que corresponden a píxeles vecinos están lejos entre sí en el nivel 0, y el conjunto de trabajo se dispara. Muestrear el nivel adecuado restaura la coherencia. Cuesta un 33 por ciento más de memoria y puede dividir el tráfico de texturas por un orden de magnitud.
3. El orden es opaco y no lo controlas. Lo único que asoma de él en la API es la obligación de que bytesPerRow sea múltiplo de 256 al copiar entre un buffer y una textura: esa conversión de lineal a zigzag es trabajo real que hace el driver.

Del lado de WGSL, las reglas de alineación mueven bytes que no usas. La que más duele: vec3<f32> se alinea a 16 bytes aunque solo ocupe 12. Un struct con dos vec3<f32> seguidos ocupa 32 bytes, de los que 8 son relleno: un 25 por ciento de tu ancho de banda transportando nada. La solución es meter escalares en los huecos, que es gratis:

struct Malo  { pos: vec3f, nor: vec3f };                        // 32 B, 8 de relleno
struct Bueno { pos: vec3f, vida: f32, nor: vec3f, tam: f32 };   // 32 B, 0 de relleno

Y una regla del espacio uniform que sorprende a todo el mundo la primera vez: en ese espacio de direcciones, el paso entre elementos de un array se redondea al múltiplo de 16 más cercano por arriba. Un array<f32, 64> declarado en un var<uniform> no ocupa 256 bytes: ocupa 1024, porque cada f32 se lleva 16 bytes de los que 12 son relleno. Cuadruplicas la memoria y el tráfico sin escribir una línea distinta. Si necesitas un array de escalares en un uniform, declara array<vec4f, 16> y desempaqueta a mano, o usa un storage buffer, donde la regla no se aplica.