Bienvenidos a este apéndice escrito con el fin de enseñar a crear números aleatorios con CUDA C.
¿Por qué no usamos cuRAND - CUDA C para la segunda parte de las notas?
Para esa pregunta hay dos respuestas. La primera es que uno de los objetivos generales de estas notas eran enseñar al lector tanto CUDA C como pyCUDA. Los autores pensamos que la parte 2 era un buen punto de inicio para pyCUDA; y es ahí donde entra la segunda respuesta. cuRAND - CUDA C tiene una sintáxis mucho más complicada que cuRAND en pyCUDA, por lo que no quisimos complicarnos (y complicarles) la vida, cuando pyCUDA puede ser muy accesible.
Así que empecemos ahora con esta parte, que es como una continuación de la parte 1 de las notas, para aquellos curiosos que quieran aprender a generar números aleatorios con CUDA C.
Como ya seguro sabrán, las computadoras no pueden en realidad generar números aleatorios per se, sino que generan números pseudoaleatorios o cuasialeatorios. Esto debido al determinismo de la computadora misma.
Sin embargo los matemáticos y computólogos han trabajado durante décadas para poder crear algoritmos que creen números lo suficientemente aleatorios. Hasta hoy, el algoritmo preferido por todos aquellos que trabajan en estadística es el Mersenne Twister creado en el Japón a finales del milenio pasado. Para un poco más de información sobre esto revisa las referencias de este notebook.
CUDA ha decidido utilizar en su librería cuRAND otros generadores tales como el XORWOW o el SOBOL. Pero no nos metamos en eso. Hasta hace unas versiones, el generador MT no estaba disponible en device API. Afortunadamente han corregido eso.
cuRAND fue creada para la utilización de CUDA en otros ambientes tales como los financieros, matemáticos, físicos donde la estadística puede ser pan de cada día.
Para todos aquellos ingenuos mal acostumbrados (como nosotros en un principio) que piensan que únicamente basta con poner en el kernel algo así como
__global__ void Aleatorio(int n, float *d_A)
{
idx = blockDim.x*blockIdx.x + threadIdx.x ;
d_A[idx] = rand() ;
}
Pues están equívocados. El proceso es más complicado que eso y hay algunos detalles extras que poner.
Uno de los primeros detalles a comentar es a la hora de compilar. al usar la librería cuRAND uno tiene escribir la bandera -lcurand en la línea nvcc.
In [ ]:
nvcc mi_primer_programa.cu -lcurand
De la misma manera, en el programa se tiene que incluir la librería escribiendo
#include <curand.h>
Una de las diferencias esenciales entre usar Host API y Device API es el dinamismo del problema. Si de entrada uno sabe la cantidad de números aleatorios y no desea controlar esta cantidad, Host API puede llegar a ser más conveniente y práctico. Sin embargo, si uno no sabe la cantidad de número aleatorios que necesitara y desea tener un mayor control sobre el programa, Device API, aunque más tedioso y complejo, traerá mayores ventajas.
El procedimiento es en realidad bastante sencillo.
curandRandGenerator_t.GPU la memoria que será destinada a los números aleatorios con cudaMalloc().curandCreateGenerator() y curandSetPseudoRandomGeneratorSeed().GPU mediante curandGenerateUniform().En este último paso hay distintas opciones para generar números aleatorios en distintas distribuciones. En nuestra lista utilizamos una distribución uniforme, pero también las hay normal bajo curandGenerateNormal() o log-normal con curandGenerateLogNormal().
Cada una de estas distribuciones tiene la opción de dar un número con precisión doble agregando Double al final del nombre de la función, como en curandGenerateLogNormalDouble().
También hay la opción de dar dúplas de números aleatorios (bastante útil) agregando un 2 al final del nombre de la función, como en curandGenerateNormal2().
cudaFree() y free(), en cuRAND tenemos curandDestroyGenerator().Mostremos un ejemplo sencillo para irnos acostumbrando a esto.
// Este programa utiliza cuRAND para generar 10 números aleatorios
#include<stdio.h>
#include<stdlib.h>
#include<cuda.h>
#include<curand.h>
int main ( int argc, char ∗argv[] ) {
int n = 10 ;
curandGenerator_t gen ; // Creamos la variable gen que será nuestro generador
float ∗devData, ∗hostData ;
hostData = (float∗)calloc(n, sizeof(float) ) ; // Alocación de 10 floats en el CPU
cudaMalloc( (void ∗∗) &devData, n∗sizeof(float) ) ; // Alocación de la memoria en el GPU
curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_MTGP32) ; // Creación de un generador MT
curandSetPseudoRandomGeneratorSeed(gen, 1234ULL) ; // Ajustamos la semilla o seed
curandGenerateUniform(gen, devData, n) ; // Generamos los números aleatorios en el GPU
cudaMemcpy( hostData, devData, n∗sizeof(float), cudaMemcpyDeviceToHost ) ;
// Unas líneas para mostrar los resultados
printf( ”Distribucion uniforme entre 0. y 1.:\n” ) ;
for( int i = 0 ; i < n ; i ++)
{
printf(”%1.4f\n”, hostData[i]) ;
}
printf( ”\n” ) ;
curandDestroyGenerator(gen) ;
cudaFree(devData) ;
free(hostData) ;
return 0 ;
}
Este código lo colocamos en un archivo llamado host_api.cu
In [ ]:
nvcc host_api.cu -lcurand -o host_api
In [ ]:
./host_api
Como mencionamos anteriormente, Device API nos da el control sobre todos los generadores y "números aleatorios" que creamos en el GPU. Pero mayores poderes conllevan mayores responsabilidades, por lo que el código se vuelve a primera vista un poco más complicado, sin embargo la idea de fondo es bastante sencilla.
A diferencia de Host API, en este caso habrá que escribir un kernel el cual se encarge de iniciar un generador distinto en cada thread que vayamos a ocupar. Este kernel tiene la forma siguiente:
__global__ void Init(curandState * state)
{
int idx = blockDim.x*blockIdx.x + threadIdx.x ;
curandinit( 1234, idx, 0, state[idx]) ;
}
Veamos. Los generadores que había en nuestro ejemplo Host API han sido remplazados por las variables tipo State. En este ejemplo hemos usado un generador de números pseudo-aleatorios (el predefinido para Device API es XORWOW, para MT lo haremos más adelante) . La variable se declara en el host como curandState. Claro hay también otro tipo de generadores, que pueden ser consultados en la documentación.
La segunda cosa que salta a la vista es la función curandinit() la cual coloca en cada thread un generador. Los cuatro argumentos de curandinit() se refieren a:
semilla o seed del generador.thread no generará la misma cadena de números aleatorios.offset, el cual para fines prácticos lo anularemos.threads.Una vez que se hayan iniciado los generadores en el GPU, ¡es tiempo de poner manos a la obra!
Ahora un kernel muy sencillo con el que únicamente hacemos un Arreglo de números aleatorios.
__global__ void Arreglo_Alea(curandState * state, float * d_randArray)
{
int idx = blockDim.x*blockIdx.x + threadIdx.x ;
curandState localState = state[idx] ;
d_randArray[idx] = curand_uniform(&localState) ;
state[idx] = localState ;
}
La función curand_uniform() es básicamente la misma que curandGenerateUniform() en el Host API, por lo que no tendría que haber ningún problema.
Creamos la variable localState para fines de eficiencias, jugando un poco con el tiempo de acceso entre la función curand_uniform() y un estado local.
La última línea nos asegura que si en un momento dado se quiere lanzar el kernel más de una vez, este de distintos números aleatorios cada vez que es lanzado.
Veamos el programa completo.
#include <stdio.h>
#include <stdlib.h>
#include <curand.h>
#include <curand_kernel.h>
#define T_BLOQUE 16
#define T_GRID 16
__global__ void Init(curandState * state)
{
int idx = blockDim.x*blockIdx.x + threadIdx.x ;
curandinit( 1234, idx, 0, state[idx]) ;
}
__global__ void Arreglo_Alea(curandState * state, float * d_randArray)
{
int idx = blockDim.x*blockIdx.x + threadIdx.x ;
curandState localState = state[idx] ;
d_randArray[idx] = curand_uniform(&localState) ;
state[idx] = localState ;
}
int main(int argc, char *argv[] )
{
float *d_Resultados, h_Resultados ;
curandState *d_Estados ;
h_Resultados = (float *) calloc( T_BLOQUE * T_GRID, sizeof(float) ) ;
cudaMalloc( (void **)&d_Resultados, T_BLOQUE*T_GRID*sizeof(float) ) ;
cudaMalloc( (void **)&d_Estados, T_BLOQUE*T_GRID*sizeof(curandState) ) ;
cudaMemset( d_Resultados, 0., T_GRID*T_BLOQUE*sizeof(float) ) ;
dim3 dimBlock(T_BLOQUE, 1, 1) ;
dim3 dimGrid(T_GRID, 1, 1) ;
Init<<<dimGrid, dimBlock>>>(d_Estados) ;
Arreglo_Alea<<<dimGrid, dimBlock>>>(d_Estados, d_Resultados) ;
cudaMemcpy(h_Resultados, d_Resultados, T_BLOQUE*T_GRID*sizeof(float), cudaMemcpyDeviceToHost) ;
for (int i = 0; i < T_BLOQUE*T_GRID ; i ++)
{
printf("%1.4f \n", h_Resultados[i]) ;
{
printf("\n") ;
cudaFree(d_Estados) ;
cudaFree(d_Resultados) ;
Free(h_Resultados) ;
return 0;
}
Un detalle muy importante es que a la hora de usar Device API hay que incluir una un paquete más al inicio de nuestro código. Este es
#include <curand_kernel.h>
Fuera de eso, también es importante notar el trato que se le da a la variable d_Estados.
curandState.GPU. En Device API no hay una función como curandDestroyGenerator() como en Host API.Aquí hemos incluido también una nueva función nativa de CUDA C: Memset(). Esto nos ahorra varias líneas a la hora de querer fijar un arreglo con valores fijos tales como $0$ o $1$.
Pararemos por ahora los detalles generales de Host API y Device API esperando que el lector se lleve una buena impresión sobre como utilizar cuRAND. A lo largo de los siguientes notebooks se ahondará en todos estos temas para que se pueda programar con más seguridad y fluidez.
Ahora pasaremos a un ejemplo concreto: el crear generadores Mersenne Twister con Device API. Esto debido a que hasta el momento el generador predefinido es XORWOW, y a diferencia de Host API, el crear un generador MT no es tan sencillo.
In [ ]:
nvcc device_api.cu -lcurand -o device_api
In [ ]:
./device_api
Para generar pseudo-números aleatorios con el algoritmo MT tendremos que importar algunas librerías extras, las cuales están señaladas en el código escrito a continuación. Se importan no sólo algunas funciones extras, sino también parámetros ya computados.
Una vez hecho esta importación, hay que escribir algunas líneas extras las cuales no daremos detalle tan extensivo, con el objetivo de que no se vuelva cansado y tedioso.
En primer lugar hay que alojar la memoria de los estados curandState y los parámetros correspondientes a MT. Una vez hecho esto, modificamos los parámetros importados. Para más información ir a la documentación de Nvidia sobre Mersenne Twister en las referencias.
Nota: Si las cosas parecían complicarse aún más, en este caso tenemos la desventaja que el máximo número de generadores activos en GPU's tipo Tesla es de 256. Para el caso de GPU's tipo Fermi, 1024 generadores pueden ser activados sin -en principio- tener mayor problema.
Hemos escrito más comentarios que lo acostumbrado en el código a continuación para una mejor comprensión del lector.
#include <stdio.h>
#include <stdlib.h>
#include <cuda.h>
#include <curand_kernel.h>
/* incluimos las funciones auxiliares del host de MTGP */
#include <curand_mtgp32_host.h>
/* incluimos los parámetros pre-computados de MTGP */
#include <curand_mtgp32dc_p_11213.h>
#define T_GRID 64
#define T_BLOQUE 256
__global__ void generate_kernel(curandStateMtgp32 *state, int *d_randArray)
{
int idx = blockDim.x*blockIdx.x + threadIdx.x ;
d_randArray[idx] = curand(&state[idx]) ;
}
int main(int argc, char *argv[])
{
curandStateMtgp32 *d_MTGPStates ;
mtgp32_kernel_params *d_KernelParams ;
int *d_Results, *h_Results ;
/* Allojamos espacio para los resultados en el CPU*/
h_Results = (int *)calloc( T_GRID*T_BLOQUE, sizeof(int)) ;
/* Alojamos espsacio para los resultados en el GPU */
cudaMalloc((void **)&d_Results, T_GRID * T_BLOQUE * sizeof(int)) ;
/* Fijamos los resultados como 0 */
cudaMemset(d_Results, 0, T_GRID * T_BLOQUE * sizeof(int)) ;
/* Alojamos espacio para los estados en el GPU */
cudaMalloc((void **)&d_MTGPStates, T_GRID * sizeof(curandStateMtgp32)) ;
/* Hasta ahora todo está igual, ahora vienen los ajustes para MTGP32 */
/* Alojamos espacio para los parámetros del kernel MTGP32 */
cudaMalloc((void**)&d_KernelParams, sizeof(mtgp32_kernel_params)) ;
/* Redefinimos los parámetros predefinios al formato del kernel */
/* y copiamos los parámetros del kernel a la memoria del GPU */
curandMakeMTGP32Constants(mtgp32dc_params_fast_11213, d_KernelParams);
/* E iniciamos un estado por bloque */
curandMakeMTGP32KernelState(d_MTGPStates, mtgp32dc_params_fast_11213, d_KernelParams, T_GRID, 1234);
/* Ajustes completos */
/* Lanzamos el kernel para generar números pseudo-aleatorios */
dim3 dimGrid(T_GRID, 1, 1) ;
dim3 dimBlock(T_BLOQUE, 1, 1) ;
generate_kernel<<<dimGrid, dimBlock>>>(d_MTGPStates, d_Results) ;
/* Copiamos la memoria de vuelta al CPU */
cudaMemcpy(h_Results, d_Results, T_GRID * T_BLOQUE * sizeof(int), cudaMemcpyDeviceToHost) ;
/* Mostramos los resultados */
for(int i = 0; i < T_GRID*T_BLOQUE ; i++) {
printf("%1.f \n", hostResults[i]) ;
}
printf("\n")
/* Limpiamos */
cudaFree(d_MTGPStates) ;
cudaFree(d_Results) ;
free(h_Results) ;
}
cuRAND no es una librería tan sencilla para generar números aleatorios como se puede observar en otros lenguajes de programación. Sin embargo detrás de todas estas líneas de código se encuentran grandes ventajas como el control simultáneo de los generadores de números pseudo-aleatorios, además de la gran eficiencia de los códigos en parelelo.
A lo largo de los siguientes notebooks iremos sacando todo el jugo que esta librería puede darnos. Sin embargo primero tenremos que pasar por algunos ejemplos clásicos como caminantes aleatorios o distribuciones de probabilidad. Así que vayamos tomando vuelo porque esto apenas empieza.
In [ ]: