Paralelismo en Julia: Threads, @distributed y GPU computing con CUDA.jl

Python tiene el GIL, que limita el paralelismo real en threads para código que no llame a extensiones C. Julia no tiene ese problema: los threads de Julia ejecutan código Julia en paralelo de verdad, con acceso compartido a la misma memoria. Eso simplifica muchos patrones pero también significa que tienes que pensar en condiciones de carrera.

Threads: paralelismo multihilo

Para usar threads hay que lanzar Julia con el número de threads deseado. Desde la terminal: julia --threads 8, o estableciendo la variable de entorno JULIA_NUM_THREADS=8.

using Base.Threads

# Ver cuántos threads hay disponibles
println("Threads: ", nthreads())

# Paralelizar un bucle con @threads
function suma_paralela(n)
    resultados = zeros(nthreads())
    @threads for i in 1:n
        tid = threadid()
        resultados[tid] += i
    end
    return sum(resultados)
end

suma_paralela(10_000_000)

# Alternativa: @threads con reducción explícita
function suma_simple(arr)
    total = Atomic{Float64}(0.0)
    @threads for x in arr
        atomic_add!(total, x)
    end
    return total[]
end

El patrón de dividir los resultados por thread (resultados[tid]) evita la necesidad de bloqueos al escribir. Cada thread escribe en su propio slot y al final sumas todos. Es más eficiente que usar un mutex o un Atomic para acumulaciones.

@spawn: tareas asíncronas

Cuando quieres lanzar una tarea y recoger el resultado más tarde, @spawn es más flexible que @threads:

# Lanzar tareas asíncronas
t1 = @spawn calcular_parte_1(datos)
t2 = @spawn calcular_parte_2(datos)

# Recoger resultados (bloquea hasta que terminan)
r1 = fetch(t1)
r2 = fetch(t2)

resultado = r1 + r2

# Ejemplo más completo: dividir array en chunks
function procesar_paralelo(datos, n_chunks)
    chunk_size = div(length(datos), n_chunks)
    tareas = [@spawn sum(datos[i:i+chunk_size-1])
              for i in 1:chunk_size:length(datos)]
    return sum(fetch.(tareas))
end

Distributed: multiproceso

Para escalar a varios procesos (o varias máquinas), Distributed.jl de la librería estándar añade workers que no comparten memoria. Cada worker es un proceso Julia independiente.

using Distributed

# Añadir workers locales
addprocs(4)   # 4 workers adicionales
println("Workers: ", workers())

# @everywhere ejecuta código en todos los workers
@everywhere function f_costosa(x)
    sum(sin(i) for i in 1:x)
end

# @distributed reduce un bucle distribuido
resultado = @distributed (+) for i in 1:1_000_000
    f_costosa(i % 100 + 1)
end

# pmap: map paralelo distribuido
datos = 1:100
resultados = pmap(f_costosa, datos)

# Workers en máquinas remotas (via SSH)
# addprocs([("servidor1", 4), ("servidor2", 4)])

CUDA.jl: computing en GPU NVIDIA

Para aprovechar GPUs NVIDIA, CUDA.jl proporciona acceso a arrays en GPU con la misma API que los arrays de CPU:

using CUDA

# Verificar disponibilidad
CUDA.has_cuda() && println("GPU: ", CUDA.name(CUDA.device()))

# Crear array en GPU
v_cpu = rand(Float32, 1_000_000)
v_gpu = cu(v_cpu)           # mover a GPU

# Las operaciones estándar funcionan en GPU
resultado = sum(v_gpu)
producto  = v_gpu .* 2.0f0

# Multiplicación de matrices en GPU (usa cuBLAS)
A_gpu = CUDA.rand(Float32, 1000, 1000)
B_gpu = CUDA.rand(Float32, 1000, 1000)
C_gpu = A_gpu * B_gpu   # cuBLAS DGEMM en GPU

# Kernels personalizados en Julia puro
function kernel_suma!(resultado, a, b)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if i <= length(resultado)
        resultado[i] = a[i] + b[i]
    end
    return
end

n = 1024
a = CUDA.rand(Float32, n)
b = CUDA.rand(Float32, n)
c = CUDA.zeros(Float32, n)
@cuda threads=256 blocks=4 kernel_suma!(c, a, b)

Metal.jl: GPU Apple Silicon

Para Macs con chips M1/M2/M3, Metal.jl proporciona acceso al GPU con una API similar a CUDA.jl:

using Metal

# API prácticamente idéntica a CUDA.jl
v_metal = MtlArray(rand(Float32, 1000))
resultado = sum(v_metal)

Cuándo usar cada enfoque

El paralelismo en Julia se integra bien con el resto del ecosistema. Flux.jl usa CUDA.jl por debajo para entrenamiento en GPU. Las operaciones de álgebra lineal ya son paralelas a nivel BLAS sin que tengas que hacer nada. Y para distribución de trabajo a gran escala, el modelo de Distributed.jl es más simple que frameworks como Spark, aunque para datos muy grandes puede valer la pena mirar Scala con Apache Spark.

Imagen: Pexels / Valentine Tanasovich