Paralelismo, SIMD y grandes volúmenes

Un bucle que recorre cincuenta millones de double hace siempre el mismo trabajo por elemento. Si tu equipo tiene diez núcleos, nueve están esperando mientras uno suma. Paralelizar consiste en repartir ese rango en trozos; usar SIMD consiste en que cada instrucción sume varios valores a la vez. Son dos formas distintas de hacer más con el mismo tiempo, y se combinan.

¿Qué es paralelizar y qué es SIMD?

El paralelismo de hilos divide los datos entre núcleos: cada hilo ejecuta el mismo código sobre una parte del rango. El SIMD (single instruction, multiple data) actúa dentro de un núcleo: una sola instrucción opera sobre un registro de 128 o 256 bits que contiene dos, cuatro u ocho valores. El primero multiplica la capacidad de trabajo; el segundo reduce el número de instrucciones.

¿Para qué sirve con grandes volúmenes?

Con cien filas, paralelizar no compensa: crear hilos cuesta microsegundos y el reparto se come la ganancia. Con cien millones de filas, la situación se invierte. Este es el terreno donde C++ se separa de Python: un for de Python con diez millones de iteraciones tarda cerca de un segundo, mientras que el mismo cálculo vectorizado en C++ tarda milisegundos.

• Normalizar o escalar columnas de una matriz antes de entrenar: cada columna es independiente.
• Calcular la pérdida sobre un lote grande: la suma de errores se puede repartir.
• Transformar características (logaritmos, cuadrados, recortes) elemento a elemento.
• Multiplicar matrices por bloques: cada bloque del resultado se calcula por separado.

La condición para que todo esto funcione es la misma en los cuatro casos: que una iteración no dependa del resultado de otra. Si el cálculo de la fila cinco necesita la fila cuatro, el bucle no se puede repartir tal cual y hay que replantearlo.

Cómo funciona por dentro: núcleos, caché y ancho de banda

Un procesador de diez núcleos no es diez veces más rápido. Cada núcleo tiene su propia caché L1 y L2; si dos núcleos escriben variables distintas que caen en la misma línea de 64 bytes, la coherencia de caché invalida la línea en el otro núcleo en cada escritura. Ese fenómeno se llama false sharing y convierte un ++ en un viaje a memoria.

El segundo límite es el ancho de banda. Sumar cincuenta millones de double son 400 MB de datos. Medido en un Apple M1 Max, la versión secuencial tarda unos 47 ms, es decir, unos 8.5 GB/s: muy por debajo de lo que la memoria puede entregar, porque el cuello de botella es la cadena de sumas, no el transporte. Con cuatro acumuladores baja a unos 12 ms, unos 34 GB/s. Cuando el cómputo deja de estorbar, el límite pasa a ser la memoria y añadir hilos ya no ayuda.

Por eso el orden de las optimizaciones importa: primero el algoritmo, después la localidad de los datos, luego SIMD y al final los hilos. Paralelizar un bucle que desperdicia caché reparte el desperdicio entre más núcleos.

Sintaxis: OpenMP, std::execution y pragmas SIMD

OpenMP es un conjunto de directivas #pragma que el compilador traduce a hilos. La forma canónica reparte un for y combina los resultados parciales con reduction.

1234567891011121314151617181920

#include <chrono>
#include <cstdio>
#include <vector>

int main() {
    std::vector<double> v(50'000'000, 1.0);

    const auto t0 = std::chrono::steady_clock::now();
    double suma = 0.0;

#pragma omp parallel for reduction(+ : suma) schedule(static)
    for (std::size_t i = 0; i < v.size(); ++i) {
        suma += v[i];
    }

    const auto t1 = std::chrono::steady_clock::now();
    const double ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
    std::printf("suma = %.0f  (%.1f ms)\n", suma, ms);
    return 0;
}

123456789

# Linux con GCC: el flag activa el runtime
g++ -std=c++20 -O2 -Wall -Wextra -fopenmp suma_omp.cpp -o suma_omp

# macOS con Apple clang: hace falta libomp instalado con Homebrew
brew install libomp
g++ -std=c++20 -O2 -Xpreprocessor -fopenmp -lomp suma_omp.cpp -o suma_omp

# Cuántos hilos se usan si no lo cambias
OMP_NUM_THREADS=4 ./suma_omp

std::execution ofrece la misma idea dentro de la biblioteca estándar: pasas std::execution::par a un algoritmo y el runtime decide cómo repartirlo. Es más legible que OpenMP, pero depende de que tu biblioteca estándar lo implemente.

12345678910111213141516

#include <algorithm>
#include <cstdio>
#include <execution>
#include <numeric>
#include <vector>

int main() {
    std::vector<double> v(10'000'000, 2.0);

    std::transform(std::execution::par, v.begin(), v.end(), v.begin(),
                   [](double x) { return x * x; });

    const double total = std::reduce(std::execution::par, v.begin(), v.end(), 0.0);
    std::printf("%.0f %.0f\n", v[0], total);  // 4 40000000
    return 0;
}

• En Linux con GCC necesitas la biblioteca TBB y enlazarla con -ltbb; sin ella el enlazado falla.
• La biblioteca estándar de Apple no incluye <execution>: ese programa no compila en macOS con clang. Compruébalo antes de diseñar alrededor de ella.
• std::reduce admite reordenar la suma; std::accumulate no. Para paralelizar una suma usa reduce.

Para SIMD no necesitas bibliotecas: el compilador vectoriza si detecta que las iteraciones son independientes. Puedes insinuarlo con #pragma omp simd y comprobar el resultado en el ensamblador, donde aparecerán instrucciones con sufijo pd o ps en x86 y v en ARM.

Ejemplo básico: la misma suma, secuencial y en paralelo

Esta es la versión secuencial del programa anterior, con cronómetro, para que puedas comparar los dos binarios en tu equipo. El resultado debe ser idéntico: 50000000.

123456789101112131415161718

#include <chrono>
#include <cstdio>
#include <vector>

int main() {
    std::vector<double> v(50'000'000, 1.0);

    const auto t0 = std::chrono::steady_clock::now();
    double suma = 0.0;
    for (double x : v) {
        suma += x;
    }
    const auto t1 = std::chrono::steady_clock::now();

    const double ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
    std::printf("suma = %.0f  (%.1f ms)\n", suma, ms);
    return 0;
}

Con -O2 en un Apple M1 Max, la versión secuencial tarda unos 47 ms. La paralela con diez hilos baja a la fracción que permita el ancho de banda. La diferencia entre las dos cifras es tu speedup real, y casi nunca coincide con el número de núcleos.

Ejemplo intermedio: cuatro carriles sin cambiar de lenguaje

La suma en coma flotante no es asociativa, así que el compilador no puede reordenarla por su cuenta: (a + b) + c no siempre da lo mismo que a + (b + c). Si escribes cuatro acumuladores independientes, le das permiso explícito para trabajar en paralelo dentro de un solo hilo.

12345678910111213141516171819202122232425262728

#include <chrono>
#include <cstdio>
#include <vector>

double suma_cuatro(const std::vector<double>& v) {
    double a = 0.0, b = 0.0, c = 0.0, d = 0.0;
    const std::size_t n = v.size();
    std::size_t i = 0;
    for (; i + 4 <= n; i += 4) {
        a += v[i];
        b += v[i + 1];
        c += v[i + 2];
        d += v[i + 3];
    }
    for (; i < n; ++i) a += v[i];  // resto cuando n no es múltiplo de 4
    return a + b + c + d;
}

int main() {
    std::vector<double> v(50'000'000, 1.0);
    const auto t0 = std::chrono::steady_clock::now();
    const double s = suma_cuatro(v);
    const auto t1 = std::chrono::steady_clock::now();

    const double ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
    std::printf("%.0f en %.1f ms\n", s, ms);
    return 0;
}

Medido con -O2 en el mismo equipo: unos 12 ms frente a los 47 de la versión con un solo acumulador. Con #pragma omp simd antes del bucle, algunos compiladores van un paso más allá y empaquetan dos double por instrucción en ARM o cuatro en x86 con AVX2.

El detalle del resto (i < n) no es adorno: si el tamaño no es múltiplo de cuatro, los últimos elementos se quedarían sin sumar. Ese error produce resultados que parecen correctos porque la diferencia es pequeña, y es exactamente el tipo de fallo que ninguna prueba con datos redondos detecta.

Visualización: un bloque de datos repartido entre dos hilos

La operación es normalizar cada columna entre 0 y 1. Cada hilo recibe un bloque de filas, calcula los mínimos y máximos parciales de su bloque y luego se combinan. Recorre los pasos y fíjate en que la normalización final se aplica sobre los extremos globales, no sobre los locales.

El patrón de la reducción es el mismo que hace OpenMP por dentro con reduction: cada hilo acumula en su copia privada y al final se combinan. La diferencia es que aquí hay dos operaciones (mínimo y máximo) en lugar de una suma, así que se hace a mano.

Caso práctico: false sharing

Dos hilos que incrementan contadores distintos no deberían estorbarse. Si esos contadores comparten línea de caché, sí lo hacen: cada escritura invalida la línea del otro. El programa siguiente mide el mismo trabajo con los contadores juntos y separados 64 bytes.

12345678910111213141516171819202122232425262728293031323334

#include <chrono>
#include <cstdio>
#include <thread>

struct Juntos {
    volatile long long a = 0;
    volatile long long b = 0;  // misma linea de caché que a
};

struct Separados {
    alignas(64) volatile long long a = 0;
    alignas(64) volatile long long b = 0;  // cada uno en su linea
};

constexpr long long N = 50'000'000;

template <typename T>
double medir(T& c) {
    const auto t0 = std::chrono::steady_clock::now();
    std::thread h1([&] { for (long long i = 0; i < N; ++i) c.a = c.a + 1; });
    std::thread h2([&] { for (long long i = 0; i < N; ++i) c.b = c.b + 1; });
    h1.join();
    h2.join();
    return std::chrono::duration<double, std::milli>(std::chrono::steady_clock::now() - t0).count();
}

int main() {
    Juntos juntos;
    Separados separados;
    std::printf("misma linea   %.1f ms\n", medir(juntos));
    std::printf("lineas aparte %.1f ms\n", medir(separados));
    std::printf("%lld %lld\n", juntos.a, separados.a);
    return 0;
}

Resultado en un Apple M1 Max, con el mejor de cinco intentos: 78 ms con los contadores juntos y 17 ms separados, casi cinco veces menos, con el mismo número de incrementos. Los dos contadores acaban en 50000000, así que la corrección era la misma desde el principio: lo que cambiaba era el precio. alignas(64) obliga a cada contador a empezar en su propia línea de caché.

El volatile está ahí para impedir que el compilador reduzca el bucle a c.a += N: sin él, medirías una sola suma en lugar de cincuenta millones. Es una herramienta de medición, no una forma de sincronizar hilos; para eso están std::atomic y los mutex.

Errores comunes al paralelizar

• Escribir en una variable compartida sin reducción: suma += v[i] desde varios hilos pierde sumas y no da error. Usa reduction o un acumulador por hilo.
• Crear un hilo por elemento: el coste de crear y juntar hilos supera el trabajo. Reparte bloques grandes.
• Olvidar el flag -fopenmp: el #pragma se ignora en silencio y crees haber paralelizado.
• Ignorar el resto del rango: si n no es múltiplo del ancho del bloque, los últimos elementos se quedan fuera.
• Medir con un solo hilo: sin OMP_NUM_THREADS tu prueba puede estar usando un núcleo y parecer que no escala.

Buenas prácticas

• Paraleliza solo bucles cuyo cuerpo no dependa de iteraciones anteriores; compruébalo antes de escribir el pragma.
• Verifica que el resultado paralelo coincide con el secuencial, con tolerancia para la coma flotante.
• Reparte por bloques contiguos: es más amable con la caché y con la memoria.
• Evita que dos hilos escriban en la misma línea de caché; separa los contadores con alignas(64).
• Mide el speedup real con varios números de hilos y quédate con el punto donde deja de mejorar.
• Deja la versión secuencial en el código, accesible con un flag, para poder comparar y depurar.

Ejercicios

double suma = 0.0;
#pragma omp  for  (+ : suma) schedule(static)
for (std::size_t i = 0; i < v.size(); ++i) {
    suma += v[i];
}

Soluciones explicadas

En el primer ejercicio, la comparación con tolerancia no es un detalle académico: la suma paralela agrupa los términos en distinto orden y el resultado puede diferir en los últimos bits. Comparar con == te habría hecho creer que la reducción está mal cuando solo cambió el redondeo.

En el segundo, la lección es que el rendimiento puede depender de la posición en memoria, no de la lógica. Ese es el motivo por el que las bibliotecas serias alinean sus estructuras internas y rellenan los objetos compartidos entre hilos: están comprando líneas de caché enteras para no compartirlas.

Resumen

• Paralelizar reparte el rango entre núcleos; SIMD procesa varios valores por instrucción. Se combinan.
• #pragma omp parallel for reduction(+ : suma) es la forma mínima de sumar en paralelo.
• Sin -fopenmp el pragma se ignora y el programa sigue funcionando en silencio, pero secuencial.
• El false sharing hace que dos hilos se invaliden la caché mutuamente: separa los datos con alignas(64).
• La suma en coma flotante no es asociativa; cuatro acumuladores dan paralelismo sin cambiar el resultado esperado.
• El speedup real lo limitan la memoria y las dependencias, no el número de núcleos.
• Comprueba siempre que el resultado paralelo coincide con el secuencial.

Evaluación

12345

import numpy as np

x = np.full(10_000_000, 2.0)
y = x * x          # bucle en C, vectorizado con SIMD
print(y[0], y.sum())

123456

std::vector<double> x(10'000'000, 2.0);
std::vector<double> y(x.size());

for (std::size_t i = 0; i < x.size(); ++i) {
    y[i] = x[i] * x[i];  // el compilador puede vectorizarlo
}