Por qué ordenar en GPU es otro problema

La ordenación es el problema mejor estudiado de la informática y su solución práctica en CPU lleva cincuenta años cerrada: un quicksort con mediana de tres y corte a inserción, o un timsort si el orden previo importa. Ninguno de los dos se puede ejecutar en un compute shader, y no por rendimiento: por imposibilidad estructural. WGSL prohíbe la recursión, el tamaño de las particiones depende de los datos, y la divergencia de la comparación destruye la ejecución en bloqueo de paso. Ordenar en GPU exige una familia de algoritmos completamente distinta, y entender por qué es entender el modelo de ejecución de golpe.

Lo que el modelo de ejecución prohíbe

No hay recursión. WGSL no la permite, ni directa ni mutua: la especificación exige que el grafo de llamadas sea acíclico. No es una limitación del compilador que se pueda rodear con una pila explícita sin coste, aunque eso sea lo que se hace en un recorrido de BVH. El motivo es de hardware: no hay pila por invocación, los registros se asignan estáticamente al compilar, y una llamada recursiva de profundidad desconocida haría imposible reservarlos. Un quicksort es recursión pura.

Los tamaños dependen de los datos. Tras una partición, el subarray izquierdo tiene un tamaño que no se conoce hasta ejecutar. En una GPU el reparto de trabajo se decide antes del dispatch, y no hay forma de que una invocación decida generar más invocaciones. Existe el dispatch indirecto, que permite que la GPU fije el tamaño del dispatch siguiente, pero eso convierte cada nivel de recursión en un dispatch, y un quicksort de un millón de elementos tendría veinte niveles con tamaños desconocidos y un desequilibrio impredecible.

La comparación diverge. El corazón de un quicksort es “si es menor que el pivote, va a la izquierda”. Dentro de un warp, unos carriles toman una rama y otros la otra, así que el hardware ejecuta las dos y enmascara. Peor: el destino de cada elemento depende de cuántos elementos anteriores fueron a la izquierda, que es un scan. La versión paralela de la partición de un quicksort existe y funciona, pero cuesta un scan por nivel y ha dejado de parecerse a un quicksort.

Hay una cuarta razón menos obvia y a menudo decisiva: el patrón de acceso a memoria. Un quicksort lee y escribe posiciones que dependen de los datos, así que las 32 invocaciones de un warp acceden a direcciones dispersas y la memoria se serializa. Los algoritmos que van bien en GPU acceden a direcciones que se pueden calcular a partir del índice de la invocación, no de los valores.

Redes de ordenación

Una red de ordenación es una secuencia fija de comparadores. Un comparador es un par de posiciones y una dirección: mira los dos valores y, si están en el orden equivocado, los intercambia. La secuencia no depende de los datos en absoluto: la red de ocho elementos ejecuta exactamente los mismos comparadores en el mismo orden tanto si la entrada está ordenada como si viene del revés.

Eso es exactamente lo que una GPU quiere. El control de flujo es uniforme, el patrón de memoria es calculable, y los comparadores de un mismo paso son independientes entre sí, así que se pueden ejecutar todos a la vez.

Un paso de red sobre 8 posiciones, con distancia 4:

  0 --+-- 0        cada linea vertical es un comparador
  1 --|+- 1        entre la posicion i y la i+4
  2 --||+ 2
  3 --|||+3
  4 --+||| 4
  5 ---+|| 5
  6 ----+| 6
  7 -----+ 7

Las dos métricas de una red son el tamaño, que es el número total de comparadores y equivale al trabajo, y la profundidad, que es el número de pasos secuenciales. Una red de ordenación de n elementos necesita al menos n·log2(n) comparadores por el argumento de la teoría de la información, y las redes prácticas quedan por encima: la bitónica usa (n/4)·log2(n)·(log2(n)+1), o sea O(n·log²n).

Ese exceso es el precio de la fijeza. Se paga con gusto porque el algoritmo óptimo en trabajo, con O(n log n) comparaciones, no se puede ejecutar en paralelo sin recursión.

El principio del cero y el uno

Verificar una red de ordenación parece exigir probarla con todas las permutaciones posibles, que para 16 elementos son más de veinte billones. Knuth demostró que basta con mucho menos:

Una red de comparadores ordena correctamente todas las secuencias de números si y solo si ordena correctamente todas las secuencias de ceros y unos.

Para n = 16 eso son 65536 casos, que se comprueban en milisegundos. La demostración se apoya en que un comparador conmuta con cualquier función monótona: si la red falla con una secuencia de números, aplicando la función que envía a cero todo lo menor que cierto umbral y a uno el resto se obtiene una secuencia de ceros y unos con la que también falla.

En la práctica esto da un test exhaustivo y barato para cualquier implementación de red:

// Verificacion exhaustiva de una red de n posiciones con el principio 0-1.
function verificarRed(n, aplicarRed) {
  for (let m = 0; m < (1 << n); m++) {
    const a = Array.from({ length: n }, (_, i) => (m >> i) & 1);
    aplicarRed(a);
    for (let i = 1; i < n; i++) {
      if (a[i - 1] > a[i]) return { fallo: m, resultado: a };
    }
  }
  return null;   // la red es correcta para cualquier entrada
}

Ejecutar esto contra una implementación en JavaScript de la misma aritmética de índices que va a llevar el shader es la forma más eficaz de depurar una red bitónica: si los índices están bien en JavaScript, el problema del shader está en otro sitio.

Las dos familias que funcionan, y la tercera opción

La primera familia son las redes de ordenación: bitónica y la de fusión par-impar de Batcher. Trabajan en sitio, no necesitan memoria adicional, no usan atómicas, y sirven para cualquier tipo con un orden total, incluidas claves compuestas con criterios de desempate. Su coste es O(n·log²n) y, sobre todo, un número de dispatches que crece con el cuadrado del logaritmo.

La segunda son los algoritmos basados en conteo: radix sort y counting sort. No comparan: clasifican por el valor de la clave en un número fijo de cubos y usan un scan para calcular dónde va cada uno. Su coste es O(n) por pasada y el número de pasadas depende del número de bits de la clave, no del tamaño del array. A cambio exigen que la clave se pueda convertir en un entero sin signo que preserve el orden, y necesitan memoria auxiliar del tamaño del array.

Y la tercera opción es la que más veces gana: no ordenar. La ordenación en gráficos casi nunca es el objetivo, es un medio. Si lo que quieres es encontrar vecinos, una rejilla espacial te lo da con un conteo y un scan, sin ordenación de propósito general. Si quieres dibujar transparencias en orden, quizá te baste con ordenar por celda de profundidad en 256 cubos, que es un counting sort de una pasada. Si quieres los k mejores de un millón, una reducción con un montículo por workgroup es mucho más barata que ordenar el millón.

Antes de implementar un sort completo conviene contestar dos preguntas: ¿cuántos valores distintos tiene la clave, y hace falta el orden total o basta con agrupar? Si la clave toma pocos valores y solo hace falta agrupar, el problema no era de ordenación.