SIMD y vectorización

Cada núcleo de tu máquina lleva registros de doscientos cincuenta y seis o quinientos doce bits capaces de sumar ocho o dieciséis floats en una sola instrucción. Si tus bucles no los usan, estás dejando sin tocar entre el ochenta y siete y el noventa y cuatro por ciento del ancho aritmético que compraste. La buena noticia es que el compilador quiere usarlos por ti. La mala es que solo lo hace cuando puede demostrar que es correcto, y la mayoría de los bucles de C no le dan material para esa demostración.

Por qué tu bucle no se vectoriza

Vectorizar significa reescribir un bucle escalar de n iteraciones como uno de n/w iteraciones que procesan w elementos de golpe. Esa transformación solo es legal si las iteraciones pueden ejecutarse en cualquier orden y a la vez. El auto-vectorizador pasa la mayor parte de su tiempo intentando probar eso y fallando por alguna de estas razones:

• Dependencias entre iteraciones. Si la iteración i lee lo que escribió la i - 1, agrupar iteraciones cambia el resultado. Una dependencia hacia atrás cierra la puerta; una hacia delante con distancia mayor que el ancho del vector no.
• Aliasing posible. Dados dos punteros a float, el compilador debe asumir que pueden solaparse. Si lo hacen, escribir ocho resultados de golpe pisaría entradas aún no leídas. Sin garantía de no solapamiento, o se rinde o emite comprobaciones en ejecución que duplican el código.
• Reducciones en coma flotante. Sumar ocho parciales y combinarlos al final cambia el orden de las sumas, y la suma en coma flotante no es asociativa. El compilador respeta la semántica de la norma y no lo hace salvo que se lo autorices.
• Control de flujo dentro del cuerpo. Una rama por iteración se puede vectorizar con máscaras, pero un break cuya condición dependa de los datos, una llamada no expandible en línea o un goto hacia fuera lo hacen imposible.
• Recorrido no unitario. Un paso constante distinto de uno obliga a instrucciones de recolección, mucho más lentas que una carga contigua. Un patrón indexado por otro array las obliga siempre.
• Cuenta desconocida o con signo confuso. El compilador necesita razonar sobre el número de iteraciones sin que la aritmética del contador pueda envolver.

flowchart TD
A[Bucle candidato] --> B[Hay dependencia entre iteraciones]
B -->|Si| Z[No vectoriza: reestructura el algoritmo]
B -->|No| C[Pueden solaparse los punteros]
C -->|Si| D[Anade restrict o acepta comprobacion en ejecucion]
C -->|No| E[Es una reduccion en coma flotante]
E -->|Si| F[Autoriza reasociacion acotada al bucle]
E -->|No| G[El acceso es contiguo]
G -->|No| H[Reorganiza los datos hacia SoA]
G -->|Si| I[Vectoriza]
style I fill:#a6e3a1,color:#11111b
style Z fill:#f38ba8,color:#11111b

El diálogo con el auto-vectorizador

Lo primero es dejar de adivinar. Los compiladores dicen exactamente qué bucle rechazaron y por qué:

# GCC: informe de bucles vectorizados y de los que no, con el motivo
gcc -std=c23 -O3 -march=native -fopt-info-vec-optimized -fopt-info-vec-missed prog.c

# Clang: el mismo diagnostico por el canal de remarks
clang -std=c23 -O3 -march=native -Rpass=loop-vectorize \
      -Rpass-missed=loop-vectorize -Rpass-analysis=loop-vectorize prog.c

Un detalle de compilación importa tanto como el código: sin -march el compilador genera para la línea base de x86-64, que solo garantiza SSE2, es decir vectores de ciento veintiocho bits. -march=native desbloquea AVX2 o AVX-512 en la máquina donde compilas; -march=x86-64-v3 es el compromiso razonable para binarios distribuibles, porque exige AVX2 y FMA sin atarse a un modelo concreto.

Con el diagnóstico en la mano, la intervención habitual es prometerle al compilador lo que no puede deducir:

#include <stddef.h>

/* Version que el compilador no puede vectorizar sin comprobaciones:
   a, b y c podrian solaparse. */
void saxpy_lento(float *c, const float *a, const float *b, float k, size_t n) {
    for (size_t i = 0; i < n; i++)
        c[i] = k * a[i] + b[i];
}

/* Con restrict prometes que los tres rangos son disjuntos.
   El compilador emite un bucle vectorial limpio, sin prologo de guardas. */
void saxpy(float *restrict c, const float *restrict a,
           const float *restrict b, float k, size_t n) {
    for (size_t i = 0; i < n; i++)
        c[i] = k * a[i] + b[i];
}

restrict es un contrato que tú garantizas y el compilador no verifica: si mientes, el comportamiento es indefinido y el fallo será silencioso y dependiente del nivel de optimización. Es el precio de convertir una prueba imposible en una premisa.

Las reducciones necesitan un permiso distinto, porque el problema no es el aliasing sino la aritmética. Autorizar la reasociación globalmente con -ffast-math es una mala idea —arrastra suposiciones sobre NaN, infinitos y ceros con signo por todo el programa—, así que la forma disciplinada es acotarla al bucle:

#include <stddef.h>

float suma(const float *restrict v, size_t n) {
    float s = 0.0f;
    #pragma omp simd reduction(+:s)     /* autoriza reasociar SOLO aqui */
    for (size_t i = 0; i < n; i++)
        s += v[i];
    return s;
}

Compilado con -fopenmp-simd esto no crea hilos: solo activa la directiva vectorial. El resultado en coma flotante diferirá del escalar en los últimos bits, y esa diferencia es precisamente lo que estabas autorizando. Si necesitas reproducibilidad exacta, no lo hagas; si necesitas precisión, la suma por parciales suele ser más exacta que la secuencial, porque acumula en cadenas más cortas.

Intrínsecos, y cómo no perder portabilidad

Cuando el algoritmo no encaja en un bucle regular —mezclas, permutaciones, búsquedas de bytes, criptografía— los intrínsecos son la herramienta correcta. Son funciones que se traducen una a una a instrucciones de la máquina, con el compilador ocupándose todavía del reparto de registros y la planificación:

#include <immintrin.h>
#include <stddef.h>

/* Producto escalar con AVX2 y FMA: cuatro acumuladores para ocultar
   la latencia de la instruccion de multiplicacion y suma fundida. */
float dot_avx2(const float *a, const float *b, size_t n) {
    __m256 s0 = _mm256_setzero_ps(), s1 = _mm256_setzero_ps();
    size_t i = 0;
    for (; i + 16 <= n; i += 16) {
        s0 = _mm256_fmadd_ps(_mm256_loadu_ps(a + i),      _mm256_loadu_ps(b + i),      s0);
        s1 = _mm256_fmadd_ps(_mm256_loadu_ps(a + i + 8),  _mm256_loadu_ps(b + i + 8),  s1);
    }
    __m256 s = _mm256_add_ps(s0, s1);
    float t[8];
    _mm256_storeu_ps(t, s);
    float r = t[0] + t[1] + t[2] + t[3] + t[4] + t[5] + t[6] + t[7];
    for (; i < n; i++) r += a[i] * b[i];   /* cola escalar, siempre necesaria */
    return r;
}

Dos detalles que separan el código de manual del código de producción. El primero son los acumuladores múltiples: _mm256_fmadd_ps tiene una latencia de cuatro ciclos y un rendimiento de uno o dos por ciclo, así que con un solo acumulador la cadena de dependencias limita el bucle a una operación cada cuatro ciclos y desperdicias tres cuartas partes de la unidad. El segundo es la cola escalar: todo bucle vectorial la necesita y todo bucle vectorial escrito con prisa la olvida.

El problema de los intrínsecos es que un binario compilado con -mavx2 muere con instrucción ilegal en una máquina sin AVX2. La solución profesional es el despacho en tiempo de ejecución: compilar varias versiones de la misma función, cada una con su conjunto de instrucciones, y elegir al arrancar.

__attribute__((target("avx2,fma")))
float dot_v3(const float *a, const float *b, size_t n) { /* ... */ }

__attribute__((target("default")))
float dot_base(const float *a, const float *b, size_t n) { /* ... */ }

static float (*dot)(const float *, const float *, size_t);

static void seleccionar(void) {
    __builtin_cpu_init();
    dot = __builtin_cpu_supports("avx2") ? dot_v3 : dot_base;
}

El atributo target permite tener AVX2 en una función concreta sin exigirlo en toda la unidad de traducción, y __builtin_cpu_supports consulta la identificación del procesador. Existe además una vía intermedia y muy infravalorada: las extensiones vectoriales genéricas de GCC y Clang, que definen tipos vectoriales con los operadores aritméticos habituales y compilan a las instrucciones nativas de cada arquitectura, sean AVX, NEON o RVV. Dan el ochenta por ciento del rendimiento de los intrínsecos con código legible que compila en ARM y en x86 sin ramificar el fuente.

Un aviso final para cerrar el círculo con la lección anterior. Vectorizar solo acelera lo que está limitado por cómputo. Si tu bucle recorre un array de cien mebibytes sumando elementos, el límite es el ancho de banda de DRAM y las instrucciones vectoriales no lo mueven ni un uno por ciento: se limitan a esperar más rápido. Antes de vectorizar, comprueba con los contadores en qué régimen estás.