El scan de Hillis-Steele

El scan, o suma de prefijos, convierte un array en otro donde cada posición contiene la combinación de todo lo que había hasta ella. Parece una curiosidad y es el segundo ladrillo del cómputo paralelo: sin scan no hay compactación de arrays, ni ordenación por radix, ni construcción de rejillas espaciales, ni asignación de memoria en paralelo. Su formulación más simple tiene un nombre y una fecha: Hillis y Steele la publicaron en 1986 para la Connection Machine, y sigue siendo la que hay que escribir la primera vez.

Inclusivo, exclusivo y para qué sirve

Dado un array x y un operador asociativo, el scan inclusivo produce y[i] = x[0] op x[1] op ... op x[i], y el exclusivo produce y[i] = x[0] op ... op x[i-1], con el elemento neutro en la posición cero.

x           = [3, 1, 7, 0, 4, 1, 6, 3]
inclusivo   = [3, 4, 11, 11, 15, 16, 22, 25]
exclusivo   = [0, 3, 4, 11, 11, 15, 16, 22]

El exclusivo es el que se usa casi siempre, y la razón es que responde a la pregunta “¿en qué posición empiezo yo?”. Si x[i] es el número de elementos que el bloque i va a escribir, el scan exclusivo te da el desplazamiento exacto donde cada bloque debe escribir, sin solaparse con nadie y sin dejar huecos. Ese es el mecanismo con el que se compacta un array, se construye el índice de una rejilla espacial o se reparte trabajo de tamaño variable.

Convertir entre los dos es trivial: el exclusivo es el inclusivo desplazado una posición a la derecha con el neutro al principio. Con la suma también vale restar el elemento original, pero el desplazamiento funciona con cualquier operador y no pierde precisión.

El algoritmo

Hillis-Steele hace log2(n) pasos. En el paso con desplazamiento d, cada posición i mayor o igual que d se combina con la posición i - d.

paso d=1:  y[i] = x[i-1] op x[i]              (rangos de 2)
paso d=2:  y[i] = y[i-2] op y[i]              (rangos de 4)
paso d=4:  y[i] = y[i-4] op y[i]              (rangos de 8)

Después de k pasos, cada posición contiene la combinación de los 2^k elementos que la preceden, ella incluida. Con log2(n) pasos, la última posición cubre el array entero.

La implementación ingenua en un solo array tiene la carrera del segundo tipo que ya vimos: alguien escribe buf[i] mientras otro todavía lee buf[i - d]. Se puede arreglar con dos barreras por paso, pero hay una solución mejor: doble buffer. Se lee siempre de una mitad y se escribe siempre en la otra, y se intercambian. Una sola barrera por paso y ninguna carrera posible.

const TAM: u32 = 256u;
const CAP: u32 = 2u * TAM;
const NEUTRO: f32 = 0.0;
fn comb(a: f32, b: f32) -> f32 { return a + b; }

struct Params { conteo: u32 };
@group(0) @binding(0) var<uniform> params: Params;
@group(0) @binding(1) var<storage, read>       entrada: array<f32>;
@group(0) @binding(2) var<storage, read_write> salida:  array<f32>;

var<workgroup> buf: array<f32, CAP>;

@compute @workgroup_size(TAM)
fn scan(@builtin(global_invocation_id) gid: vec3u,
        @builtin(local_invocation_index) li: u32) {

  let i = gid.x;

  // Doble buffer por desplazamiento: dos mitades del mismo array.
  var lee:    u32 = 0u;
  var escribe: u32 = TAM;

  var v = NEUTRO;
  if (i < params.conteo) { v = entrada[i]; }
  buf[lee + li] = v;
  workgroupBarrier();

  for (var d: u32 = 1u; d < TAM; d <<= 1u) {
    var r = buf[lee + li];
    if (li >= d) { r = comb(buf[lee + li - d], r); }
    buf[escribe + li] = r;
    workgroupBarrier();
    let t = lee; lee = escribe; escribe = t;
  }

  // Inclusivo en buf[lee + li]. Exclusivo, desplazando uno.
  if (i < params.conteo) {
    var e = NEUTRO;
    if (li > 0u) { e = buf[lee + li - 1u]; }
    salida[i] = e;
  }
}

El orden de los argumentos de comb importa. Se escribe comb(anterior, actual), no al revés: con la suma da igual, pero con un operador asociativo no conmutativo —concatenar, componer transformaciones, tomar el último no nulo— el orden es el resultado.

El relleno neutro vuelve a ser lo que hace correcto el kernel para cualquier conteo. El bucle recorre siempre TAM posiciones y el sobrante aporta el neutro, que no altera ninguna combinación. Y ninguna invocación sale antes de tiempo, así que todas llegan a todas las barreras.

El coste, y por qué no es óptimo

En el paso d, las invocaciones que hacen una combinación son las que cumplen li >= d, o sea n - d de ellas. El trabajo total es la suma sobre d = 1, 2, 4, ..., n/2:

W = (n-1) + (n-2) + (n-4) + ... + (n - n/2)
  = n·log2(n) - (n - 1)

Para n = 256: 256·8 − 255 = 1793 operaciones. El bucle secuencial en CPU hace 255. Hillis-Steele gasta siete veces más operaciones que el algoritmo secuencial, y ese factor crece como log2(n).

La profundidad, en cambio, es log2(n) = 8 pasos, contra los 255 del bucle secuencial. Y son 8 barreras, una por paso.

Un algoritmo paralelo cuyo trabajo total supera al del mejor secuencial se llama no óptimo en trabajo. No es un defecto fatal: es un intercambio. Si tienes 256 invocaciones ociosas y el cuello de botella es la latencia, gastar siete veces más operaciones para acabar 32 veces antes es un negocio excelente. Se convierte en un mal negocio cuando el kernel está limitado por ancho de banda o por energía, o cuando n es grande, porque log2(n) crece.

La alternativa óptima en trabajo es el scan de Blelloch, que hace 2n operaciones a cambio del doble de profundidad. Cuál gana depende del tamaño del bloque y del hardware, y la respuesta honesta es que hay que medirlo.

Un apunte de rendimiento que se olvida: Hillis-Steele con doble buffer consume el doble de memoria compartida. Con f32 y bloques de 256 son 2048 bytes, irrelevante. Con bloques de 1024 y vec4f serían 32768 bytes, que ya no caben en el límite de 16384.

Cuándo Hillis-Steele es la respuesta correcta

Tres situaciones, y son frecuentes.

Bloques pequeños. Con 64 o 128 elementos, log2(n) vale 6 o 7 y el exceso de trabajo es modesto. La simplicidad del código, que no tiene la coreografía de índices de Blelloch, se paga sola en tiempo de depuración.

Scan dentro de un grupo SIMD. Si el bloque cabe en un warp, las barreras se pueden omitir y la versión de Hillis-Steele se convierte en una secuencia de instrucciones sin sincronización. Con la extensión de subgrupos, disponible en Chromium desde principios de 2025, ni siquiera hay que escribirla: subgroupInclusiveAdd la hace el hardware. Como la extensión no está en los tres motores, el camino portable sigue siendo necesario.

Operadores caros de memoria pero baratos de cálculo. Si el operador es una suma de flotantes, el trabajo extra son instrucciones aritméticas que la GPU tiene de sobra, mientras que el tiempo se lo lleva la memoria compartida. En ese régimen, los dos algoritmos rinden casi igual y gana el simple.

El caso donde claramente no es la respuesta es el scan de un array grande, porque ahí el bloque conviene que sea lo mayor posible para minimizar el número de bloques, y con bloques grandes el factor log2(n) pesa. Ese caso se resuelve con la descomposición en varios bloques.