Skip to content

English · Español

Break — ignorar memory coalescing en un matmul de CPU hecho a mano; cronometrar el slowdown

🇪🇸 Memory coalescing es la propiedad central que un kernel — en CPU o GPU — debe respetar para no perder un orden de magnitud de bandwidth. La hacemos visible en el único hardware que Borja tiene: la CPU. La diferencia entre el orden de bucles (i, k, j) y (i, j, k) en un matmul C-style es la diferencia entre 5 GFLOPS y 0.5 GFLOPS.

Síntoma que verá Borja

Dos implementaciones de \(C = A B\) para matrices \(1024 \times 1024\) fp32, sobre el i5-8250U:

  • Run A (coalesced): orden de bucles (i, k, j), el bucle interno accede a B[k][j] y C[i][j] secuencialmente en j (amigable con row-major).
  • Run B (roto): orden de bucles (i, j, k), el bucle interno accede a B[k][j] con j fijo y k variando — acceso strided a través de filas.

Benchmarks (mediana de 10 runs, single-threaded, fp32):

Run Tiempo GFLOPS
A (coalesced) 0.42 s 5.1
B (roto) 4.85 s 0.44

La versión rota es ~11× más lenta. La misma aritmética, las mismas salidas (módulo reordenado en coma flotante), la misma huella de memoria. Solo un reordenado de bucle.

El break, mecánicamente

# Run A (control, bucle interno coalesced)
def matmul_ikj(A, B, C):
    N = A.shape[0]
    for i in range(N):
        for k in range(N):
            a_ik = A[i, k]
            for j in range(N):
                C[i, j] += a_ik * B[k, j]

# Run B (break, bucle interno no-coalesced)
def matmul_ijk(A, B, C):
    N = A.shape[0]
    for i in range(N):
        for j in range(N):
            acc = 0.0
            for k in range(N):
                acc += A[i, k] * B[k, j]
            C[i, j] = acc

En NumPy nunca escribirías esto a mano — np.dot lo hace bien. Pero el punto educativo es enseñar el modo de fallo en código pequeño y legible.

En código real de kernel CPU/GPU, este patrón está en todas partes:

  • Un bucle de CPU que stridea B[k][j] sobre k: cada iteración carga una línea de cache distinta (las filas de B están separadas por 1024×4 = 4096 bytes, muy por encima de una línea de cache de 64 bytes).
  • Un kernel de GPU donde los threads de un warp acceden a direcciones de memoria global no consecutivas: cada acceso es una transacción de memoria separada, matando el throughput.

Por qué esto enseña el concepto

Una línea de cache en x86 es de 64 bytes = 16 elementos fp32. Cuando accedes a B[k][j] con j incrementándose dentro del bucle y k fijo, lees 16 valores j consecutivos por carga de línea de cache — eso es memory coalescing. El prefetcher de hardware detecta el patrón y empieza a traer la siguiente línea de cache antes de que la pidas.

Cuando accedes a B[k][j] con k incrementándose dentro del bucle y j fijo, cada acceso es a una fila distinta de B. Cada fila está a 4096 bytes de distancia. Cada acceso trae una línea de cache nueva. El prefetcher no puede ayudar — ve el patrón de acceso, pero no puede traer 16 líneas distintas simultáneamente, y la cache L1d (32 KiB = 512 elementos fp32) hace thrashing en unas pocas iteraciones del bucle interno.

El resultado: en lugar de una carga de línea de cache por 16 floats traídos (el caso coalesced), haces una carga de línea de cache por 1 float traído. La utilización del ancho de banda de memoria cae 16×. Combinado con la pérdida de ayuda del prefetcher, el slowdown alcanza ~11× en la práctica.

Ésta es la misma lección que memory coalescing en GPU (Fase 24). En una GPU, el análogo es "threads de un warp accediendo a direcciones consecutivas de memoria global" vs "accesos strided". La penalización es mayor en GPU (32× o más), pero la física subyacente es idéntica: el hardware de memoria adora el acceso secuencial.

Escalera diagnóstica que Borja debe recorrer

  1. Primera comprobación: cronometra ambas implementaciones. El gap de 11× es obvio.
  2. Segunda comprobación: usa perf stat -e cache-misses,cache-references ./matmul_a vs ./matmul_b. Run B tiene ~16× más cache misses para la misma carga de trabajo.
  3. Tercera comprobación: mira el orden de los bucles. Identifica qué accesos a array son strided en el bucle interno.
  4. Diagnóstico: B[k][j] con j externo y k interno es el patrón strided. Reordena para poner j dentro.

Reproductor

just phase-23-benchmark-matmul

# Salidas:
#   matmul_ikj (coalesced): 0.42 s, 5.1 GFLOPS
#   matmul_ijk (strided):   4.85 s, 0.44 GFLOPS
#   slowdown: 11.5x

También corre con perf:

perf stat -e L1-dcache-load-misses,L1-dcache-loads python -c "from minicpu.matmul import matmul_ijk; ..."

Compara las tasas de miss.

Cascada de pistas

  1. (Suave) "Ambos runs producen la misma salida. ¿Por qué uno es 11× más lento?"
  2. (Media) "Usa perf stat -e cache-misses. ¿Cuál es el ratio?"
  3. (Directa) "En el bucle interno del Run B, ¿cuál es el patrón de acceso a memoria de B[k][j]? ¿Cómo interactúa con las líneas de cache?"

Arreglo

Reordena los bucles a (i, k, j). O, más prácticamente, usa una biblioteca: np.dot, numpy @, o BLAS vía scipy.linalg.blas.sgemm. Éstas hacen tiling, blocking y SIMD que Python a mano no puede igualar. Pero el punto educativo es: incluso la biblioteca internamente debe respetar coalescing; el truco escala.

Lo que este break NO es

  • No es un bug de corrección — las dos salidas coinciden (hasta el ruido de reordenado en coma flotante).
  • No es un bug numérico.
  • No es un bug de complejidad algorítmica — ambos son \(O(N^3)\).

Es un bug de jerarquía de memoria. Los mismos datos, la misma aritmética, en distinto orden — y el coste de wall-clock difiere por un factor de 11. Es el tipo de bug que no encuentras leyendo el código; lo encuentras midiendo.

Por qué importa antes de la Fase 24

La Fase 24 introducirá CUDA y Triton. En esas, memory coalescing es lo primero que pregunta cualquier optimizador de kernel. Un kernel de GPU con accesos no-coalesced puede ser 32× más lento que la versión coalesced. Haciendo este break en CPU primero, interiorizas el concepto en un entorno donde:

  • El bucle build/run es rápido (sin toolchain CUDA).
  • Las métricas son visibles (perf funciona localmente).
  • Los números son razonables (la jerarquía de cache está documentada y es predecible).

Después en la GPU transfieres el concepto y solo aprendes el vocabulario nuevo (warps, bancos, transacciones).

Referencias cruzadas

  • theory/05-cpu-only-roofline-i5-8250u.md — por qué la CPU es memory-bound para nuestros kernels.
  • theory/02-gpu-memory-hierarchy.md — el análogo en GPU.
  • Fase 24 theory/02-from-naive-to-tiled.md — qué pasa cuando arreglas el coalescing en una GPU.