· Grilla bidimensional M de tamaño total N = L×L.
· Cada celda contiene un spin, M(i,j) ∈ (0..q].
· Los valores de q van desde 2 hasta algunos cientos.
· Para q=2, es el Modelo de Ising.

Energía

La energía H está dada por el Hamilitoniano:

Esto es:
· Para cada celda i,j acumular los primeros vecinos.
· Se acumula el delta de Kronecker δ(i,j) = i==j.
· Esto sobre todas las casillas.

Dinámica probabilística, usando el algoritmo Metrópolis.

Un Paso de Monte Carlo (MCS) está dado por el mapeo a todas las casillas i,j de:

• Sortear un nuevo spin q' ≠ M(i,j).
• Calcular la energía local vieja y nueva: H(i,j), H'(i,j).
• Calcular el ΔE = H'(i,j) - H(i,j).
• Si la energía disminuye, se toma el cambio.
• Si la energía aumenta:
  • Se toma el cambio de manera probabilística según ΔE y temp.

Notar que:

• H(i,j) ∈ [0,4].
• ΔE ∈ [-4,4].

Ejemplo evolución de la energía: q=24, L=200. Notar como el sistema pasa de totalmente desordenado (máxima energía) a ordenado (mínima energía) de manera muy rápida: transición de fase.

Para q≤4 esta discontinuidad no existe. Da la pauta que es un modelo muy rico a pesar de su simplicidad.

L=2048, q=9.

• Repetir samples=1 veces:
  • Comenzar desde un estado totalmente ordenado M(i,j) = 0.
  • Barrer desde temperatura Tmin=0.7212 hasta Tmax=0.721347 en δt=0.00001 (3):
    • Correr ttran=100000 MCS para lograr equilibrio.
    • Correr tmax=10000 MCS, tomando mediciones cada δT=500 MCS.

Total:

• 1650000 MCS.
• 300 mediciones.
• 0 movimientos de matrices host <-> device (wink)

Por cada celda:

• Un entero sorteado uniformemente en [0,q).
• Un fp32 sorteado uniformemente en [0,1).

Para el ejemplo anterior totaliza Doce Teras de números aleatorios.

>>> 1650000*2048*2048*2.0 / (1<<40)
12.5885009765625

Necesitamos un RNG:

• Rápido.
• Período grande.
• Paralelo.
  • Gran cantidad de secuencias independiente.
  • Estado pequeño.

Procesar las celdas puede tener diversos órdendes.

• Secuencial: grilla vieja, grilla nueva, el orden i,j es irrelevante.
• Typewriter: una grilla, orden i,j secuencial.
• Checkerboard: una grilla, celdas blancas y negras, orden i,j irrelevante dentro del mismo color.

Optamos por checkerboard.

Estado: 64 bits.

Siguiente elemento de secuencia: x{n+1} = (xn × a + cn ) mod b.

Donde:
· a es el multiplicador,
· b la base, fijada a (1<<32) para operar bitwise
· cn el acarreo del módulo anterior.

Cada multiplicador bueno genera un secuencia independiente de RNG, pero los multiplicadores a tienen que satisfacer.

goodmult(a) ≡ prime(a × 2^32 − 1) ∧ prime((a × 2^32 − 2)/2)

Obtener estos goodmult es computacionalmente intensivo: computarlos y almacenarlos.

¿Cuántos RNG independientes aka goodmult(a) se puede tener?

Tensión entre:
· Paralelismo.
· Memoria ocupada por el estado de los RNG y calidad de las secuencias.

Para esta implementación fijamos.

• 512×(512/2) = 131072 hilos, secuencias independientes de RNG.
• 64 bits estado, 32 bits multiplicador: 12 bytes. Total 1.5MB.
• ¿Calidad?

Para trabajar tamaños mayores 512×512 ...

Framing

Un hilo computa múltiples casillas, una por frame.

Adherimos a la práctica "Compute multiple outputs per thread": Volkov y Baxter.

Primeros vecinos divide en dos subconjuntos independientes: blancas y negras.

Dentro de cada subconjunto no hay dependencia de datos.

Se pueden procesar en paralelo.

Conocido como Red-Black Gauss-Seidel.

• Tablero dividido en 4 frames.
• Se marca con t0 todas las casillas que procesa el hilo 0. Una por frame.
• Mapeo bijectivo para dividir en casillas blancas y negras: compaction.
  Un poco más de localidad en lecturas y escrituras.

• Usar registros para mantener el estado y multiplicador del RNG.
  • Aumentó la performance un 30%, suavizó las curvas de performance vs L.
    Era obvio, pero no lo vimos.
• El cálculo de la casilla sur tiene su complejidad (mental, no computacional).
• Para obtener un spin diferente se le suma al actual un número en [1,q-1].
• Probamos muchas variaciones para calcular la suma de los delta de Kronecker.
  • Quedó la versión sencilla.
• expf(-delta_E/temp) se puede tabular para cada temperatura, delta_E tiene 9 valores posibles.
  • (obviamente) No se gana nada.
• La escritura a memoria es probabilísica. Para medir performance la hicimos determinística.
  int change = delta_E<=0 || p<=expf(-delta_E/temp);
write[i*L + j] = (change)*spin_new + (1-change)*spin_old;

Notar que hay barras de error, pero son tan chicas que se ocultan en el punto.

• GTX 280, GTX 470, GTX 480, con nvcc-3.2 -O3 (circa 2010).
  Intel Core 2 Duo E8400@3.0 GHz CPU, con gcc-4.4.5 -O3 -ffast-math -march=native -funroll-loops.
• La versión CPU no tiene dependencias con el tamaño.
• Variaciones de q no afectan la performance.
• Hay alguna explicación para la "panza" que forman las curvas.
• La curva CC 2.0 se logró definiendo bloques de 32×32 en la arquitectura Fermi (1024 hilos por bloque y no 512 como en Fermi).
• El mismo .cu corriendo en una placa Fermi, no tiene diferencias si lo compilamos para CC 1.3 y CC 2.0.
• Barrimos todas los flags de CC 2.0 para caché. Ganó -Xptxas -dlcm=cg -Xptxas -dscm=cg: deshabilitar la cache L1. (ver manual de PTX).
• La máxima diferencia CPU-GPU es 22.8ns vs. 0.147ns = 155x. ¿Pusimos un esfuerzo comparable en CPU-GPU?
• Probamos otros framings, más chicos y más grandes: 512×512 es una buena relación performance-memoria ocupada.

Corrida actual sobre Tesla C2070. Probamos algunas cosas, como por ejemplo compilar con nvcc-4.0 (Open64 backend) vs. nvcc-4.2 (LLVM backend).

CUDA-4.0

Q L   spinflip stddev
9 512 0.422714 0.00286214
9 1024 0.294745 0.00165411
9 2048 0.259109 0.00131451
9 4096 0.244878 0.00108468
9 8192 0.236913 0.000977784
9 16384 0.233761 0.00104113
9 32768 0.238528 0.00123815

CUDA-4.2

9 512 0.423447 0.00246027
9 1024 0.282109 0.000991892
9 2048 0.25659 0.000240049
9 4096 0.243577 0.000213649
9 8192 0.236021 0.000164541
9 16384 0.234637 0.000110928
9 32768 0.239124 0.000582012

· Se comparta casi igual en cuanto a tiempos. Notar que si mejora la desviación estándar en nvcc-4.2.
· Los tiempos son mayores que los reportados para GTX 480 (0.147ns). Las Tesla tienen menor frecuencia para los SM y para la memoria.

Intentamos cambiar todos los unsigned int a int. En scan_warp producía efectos benéficos.

CUDA-4.0

Q L   spinflip stddev
9 512 0.4268 0.0025893
9 1024 0.327349 0.00150298
9 2048 0.290854 0.00140329
9 4096 0.27214 0.00111708
9 8192 0.264136 0.00108918
9 16384 0.260818 0.00106291
9 32768 0.261503 0.000845782

CUDA-4.2

9 512 0.442921 0.00349175
9 1024 0.305273 0.00077226
9 2048 0.297691 0.00022777
9 4096 0.288527 0.000144492
9 8192 0.283241 8.79968e-05
9 16384 0.28131 8.16863e-05
9 32768 0.282043 0.000379256

· Solo empeora.
· Acá si se nota algo de diferencia entre los compiladores.

Versión "de la noche después" de Carlos Bederián.

CUDA-4.0

Q L   spinflip stddev
9 512 0.320011 0.00632055
9 1024 0.255473 0.00358806
9 2048 0.224165 0.00214458
9 4096 0.214321 0.00229368
9 8192 0.210495 0.002451
9 16384 0.208602 0.00248176
9 32768 0.208539 0.00221192

CUDA-4.2

9 512 0.327076 0.00926237
9 1024 0.226146 0.00335537
9 2048 0.221503 0.000969065
9 4096 0.213315 0.000593172
9 8192 0.210459 0.000496381
9 16384 0.209569 0.000602072
9 32768 0.21212 0.000555135

· 15% de ganancia, el mejor código hasta ahora.
· Nota que sigue siendo consistente la menor stddev de nvcc-4.2.
· No hay diferencia entre los compiladores.

Algunos cambios menores.

CUDA-4.2 (int)

9 512 0.352479 0.00793455
9 1024 0.243674 0.00314107
9 2048 0.233309 0.0018368
9 4096 0.221555 0.000589989
9 8192 0.217872 0.00119905
9 16384 0.215092 0.000878846
9 32768 0.216102 0.000658448

CUDA-4.2 (-Xptxas -dlcm=cg -Xptxas -dscm=cg)

9 512 0.342793 0.0105217
9 1024 0.22178 0.00296488
9 2048 0.215611 0.000488704
9 4096 0.209739 0.000243782
9 8192 0.206776 0.000155276
9 16384 0.20567 0.000127431
9 32768 0.205315 0.000125488

· Nuevamente cambiar a int no mejora.
· 2% de mejora si deshabilitamos la cache L1.

Finalmente comprobamos que una versión sin compactar produce una menor utilización del BW de memoria.

unsigned int

9 512 0.44261 0.00905423
9 1024 0.334046 0.0046096
9 2048 0.291685 0.00199911
9 4096 0.272963 0.00220491
9 8192 0.265586 0.00233615
9 16384 0.264296 0.00232231
9 32768 0.265268 0.0021042

int

9 512 0.441877 0.00915216
9 1024 0.334246 0.00464585
9 2048 0.29183 0.00292973
9 4096 0.27297 0.00219631
9 8192 0.265794 0.00242677
9 16384 0.26442 0.00271465
9 32768 0.265333 0.00218667

· El código es nuevo y no fue verificado aun.
· Pierde un 10% de performance.

(si, esto es una invitación a un trabajo final)

• Revisar otros generadores de números aleatorios.
  · M. Manssen, M. Weigel, A. Hartmann, Random number generators for massively parallel simulations on GPU, 2012.
  · Alguno es stateless!
• Un byte por int es como poco. 1/4 del BW de memoria real.
  · Empaquetar 4 spines por uint.
  · Usar SIMD Video Instructions de CC 2.0.
  · ¿4x a los ya ~100x de la versión secuencial?
• Hacer una versión CPU SSE, aunque el problema sea memory BW bound.
• Aplicar las técnicas de Volkov:
  V. Volkov. Programming inverse memory hierarchy: case of stencils on GPUs, International Conference on Parallel Computational Fluid Dynamics 2010 (ParCFD 2010).
• Hacer real-time rendering:
  R. Farber, Chapter 9: Mixing CUDA and Rendering, CUDA Application Design and Development, Morgan Kaufmann, 2011.

• E. Ferrero, J. P. De Francesco, N. Wolovick, S. Cannas, q-state Potts model metastability study using optimized GPU-based Monte Carlo algorithms, arXiv, on Computer Physics Communications, 183 (2012) 1578–1587.
• E. Ferrero, Introducción al Cálculo Numérico en Procesadores Gráficos, Clase 9, CAB, 2012.

• No hay clase que viene!
• Gracias a todos los que llegaron hasta este punto.
• Seguimos en la lista de correo.