Nota generada a partir de liga

Notas para contenedor de docker:

Comando de docker para ejecución de la nota de forma local:

nota: cambiar <ruta a mi directorio> por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker.

docker run --gpus all --rm -v <ruta a mi directorio>:/datos --name jupyterlab_nvidia_cuda_c_container -p 8888:8888 -d palmoreck/jupyterlab_nvidia_cuda_c:1.1.0_10.2

password para jupyterlab: qwerty

Detener el contenedor de docker:

docker stop jupyterlab_nvidia_cuda_c_container

Documentación de la imagen de docker palmoreck/jupyterlab_nvidia_cuda_c:1.1.0_10.2 en liga.


Esta nota utiliza métodos vistos en 1.5.Integracion_numerica

Nota: si se desean ejecutar los ejemplos que se presentan a continuación de forma local, es necesario tener una tarjeta gráfica NVIDIA.

Para ejecutar el notebook si están usando google colab:

  1. Elegir en el menú las opciones Entorno de ejecución > Cambiar tipo de entorno de ejecución > Acelerador por Hardware > GPU Runtime > Change > Hardware Accelerator > GPU .

Ver googlecolab github para información sobre google colab.

CUDA C y generalidades de CUDA y GPU

CUDA C

Consiste en extensiones al lenguaje C y en una runtime library. Ver 2.3.CUDA para más información.

Kernel

  • En CUDA C se define una función que se ejecuta en el device* y que se le nombra kernel. El kernel inicia con la sintaxis:
__global__ void mifun(int param){
...
}

y siempre es tipo void (no hay return).

  • El llamado al kernel se realiza desde el host* y con una sintaxis en la que se define el número de threads, llamados CUDA threads (que son distintos a los CPU threads), y bloques, nombrados CUDA blocks, que serán utilizados para la ejecución del kernel. La sintaxis que se utiliza es con <<< >>> y en la primera entrada se coloca el número de CUDA blocks y en la segunda entrada el número de CUDA threads:
__global__ void mifun(int param){
...
}

int main(){
    int par=0;
    mifun<<<N,5>>> (par); //N bloques de 5 threads
}

* Los nombres de host y device hacen referencia a la CPU y a la GPU.

Ejemplos

1) Programa de hello world


In [1]:
%%file hello_world.cu
#include<stdio.h>
__global__ void func(void){
}
int main(void){
    func<<<1,1>>>(); //1 bloque de 1 thread
    printf("Hello world!\n");
return 0;
}


Writing hello_world.cu

Compilación:


In [2]:
%%bash
nvcc --compiler-options -Wall hello_world.cu -o hello_world.out

Ejecución:


In [3]:
%%bash
./hello_world.out


Hello world!

Comentarios:

  • La función main se ejecuta en la CPU.

  • func es un kernel y es ejecutada por los CUDA threads en el device (GPU). Obsérvese que tal función inicia con la sintaxis __global__. En este caso el CUDA thread que fue lanzado no realiza ninguna acción pues el cuerpo del kernel está vacío.

  • El kernel sólo puede tener un return tipo void: __global__ void func por lo que el kernel debe regresar sus resultados a través de sus argumentos.

  • nvcc es un wrapper para el compilador de programas escritos en C. El compilador instalado en el contenedor de docker descrito al inicio de ésta nota es gcc.

  • La extensión del archivo debe ser .cu aunque esto puede modificarse al compilar con nvcc:

$nvcc -x cu hello_world.c -o hello_world.out

  • En ocasiones para tener funcionalidad de un determinado compute capability se especifica la flag de -arch=sm_11 en la línea de nvcc. En este caso se le indica al compilador que compile el programa para un compute capability de $1.1$. Ver liga y liga2 para más sobre esto.

¿Bloques de threads?

Los CUDA threads son divididos los CUDA blocks y éstos se encuentran en un grid. En el lanzamiento del kernel se debe especificar al hardware cuántos CUDA blocks tendrá nuestro grid y cuántos CUDA threads estarán en cada bloque.

2) Programa de hello world 2


In [4]:
%%file hello_world_2.cu
#include<stdio.h>
__global__ void func(void){
    printf("Hello world del bloque %d del thread %d!\n", blockIdx.x, threadIdx.x);
}
int main(void){
    func<<<2,3>>>(); //2 bloques de 3 threads cada uno
    cudaDeviceSynchronize();
    printf("Hola del cpu thread\n");
    return 0;
}


Writing hello_world_2.cu

In [5]:
%%bash 
nvcc --compiler-options -Wall hello_world_2.cu -o hello_world_2.out

In [6]:
%%bash
./hello_world_2.out


Hello world del bloque 0 del thread 0!
Hello world del bloque 0 del thread 1!
Hello world del bloque 0 del thread 2!
Hello world del bloque 1 del thread 0!
Hello world del bloque 1 del thread 1!
Hello world del bloque 1 del thread 2!
Hola del cpu thread

Comentarios:

  • En lo que continúa el nombre thread hace referencia a CUDA thread y el nombre bloque a CUDA block.
  • El llamado a la ejecución del kernel se realizó en el host (CPU) y se lanzaron $2$ bloques (primera posición en la sintaxis <<<>>>), cada uno con $3$ threads.

  • Se utiliza la función cudaDeviceSynchronize para que el cpu-thread espere la finalización de la ejecución del kernel.

  • En el ejemplo anterior, las variables blockIdx y threadIdx hacen referencia a los id's que tienen los bloques y los threads: el id del bloque dentro del grid y el id del thread dentro del bloque. La parte .x de las variables: blockIdx.x y threadIdx.x refieren a la primera coordenada del bloque en el grid y a la primera coordenada del thread en en el bloque.

  • La elección del número de bloques en un grid o el número de threads en un bloque no corresponde a alguna disposición del hardware, esto es, si se lanza un kernel con <<< 1, 3 >>> no implica que la GPU tenga en su hardware un bloque o 3 threads. Asimismo, las coordenadas que se obtienen vía blockIdx o threadIdx son meras abstracciones, no corresponden a algún ordenamiento en el hardware de la GPU.

  • Todos los threads de un bloque ejecutan el kernel por lo que se tienen tantas copias del kernel como número de bloques sean lanzados. Aquí encontramos que en una GPU se tiene el modelo Single Instruction Multiple Threads SIMT.

¿Grid's y bloques 3-dimensionales?

En una GPU podemos definir el grid de bloques y el bloque de threads utilizando el tipo de dato dim3 el cual también es parte de CUDA C:

3) Ejemplo


In [7]:
%%file hello_world_3.cu
#include<stdio.h>
__global__ void func(void){
    printf("Hello world del bloque %d del thread %d!\n", blockIdx.y, threadIdx.z);
}
int main(void){
    dim3 dimGrid(1,2,1); //2 bloques en el grid
    dim3 dimBlock(1,1,3); //3 threads por bloque
    func<<<dimGrid,dimBlock>>>(); 
    cudaDeviceSynchronize();
    printf("Hola del cpu thread\n");
    return 0;
}


Writing hello_world_3.cu

In [8]:
%%bash 
nvcc --compiler-options -Wall hello_world_3.cu -o hello_world_3.out

In [9]:
%%bash
./hello_world_3.out


Hello world del bloque 0 del thread 0!
Hello world del bloque 0 del thread 1!
Hello world del bloque 0 del thread 2!
Hello world del bloque 1 del thread 0!
Hello world del bloque 1 del thread 1!
Hello world del bloque 1 del thread 2!
Hola del cpu thread

4) Ejemplo


In [10]:
%%file thread_idxs.cu
#include<stdio.h>
__global__ void func(void){
    if(threadIdx.x==0 && threadIdx.y==0 && threadIdx.z==0){
        printf("blockIdx.x:%d\n",blockIdx.x);
    }
    printf("thread idx.x:%d\n",threadIdx.x);
    printf("thread idx.y:%d\n",threadIdx.y);
    printf("thread idx.z:%d\n",threadIdx.z);
}
int main(void){
    dim3 dimGrid(1,1,1); //1 bloque en el grid
    dim3 dimBlock(1,3,1); //3 threads por bloque
    func<<<dimGrid,dimBlock>>>(); 
    cudaDeviceSynchronize();
    return 0;
}


Writing thread_idxs.cu

In [11]:
%%bash 
nvcc --compiler-options -Wall thread_idxs.cu -o thread_idxs.out

In [12]:
%%bash
./thread_idxs.out


blockIdx.x:0
thread idx.x:0
thread idx.x:0
thread idx.x:0
thread idx.y:0
thread idx.y:1
thread idx.y:2
thread idx.z:0
thread idx.z:0
thread idx.z:0

5) Ejemplo


In [13]:
%%file block_idxs.cu
#include<stdio.h>
__global__ void func(void){
    printf("blockIdx.x:%d\n",blockIdx.x);
    printf("blockIdx.y:%d\n",blockIdx.y);
    printf("blockIdx.z:%d\n",blockIdx.z);

}
int main(void){
    dim3 dimGrid(1,2,2); //4 bloques en el grid
    dim3 dimBlock(1,1,1); //1 thread por bloque
    func<<<dimGrid,dimBlock>>>(); 
    cudaDeviceSynchronize();
    return 0;
}


Writing block_idxs.cu

In [14]:
%%bash 
nvcc --compiler-options -Wall block_idxs.cu -o block_idxs.out

In [15]:
%%bash
./block_idxs.out


blockIdx.x:0
blockIdx.x:0
blockIdx.x:0
blockIdx.x:0
blockIdx.y:1
blockIdx.y:0
blockIdx.y:0
blockIdx.y:1
blockIdx.z:0
blockIdx.z:0
blockIdx.z:1
blockIdx.z:1

6) Ejemplo

Podemos usar la variable blockDim para cada coordenada x, y o z y obtener la dimensión de los bloques:


In [16]:
%%file block_dims.cu
#include<stdio.h>
__global__ void func(void){
    if(threadIdx.x==0 && threadIdx.y==0 && threadIdx.z==0 && blockIdx.z==1){
    printf("blockDim.x:%d\n",blockDim.x);
    printf("blockDim.y:%d\n",blockDim.y);
    printf("blockDim.z:%d\n",blockDim.z);
    }

}
int main(void){
    dim3 dimGrid(2,2,2); //8 bloques en el grid
    dim3 dimBlock(3,1,2); //6 threads por bloque
    func<<<dimGrid,dimBlock>>>(); 
    cudaDeviceSynchronize();
    return 0;
}


Writing block_dims.cu

In [17]:
%%bash 
nvcc --compiler-options -Wall block_dims.cu -o block_dims.out

In [18]:
%%bash
./block_dims.out


blockDim.x:3
blockDim.x:3
blockDim.x:3
blockDim.x:3
blockDim.y:1
blockDim.y:1
blockDim.y:1
blockDim.y:1
blockDim.z:2
blockDim.z:2
blockDim.z:2
blockDim.z:2

Alojamiento de memoria en el device

Para alojar memoria en el device se utiliza el llamado a cudaMalloc y para transferir datos del host al device o viceversa se llama a lafunción cudaMemcpy con respectivos parámetros como cudaMemcpyHostToDevice o cudaMemcpyDeviceToHost.

Para desalojar memoria del device se utiliza el llamado a cudaFree

7) Programa de suma vectorial

N bloques de 1 thread


In [19]:
%%file suma_vectorial.cu
#include<stdio.h>
#define N 10
__global__ void suma_vect(int *a, int *b, int *c){
    int block_id_x = blockIdx.x;
    if(block_id_x<N) //aquí se asume que el valor de N 
                     //es menor al número máximo de bloques que se pueden lanzar
                     //si fuese mayor, hay que hacer un ajuste
        c[block_id_x] = a[block_id_x]+b[block_id_x];
}
int main(void){
    int a[N], b[N],c[N];
    int *device_a, *device_b, *device_c;
    int i;
    dim3 dimGrid(N,1,1); //N bloques en el grid
    dim3 dimBlock(1,1,1); //1 threads por bloque 
    //alojando en device
    cudaMalloc((void **)&device_a, sizeof(int)*N); 
    cudaMalloc((void **)&device_b, sizeof(int)*N);
    cudaMalloc((void **)&device_c, sizeof(int)*N);
    //llenando los arreglos con datos dummy:
    for(i=0;i<N;i++){
        a[i]=i;
        b[i]=i*i;
    }
    //copiamos arreglos a, b a la GPU
    cudaMemcpy(device_a,a,N*sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(device_b,b,N*sizeof(int), cudaMemcpyHostToDevice);
    //mandamos a llamar a suma_vect:
    suma_vect<<<dimGrid,dimBlock>>>(device_a,device_b,device_c); //N bloques de 1 thread
    cudaDeviceSynchronize();
    //copia del resultado al arreglo c:
    cudaMemcpy(c,device_c,N*sizeof(int),cudaMemcpyDeviceToHost);
    for(i=0;i<N;i++)
        printf("%d+%d = %d\n",a[i],b[i],c[i]);
    cudaFree(device_a);
    cudaFree(device_b);
    cudaFree(device_c);
    return 0;
}


Writing suma_vectorial.cu

In [20]:
%%bash
nvcc --compiler-options -Wall suma_vectorial.cu -o suma_vectorial.out

In [21]:
%%bash
./suma_vectorial.out


0+0 = 0
1+1 = 2
2+4 = 6
3+9 = 12
4+16 = 20
5+25 = 30
6+36 = 42
7+49 = 56
8+64 = 72
9+81 = 90

Comentarios:

  • Obsérvese que se están utilizando apuntadores en la línea:
    int *device_a, *device_b, *device_c;

pero estos apuntadores no apuntan a una dirección de memoria en el device pues aunque NVIDIA añadió el feature de Unified Memory (un espacio de memoria accesible para el host y el device) aquí no se está usando tal feature. Más bien se están utilizando los apuntadores anteriores para apuntar a un struct) de C en el que uno de sus tipos de datos es una dirección de memoria en el device.

  • Obsérvese el uso de (void **) por la definición de la función cudaMalloc.

  • Obsérvese que en el programa anterior se coloca en comentario que se asume que $N$ el número de datos en el arreglo es menor al número de bloques que es posible lanzar. Esto como veremos más adelante es importante considerar pues aunque en un device se pueden lanzar muchos bloques y muchos threads, se tienen límites en el número de éstos que es posible lanzar.

¿Perfilamiento en CUDA?

Al instalar el CUDA toolkit en sus máquinas o bien si utilizan el contenedor de docker (descrito al inicio de la nota) se instala la línea de comando nvprof para perfilamiento. Se puede ejecutar con:


In [22]:
%%bash 
nvprof --normalized-time-unit s ./suma_vectorial.out


0+0 = 0
1+1 = 2
2+4 = 6
3+9 = 12
4+16 = 20
5+25 = 30
6+36 = 42
7+49 = 56
8+64 = 72
9+81 = 90
==342== NVPROF is profiling process 342, command: ./suma_vectorial.out
==342== Profiling application: ./suma_vectorial.out
==342== Profiling result:
            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
                        %         s                   s         s         s
 GPU activities:    36.80  1.47e-06         1  1.47e-06  1.47e-06  1.47e-06  suma_vect(int*, int*, int*)
                    35.20  1.41e-06         2  7.04e-07  5.44e-07  8.64e-07  [CUDA memcpy HtoD]
                    28.00  1.12e-06         1  1.12e-06  1.12e-06  1.12e-06  [CUDA memcpy DtoH]
      API calls:    99.55  0.087959         3  0.029320  4.69e-06  0.087948  cudaMalloc
                     0.25  2.18e-04        97  2.25e-06  1.79e-07  1.36e-04  cuDeviceGetAttribute
                     0.06  5.50e-05         3  1.83e-05  3.40e-06  4.55e-05  cudaFree
                     0.05  4.70e-05         1  4.70e-05  4.70e-05  4.70e-05  cuDeviceTotalMem
                     0.03  2.69e-05         3  8.96e-06  5.04e-06  1.11e-05  cudaMemcpy
                     0.02  1.92e-05         1  1.92e-05  1.92e-05  1.92e-05  cudaLaunchKernel
                     0.02  1.86e-05         1  1.86e-05  1.86e-05  1.86e-05  cuDeviceGetName
                     0.00  4.01e-06         1  4.01e-06  4.01e-06  4.01e-06  cudaDeviceSynchronize
                     0.00  2.45e-06         1  2.45e-06  2.45e-06  2.45e-06  cuDeviceGetPCIBusId
                     0.00  1.65e-06         3  5.49e-07  2.53e-07  1.03e-06  cuDeviceGetCount
                     0.00  9.27e-07         2  4.63e-07  1.84e-07  7.43e-07  cuDeviceGet
                     0.00  2.72e-07         1  2.72e-07  2.72e-07  2.72e-07  cuDeviceGetUuid

Comentarios:

  • Las unidades en las que se reporta son s: second, ms: millisecond, us: microsecond, ns: nanosecond.

  • En la documentación de NVIDIA se menciona que nvprof será reemplazada próximamente por NVIDIA Nsight Compute y NVIDIA Nsight Systems.

En el ejemplo anterior se lanzaron $N$ bloques con $1$ thread cada uno y a continuación se lanza $1$ bloque con $N$ threads.


In [23]:
%%file suma_vectorial_2.cu
#include<stdio.h>
#define N 10
__global__ void suma_vect(int *a, int *b, int *c){
    int thread_id_x = threadIdx.x;
    if(thread_id_x<N) //aquí se asume que el valor de N 
                     //es menor al número máximo de threads que se pueden lanzar
                    //si fuese mayor, hay que hacer un ajuste
        c[thread_id_x] = a[thread_id_x]+b[thread_id_x];
}
int main(void){
    int *device_a, *device_b, *device_c;
    int i;
    dim3 dimGrid(1,1,1); //1 bloques en el grid
    dim3 dimBlock(N,1,1); //N threads por bloque 
    //alojando en device con Unified Memory
    cudaMallocManaged(&device_a, sizeof(int)*N);
    cudaMallocManaged(&device_b, sizeof(int)*N);
    cudaMallocManaged(&device_c, sizeof(int)*N);
    //llenando los arreglos:
    for(i=0;i<N;i++){
        device_a[i]=i;
        device_b[i]=i*i;
    }
    suma_vect<<<dimGrid,dimBlock>>>(device_a,device_b,device_c); //1 bloque con N threads
    cudaDeviceSynchronize();
    for(i=0;i<N;i++)
        printf("%d+%d = %d\n",device_a[i],device_b[i],device_c[i]);
    cudaFree(device_a);
    cudaFree(device_b);
    cudaFree(device_c);
    return 0;
}


Writing suma_vectorial_2.cu

In [24]:
%%bash
nvcc --compiler-options -Wall suma_vectorial_2.cu -o suma_vectorial_2.out

In [25]:
%%bash 
nvprof --normalized-time-unit s ./suma_vectorial_2.out


0+0 = 0
1+1 = 2
2+4 = 6
3+9 = 12
4+16 = 20
5+25 = 30
6+36 = 42
7+49 = 56
8+64 = 72
9+81 = 90
==389== NVPROF is profiling process 389, command: ./suma_vectorial_2.out
==389== Profiling application: ./suma_vectorial_2.out
==389== Profiling result:
            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
                        %         s                   s         s         s
 GPU activities:   100.00  2.05e-06         1  2.05e-06  2.05e-06  2.05e-06  suma_vect(int*, int*, int*)
      API calls:    99.40  0.151464         3  0.050488  1.09e-05  0.151432  cudaMallocManaged
                     0.27  4.08e-04         1  4.08e-04  4.08e-04  4.08e-04  cudaLaunchKernel
                     0.19  2.82e-04        97  2.91e-06  3.09e-07  1.61e-04  cuDeviceGetAttribute
                     0.06  9.38e-05         3  3.13e-05  1.15e-05  5.90e-05  cudaFree
                     0.05  7.66e-05         1  7.66e-05  7.66e-05  7.66e-05  cuDeviceTotalMem
                     0.02  2.95e-05         1  2.95e-05  2.95e-05  2.95e-05  cuDeviceGetName
                     0.01  1.31e-05         1  1.31e-05  1.31e-05  1.31e-05  cudaDeviceSynchronize
                     0.00  2.76e-06         3  9.18e-07  3.62e-07  1.29e-06  cuDeviceGetCount
                     0.00  1.91e-06         1  1.91e-06  1.91e-06  1.91e-06  cuDeviceGetPCIBusId
                     0.00  1.46e-06         2  7.30e-07  4.97e-07  9.63e-07  cuDeviceGet
                     0.00  4.90e-07         1  4.90e-07  4.90e-07  4.90e-07  cuDeviceGetUuid

==389== Unified Memory profiling result:
Device "GeForce GTX 750 (0)"
   Count  Avg Size  Min Size  Max Size  Total Size  Total Time  Name
       1  8.0000KB  8.0000KB  8.0000KB  8.000000KB  1.6960e-06s  Host To Device
       5  25.600KB  4.0000KB  60.000KB  128.0000KB  1.2960e-05s  Device To Host
Total CPU Page faults: 2

Comentarios:

  • El programa anterior utiliza la Unified Memory con la función cudaMallocManaged. La Unified Memory es un feature que se añadió a CUDA desde las arquitecturas de Kepler y Maxwell pero que ha ido mejorando (por ejemplo añadiendo page faulting and migration) en las arquitecturas siguientes a la de Kepler: la arquitectura Pascal y Volta. Por esto en el output anterior de nvprof aparece una sección de page fault.

  • Obsérvese que en el programa anterior se comenta que se asume que $N$ el número de datos en el arreglo es menor al número de threads que es posible lanzar. Esto como veremos más adelante es importante considerar pues aunque en el device se pueden lanzar muchos bloques y muchos threads, se tienen límites en el número de éstos que es posible lanzar.

¿Tenemos que inicializar los datos en la CPU y copiarlos hacia la GPU?

En realidad no tenemos que realizarlo para el ejemplo de suma_vectorial.cu anterior. Por ejemplo:


In [26]:
%%file suma_vectorial_3.cu
#include<stdio.h>
#define N 10
__global__ void init(int *a, int *b){
    int thread_id_x = threadIdx.x;
    a[thread_id_x]=thread_id_x;
    b[thread_id_x]=thread_id_x*thread_id_x;
}

__global__ void suma_vect(int *a, int *b, int *c){
    int thread_id_x = threadIdx.x;
    if(thread_id_x<N) //aquí se asume que el valor de N 
                     //es menor al número máximo de threads que se pueden lanzar
                    //si fuese mayor, hay que hacer un ajuste
        c[thread_id_x] = a[thread_id_x]+b[thread_id_x];
}

int main(void){
    int *device_a, *device_b, *device_c;
    int i;
    dim3 dimGrid(1,1,1); //1 bloques en el grid
    dim3 dimBlock(N,1,1); //N threads por bloque 
    //alojando en device con Unified Memory
    cudaMallocManaged(&device_a, sizeof(int)*N);
    cudaMallocManaged(&device_b, sizeof(int)*N);
    cudaMallocManaged(&device_c, sizeof(int)*N);
    //llenando los arreglos:
    init<<<dimGrid,dimBlock>>>(device_a,device_b);
    cudaDeviceSynchronize();
    //llamando al kernel suma_vect
    suma_vect<<<dimGrid,dimBlock>>>(device_a,device_b,device_c); //1 bloque con N threads
    cudaDeviceSynchronize();
    for(i=0;i<N;i++)
        printf("%d+%d = %d\n",device_a[i],device_b[i],device_c[i]);
    cudaFree(device_a);
    cudaFree(device_b);
    cudaFree(device_c);
    return 0;
}


Writing suma_vectorial_3.cu

In [27]:
%%bash
nvcc --compiler-options -Wall suma_vectorial_3.cu -o suma_vectorial_3.out

In [28]:
%%bash
nvprof --normalized-time-unit s ./suma_vectorial_3.out


0+0 = 0
1+1 = 2
2+4 = 6
3+9 = 12
4+16 = 20
5+25 = 30
6+36 = 42
7+49 = 56
8+64 = 72
9+81 = 90
==436== NVPROF is profiling process 436, command: ./suma_vectorial_3.out
==436== Profiling application: ./suma_vectorial_3.out
==436== Profiling result:
            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
                        %         s                   s         s         s
 GPU activities:    51.28  1.92e-06         1  1.92e-06  1.92e-06  1.92e-06  suma_vect(int*, int*, int*)
                    48.72  1.82e-06         1  1.82e-06  1.82e-06  1.82e-06  init(int*, int*)
      API calls:    99.37  0.149377         3  0.049792  1.17e-05  0.149330  cudaMallocManaged
                     0.28  4.28e-04         2  2.14e-04  1.33e-05  4.15e-04  cudaLaunchKernel
                     0.19  2.93e-04        97  3.02e-06  3.01e-07  1.60e-04  cuDeviceGetAttribute
                     0.07  1.05e-04         3  3.50e-05  1.25e-05  6.47e-05  cudaFree
                     0.05  7.66e-05         1  7.66e-05  7.66e-05  7.66e-05  cuDeviceTotalMem
                     0.02  2.88e-05         1  2.88e-05  2.88e-05  2.88e-05  cuDeviceGetName
                     0.01  1.25e-05         1  1.25e-05  1.25e-05  1.25e-05  cudaDeviceSynchronize
                     0.00  1.99e-06         3  6.63e-07  3.55e-07  1.18e-06  cuDeviceGetCount
                     0.00  1.96e-06         1  1.96e-06  1.96e-06  1.96e-06  cuDeviceGetPCIBusId
                     0.00  1.20e-06         2  6.00e-07  3.18e-07  8.83e-07  cuDeviceGet
                     0.00  5.35e-07         1  5.35e-07  5.35e-07  5.35e-07  cuDeviceGetUuid

==436== Unified Memory profiling result:
Device "GeForce GTX 750 (0)"
   Count  Avg Size  Min Size  Max Size  Total Size  Total Time  Name
       3  21.333KB  4.0000KB  52.000KB  64.00000KB  6.8480e-06s  Device To Host
Total CPU Page faults: 1

Arquitectura de una GPU y límites en número de threads y bloques que podemos lanzar en el kernel

Un device está compuesto por arreglos de streaming multiprocessors SM's (también denotados como MP's) y en cada SM encontramos un número (determinado por la arquitectura del device) de streaming processors SP's que comparten el caché y unidades de control (que están dentro de cada SM):

En el dibujo anterior se muestran las SM's en color rojo y los SP's en morado. Hay dos SM's por cada bloque anaranjado y ocho SP's por cada SM. Así, una GPU es una máquina multicore. Aunque cada SM ejecuta las instrucciones de forma independiente a otra SM, comparten la memoria global.

Los bloques de threads son asignados a cada SM por el CUDA runtime system, puede asignar más de un bloque a una SM pero hay un límite de bloques que pueden ser asignados a cada SM. Ver maximum number of blocks per multiprocessor.

Comentarios:

  • Por ejemplo para el modelo GT200 el máximo número de bloques que podían asignarse a cada SM eran de $8$ bloques. Tal modelo tenía $30$ SM's lo que resultaban en $240$ bloques que en un instante podían asignarse al device para su ejecución simultánea (asignándose en cualquier orden en alguna SM disponible). Por supuesto que un grid podía contener más de $240$ bloques en este modelo y en este caso el CUDA runtime system lleva una lista de bloques que va asignando a cada SM y conforme cada SM terminan la ejecución, nuevos bloques son asignados a tales SM que finalizaron. Para visualizar esta situación, considérese una simplificación de lo anterior en donde se tiene un device con $2$ SM's y con un kernel se han lanzado $6$ bloques. El CUDA runtime system ha asignado $3$ bloques a cada SM, entonces se tiene un dibujo como el siguiente:

  • Los bloques asignados a una SM comparten recursos (por ejemplo memoria) y su ejecución es independiente entre ellos, no es posible sincronizar al "bloque 1" con el "bloque 0". También no es posible sincronizar a los threads de diferentes SM's pero sí es posible sincronizar a los threads dentro de un mismo bloque.

¿Qué otros límites puedo encontrar en mi(s) device(s) de mi sistema?

Para responder lo anterior se puede utilizar el siguiente programa que está basado en liga y cudaDeviceProp Struct Reference:


In [29]:
%%file device_properties.cu

#include<stdio.h>

int main(void){
    cudaDeviceProp properties;
    int count;
    int i;
    cudaGetDeviceCount(&count);
    for(i=0;i<count;i++){
        printf("----------------------\n");
        cudaGetDeviceProperties(&properties, i);
        printf("----device %d ----\n",i); 
        printf("Device Name: %s\n", properties.name);
        printf("Compute capability: %d.%d\n", properties.major, properties.minor);
        printf("Clock rate: %d\n", properties.clockRate);
        printf("Unified memory: %d\n", properties.unifiedAddressing);
        printf(" ---Memory Information for device %d (results on bytes)---\n", i);
        printf("Total global mem: %ld\n", properties.totalGlobalMem); 
        printf("Total constant Mem: %ld\n", properties.totalConstMem);
        printf("Shared memory per thread block: %ld\n", properties.sharedMemPerBlock);
        printf("Shared memory per SM: %ld\n",properties.sharedMemPerMultiprocessor );
        printf(" ---MP Information for device %d ---\n", i);
        printf("SM count: %d\n", properties.multiProcessorCount);
        printf("Threads in warp: %d\n", properties.warpSize);
        printf("Max threads per SM: %d\n", properties.maxThreadsPerMultiProcessor);
        printf("Max warps per SM: %d\n",properties.maxThreadsPerMultiProcessor/properties.warpSize);
        printf("Max threads per block: %d\n", properties.maxThreadsPerBlock);
        printf("Max thread dimensions: (%d, %d, %d)\n", properties.maxThreadsDim[0], properties.maxThreadsDim[1], properties.maxThreadsDim[2]);
        printf("Max grid dimensions: (%d, %d, %d)\n", properties.maxGridSize[0], properties.maxGridSize[1], properties.maxGridSize[2]); 
    }
    return 0;
    
}


Writing device_properties.cu

In [30]:
%%bash

nvcc --compiler-options -Wall device_properties.cu -o device_properties.out

In [20]:
%%bash
./device_properties.out


----------------------
----device 0 ----
Device Name: Tesla K20Xm
Compute capability: 3.5
Clock rate: 732000
Unified memory: 1
 ---Memory Information for device 0 (results on bytes)---
Total global mem: 5977800704
Total constant Mem: 65536
Shared memory per thread block: 49152
Shared memory per SM: 49152
 ---MP Information for device 0 ---
SM count: 14
Threads in warp: 32
Max threads per SM: 2048
Max warps per SM: 64
Max threads per block: 1024
Max thread dimensions: (1024, 1024, 64)
Max grid dimensions: (2147483647, 65535, 65535)
----------------------
----device 1 ----
Device Name: Tesla K20Xm
Compute capability: 3.5
Clock rate: 732000
Unified memory: 1
 ---Memory Information for device 1 (results on bytes)---
Total global mem: 5977800704
Total constant Mem: 65536
Shared memory per thread block: 49152
Shared memory per SM: 49152
 ---MP Information for device 1 ---
SM count: 14
Threads in warp: 32
Max threads per SM: 2048
Max warps per SM: 64
Max threads per block: 1024
Max thread dimensions: (1024, 1024, 64)
Max grid dimensions: (2147483647, 65535, 65535)

Comentarios:

  • También en la documentación oficial de NVIDIA dentro de compute-capabilities se pueden revisar los valores anteriores y muchos más.

  • Obsérvese del output anterior que el sistema tiene dos devices con las mismas capacidades.

  • En un device encontramos diferentes tipos de memoria: global, constante, shared y texture. En esta nota únicamente trabajamos con la memoria global. Tenemos funciones en CUDA para poder comunicar/coordinar a los threads en un bloque por medio de la shared memory. Ver por ejemplo Using Shared Memory in CUDA C/C++ para un pequeño post del $2013$ sobre shared memory.

  • Los bloques de threads que son asignados a una SM son divididos en warps que es la unidad de thread scheduling* que tiene el CUDA run time system. El output anterior indica que son divisiones de $32$ threads: un warp consta de $32$ threads.

* El thread scheduling se puede pensar a la funcionalidad del hardware para seleccionar una instrucción del programa y asginar su ejecución por los threads en un warp (SIMT). Otro ejemplo es tener una instrucción que indica que se debe realizar lectura o escritura, entonces el hardware del device utiliza un warp de threads para tal operación mientras selecciona un warp de threads distinto para seleccionar otra instrucción diferente a la de I/O.

  • El número máximo de threads que pueden iniciarse de forma simultánea o en un instante por SM es de $2048$ o bien $2048/32 = 64$ warps.

  • El output anterior muestra los límites para número de bloques en las tres dimensiones de un grid y el número de threads en las tres dimensiones en un bloque.

  • Un bloque puede tener como máximo $1024$ threads en cualquier configuración: por ejemplo $(1024,1,1), (32,1,32), (4,4,64)$.

  • Por los puntos anteriores si lanzamos bloques de $1024$ threads entonces sólo $2$ bloques pueden residir en una SM en un instante. Con esta configuración alcanzaríamos $1024/32=32$ warps por cada bloque y como lanzamos $2$ bloques alcanzaríamos $64$ warps (que es el máximo de warps por SM que podemos tener en un instante). Otra configuración para alcanzar el máximo número de warps en un instante, es considerar $4$ bloques de $512$ threads pues tendríamos $512/32=16$ warps por bloque y en total serían $16*4$ (warps $\times$ bloques) $=64$ warps. Entre los datos que hay que elegir en los programas de CUDA C se tienen las configuraciones en el número de threads y el número de bloques a lanzar. La idea es alcanzar o rebasar el máximo número de warps en cada SM que soporta nuestro device en un instante.

  • Por ejemplo para el dibujo en el que se asumió que el CUDA runtime system había asignado $3$ bloques a cada SM, se tendría una división de cada bloque en un warp de $32$ threads como sigue:

Grid Configuration Choices?

Los programas de CUDA C tienen la opción de elegir el número de threads y de bloques a ser lanzados. En la referencia "Parallel Computing for Data Science. With Examples in R, C++ and CUDA" de N. Matloff se enlistan algunas consideraciones para elegir tales parámetros:

  • Given that scheduling is done on a warp basis, block size should be a multiple of the warp size (32).

  • One wants to utilize all the SMs. If one sets the block size too large, not all will be used, as a block cannot be split across SM's.

  • ..., barrier synchronization can be done effectively only at the block level. The larger the block, the more the barrier delay, so one might want smaller blocks.

  • On the other hand, if one is using shared memory, this can only be done at the block level, and efficient use may indicate using a larger block.

  • Two threads doing unrelated work, or the same work but with many if/elses, would cause a lot of thread divergence if they were in the same block. In some cases, it may be known in advance which threads will do the "ifs" and which will do the "elses", in which case they should be placed in different blocks if possible.

  • A commonly-cited rule of thumb is to have between $128$ and $256$ threads per block.

Regla compuesta del rectángulo

En el uso de CUDA se recomienda que:

  • Users escriban código de CUDA C simple.

  • Utilicen las librerías ya hechas por NVIDIA o terceros para mantener simplicidad y eficiencia en el código.

Lo anterior para disminuir el tiempo y la cantidad de código que users tengan que hacer (o rehacer) y puesto que dominar la programación de CUDA C requiere una buena inversión de tiempo.

Así, tenemos a Thrust una template library basada en la Standard Template Library (STL) de C++ construída por NVIDIA que de acuerdo a su documentación:

"Thrust provides a rich collection of data parallel primitives such as scan, sort, and reduce, which can be composed together to implement complex algorithms with concise, readable source code. By describing your computation in terms of these high-level abstractions you provide Thrust with the freedom to select the most efficient implementation automatically. As a result, Thrust can be utilized in rapid prototyping of CUDA applications, where programmer productivity matters most, as well as in production, where robustness and absolute performance are crucial."

Thrust tiene la opción de utilizarse con OpenMP, Thread Building Blocks (TBB) y con CUDA-C++. Ver por ejemplo Device Backends para conocer cómo cambiar entre OpenMP y CUDA-C++, lo cual se realiza en la compilación y sin hacer cambios en el código!. A los sistemas de software que tienen este feature se les llama Heterogeneous computing.

Si se instala el CUDA toolkit o se utiliza la imagen de docker descrita al inicio de la nota, los headers en la librería template de Thrust estarán disponibles para su uso.

En el siguiente ejemplo de la regla del rectángulo compuesta se utiliza:

y se hace explícito el uso de la política de ejecucion thrust::device.

Referencias para el programa siguiente se encuentran en thrust inside user written kernels y cuda how to sum all elements of an array into one number within the gpu de stackoverflow.

Los siguientes tiempos se calcularon con una GPU de las siguientes características:


In [31]:
%%bash

./device_properties.out


----------------------
----device 0 ----
Device Name: GeForce GTX 750
Compute capability: 5.0
Clock rate: 1293500
Unified memory: 1
 ---Memory Information for device 0 (results on bytes)---
Total global mem: 1025769472
Total constant Mem: 65536
Shared memory per thread block: 49152
Shared memory per SM: 65536
 ---MP Information for device 0 ---
SM count: 4
Threads in warp: 32
Max threads per SM: 2048
Max warps per SM: 64
Max threads per block: 1024
Max thread dimensions: (1024, 1024, 64)
Max grid dimensions: (2147483647, 65535, 65535)

$n=10^3$ subintervalos


In [1]:
%%file Rcf.cu
#include<stdio.h>
#include <thrust/reduce.h>
#include <thrust/execution_policy.h>

__global__ void Rcf(double *data, double a, double h_hat, int n, double *sum ) {
    /*
    Compute numerical approximation using rectangle or mid-point method in 
    an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
        a (int): left point of interval
        n (int): number of subintervals
        data (double): array that will hold values evaluated in function
    Returns:
        sum (double): pointer to result
    */
    double x=0.0;
    if(threadIdx.x<=n-1){
        x=a+(threadIdx.x+1/2.0)*h_hat;
        data[threadIdx.x]=std::exp(-std::pow(x,2));
    }
    if(threadIdx.x==0){
        *sum = thrust::reduce(thrust::device, data , data + n, (double)0, thrust::plus<double>());
    }
}

int main(int argc, char *argv[]){
    double suma=0.0;
    double *d_data;
    double *d_suma;
    double a=0.0, b=1.0;
    double h_hat;
    int n=1e3; //número de subintervalos
    double objetivo=0.7468241328124271;
    double time_spent;
    clock_t begin,end;
    cudaMalloc((void **)&d_data,sizeof(double)*n);
    cudaMalloc((void**)&d_suma,sizeof(double));
    h_hat=(b-a)/n;
    begin=clock();
    Rcf<<<1,n>>>(d_data, a,h_hat,n,d_suma); //1 bloque de n threads
    cudaDeviceSynchronize();
    end=clock();
    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    cudaMemcpy(&suma, d_suma, sizeof(double), cudaMemcpyDeviceToHost);
    suma=h_hat*suma;
    cudaFree(d_data) ;
    cudaFree(d_suma) ;
    printf("Integral de %f a %f = %1.15e\n", a,b,suma);
    printf("Error relativo de la solución: %1.15e\n", fabs(suma-objetivo)/fabs(objetivo));
    printf("Tiempo de cálculo en la gpu %.5f\n", time_spent);
    return 0;
}


Writing Rcf.cu

In [33]:
%%bash
nvcc --compiler-options -Wall Rcf.cu -o Rcf.out

In [34]:
%%bash
./Rcf.out


Integral de 0.000000 a 1.000000 = 7.468241634690490e-01
Error relativo de la solución: 4.104931878976858e-08
Tiempo de cálculo en la gpu 0.00009

In [35]:
%%bash
nvprof --normalized-time-unit s ./Rcf.out


Integral de 0.000000 a 1.000000 = 7.468241634690490e-01
Error relativo de la solución: 4.104931878976858e-08
Tiempo de cálculo en la gpu 0.00010
==525== NVPROF is profiling process 525, command: ./Rcf.out
==525== Profiling application: ./Rcf.out
==525== Profiling result:
            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
                        %         s                   s         s         s
 GPU activities:    98.00  6.89e-05         1  6.89e-05  6.89e-05  6.89e-05  Rcf(double*, double, double, int, double*)
                     2.00  1.41e-06         1  1.41e-06  1.41e-06  1.41e-06  [CUDA memcpy DtoH]
      API calls:    99.49  0.088044         2  0.044022  6.13e-06  0.088038  cudaMalloc
                     0.21  1.83e-04        97  1.89e-06  3.26e-07  7.10e-05  cuDeviceGetAttribute
                     0.09  7.86e-05         1  7.86e-05  7.86e-05  7.86e-05  cuDeviceTotalMem
                     0.08  7.15e-05         1  7.15e-05  7.15e-05  7.15e-05  cudaDeviceSynchronize
                     0.06  5.31e-05         2  2.65e-05  6.30e-06  4.68e-05  cudaFree
                     0.03  2.80e-05         1  2.80e-05  2.80e-05  2.80e-05  cuDeviceGetName
                     0.02  1.97e-05         1  1.97e-05  1.97e-05  1.97e-05  cudaLaunchKernel
                     0.02  1.44e-05         1  1.44e-05  1.44e-05  1.44e-05  cudaMemcpy
                     0.00  2.11e-06         3  7.03e-07  3.96e-07  1.29e-06  cuDeviceGetCount
                     0.00  2.10e-06         1  2.10e-06  2.10e-06  2.10e-06  cuDeviceGetPCIBusId
                     0.00  1.23e-06         2  6.17e-07  3.24e-07  9.10e-07  cuDeviceGet
                     0.00  5.00e-07         1  5.00e-07  5.00e-07  5.00e-07  cuDeviceGetUuid

$n=1025$ subintervalos


In [36]:
%%file Rcf2.cu
#include<stdio.h>
#include <thrust/reduce.h>
#include <thrust/execution_policy.h>

__global__ void Rcf(double *data, double a, double h_hat, int n, double *sum ) {
    /*
    Compute numerical approximation using rectangle or mid-point method in 
    an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
        a (int): left point of interval
        n (int): number of subintervals
        data (double): array that will hold values evaluated in function
    Returns:
        sum (double): pointer to result
    */
    double x=0.0;
    if(threadIdx.x<=n-1){
        x=a+(threadIdx.x+1/2.0)*h_hat;
        data[threadIdx.x]=std::exp(-std::pow(x,2));
    }
    if(threadIdx.x==0){
        *sum = thrust::reduce(thrust::device, data , data + n, (double)0, thrust::plus<double>());
    }
}

int main(int argc, char *argv[]){
    double suma=0.0;
    double *d_data;
    double *d_suma;
    double a=0.0, b=1.0;
    double h_hat;
    int n=1025; //número de subintervalos
    double objetivo=0.7468241328124271;
    double time_spent;
    clock_t begin,end;
    cudaMalloc((void **)&d_data,sizeof(double)*n);
    cudaMalloc((void**)&d_suma,sizeof(double));
    h_hat=(b-a)/n;
    begin=clock();
    Rcf<<<1,n>>>(d_data, a,h_hat,n,d_suma); //1 bloque de n threads
    cudaDeviceSynchronize();
    end=clock();
    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    cudaMemcpy(&suma, d_suma, sizeof(double), cudaMemcpyDeviceToHost);
    suma=h_hat*suma;
    cudaFree(d_data) ;
    cudaFree(d_suma) ;
    printf("Integral de %f a %f = %1.15e\n", a,b,suma);
    printf("Error relativo de la solución: %1.15e\n", fabs(suma-objetivo)/fabs(objetivo));
    printf("Tiempo de cálculo en la gpu %.5f\n", time_spent);
    return 0;
}


Writing Rcf2.cu

In [37]:
%%bash
nvcc --compiler-options -Wall Rcf2.cu -o Rcf2.out

In [38]:
%%bash
./Rcf2.out


Integral de 0.000000 a 1.000000 = 0.000000000000000e+00
Error relativo de la solución: 1.000000000000000e+00
Tiempo de cálculo en la gpu 0.00001

Obsérvese error relativo de $100\%$

¿Cómo lo arreglamos?


In [39]:
%%file Rcf3.cu
#include<stdio.h>
#include <thrust/reduce.h>
#include <thrust/execution_policy.h>

__global__ void Rcf(double *data, double a, double h_hat, int n, double *sum) {
    /*
    Compute numerical approximation using rectangle or mid-point method in 
    an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
        a (int): left point of interval
        n (int): number of subintervals
        data (double): array that will hold values evaluated in function
    Returns:
        sum (double): pointer to result
    */
    double x=0.0;
    int stride=0;
    if(threadIdx.x<=n-1){
        x=a+(threadIdx.x+1/2.0)*h_hat;
        data[threadIdx.x]=std::exp(-std::pow(x,2));
    }
    if(threadIdx.x==0){
        stride=blockDim.x;
        x=a+(threadIdx.x+stride+1/2.0)*h_hat;
        data[threadIdx.x+stride]=std::exp(-std::pow(x,2));
        *sum = thrust::reduce(thrust::device, data , data + n, (double)0, thrust::plus<double>());
    }
}

int main(int argc, char *argv[]){
    double suma=0.0;
    double *d_data;
    double *d_suma;
    double a=0.0, b=1.0;
    double h_hat;
    int n_threads_per_block=1024; 
    int n_bloques=2;
    int n=1025;//número de subintervalos
    double objetivo=0.7468241328124271;
    double time_spent;
    clock_t begin,end;
    cudaMalloc((void **)&d_data,sizeof(double)*n);
    cudaMalloc((void**)&d_suma,sizeof(double));
    h_hat=(b-a)/n;
    begin=clock();
    Rcf<<<n_bloques,n_threads_per_block>>>(d_data, a,h_hat,n,d_suma); 
    cudaDeviceSynchronize();
    end=clock();
    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    cudaMemcpy(&suma, d_suma, sizeof(double), cudaMemcpyDeviceToHost);
    suma=h_hat*suma;
    cudaFree(d_data) ;
    cudaFree(d_suma) ;
    printf("Integral de %f a %f = %1.15e\n", a,b,suma);
    printf("Error relativo de la solución: %1.15e\n", fabs(suma-objetivo)/fabs(objetivo));
    printf("Tiempo de cálculo en la gpu %.5f\n", time_spent);
    return 0;
}


Writing Rcf3.cu

In [40]:
%%bash
nvcc --compiler-options -Wall Rcf3.cu -o Rcf3.out

In [41]:
%%bash
./Rcf3.out


Integral de 0.000000 a 1.000000 = 7.468241619918411e-01
Error relativo de la solución: 3.907133247860604e-08
Tiempo de cálculo en la gpu 0.00009

Pero en la propuesta anterior lanzamos $2*1024$ (bloques $\times$ número de threads) $=2048$ threads y sólo ocupamos $1025$ threads. Entonces podemos cambiar el código anterior para aprovechar los $2048$ threads como sigue:


In [42]:
%%file Rcf4.cu
#include<stdio.h>
#include <thrust/reduce.h>
#include <thrust/execution_policy.h>

__global__ void Rcf(double *data, double a, double h_hat, int n, double *sum) {
    /*
    Compute numerical approximation using rectangle or mid-point method in 
    an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
        a (int): left point of interval
        n (int): number of subintervals
        data (double): array that will hold values evaluated in function
    Returns:
        sum (double): pointer to result
    */
    double x=0.0;
    int stride=0;
    int i;
    stride=blockDim.x;
    for(i=threadIdx.x;i<=n-1;i+=stride){
        if(i<=n-1){
            x=a+(i+1/2.0)*h_hat;
            data[i]=std::exp(-std::pow(x,2));
        }
    }
    if(threadIdx.x==0){
        *sum = thrust::reduce(thrust::device, data , data + n, (double)0, thrust::plus<double>());
    }
}

int main(int argc, char *argv[]){
    double suma=0.0;
    double *d_data;
    double *d_suma;
    double a=0.0, b=1.0;
    double h_hat;
    int n_threads_per_block=1024; 
    int n_bloques=2;
    int n=n_threads_per_block*n_bloques;//número de subintervalos
    double objetivo=0.7468241328124271;
    double time_spent;
    clock_t begin,end;
    cudaMalloc((void **)&d_data,sizeof(double)*n);
    cudaMalloc((void**)&d_suma,sizeof(double));
    h_hat=(b-a)/n;
    begin=clock();
    Rcf<<<n_bloques,n_threads_per_block>>>(d_data, a,h_hat,n,d_suma); 
    cudaDeviceSynchronize();
    end=clock();
    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    cudaMemcpy(&suma, d_suma, sizeof(double), cudaMemcpyDeviceToHost);
    suma=h_hat*suma;
    cudaFree(d_data) ;
    cudaFree(d_suma) ;
    printf("Integral de %f a %f = %1.15e\n", a,b,suma);
    printf("Error relativo de la solución: %1.15e\n", fabs(suma-objetivo)/fabs(objetivo));
    printf("Tiempo de cálculo en la gpu %.5f\n", time_spent);
    return 0;
}


Writing Rcf4.cu

In [43]:
%%bash
nvcc --compiler-options -Wall Rcf4.cu -o Rcf4.out

In [44]:
%%bash
./Rcf4.out


Integral de 0.000000 a 1.000000 = 7.468241401215338e-01
Error relativo de la solución: 9.786918140590463e-09
Tiempo de cálculo en la gpu 0.00016

Pero todavía podemos hacer aún mejor sin el for:


In [45]:
%%file Rcf5.cu
#include<stdio.h>
#include <thrust/reduce.h>
#include <thrust/execution_policy.h>

__global__ void Rcf(double *data, double a, double h_hat, int n, double *sum ) {
    /*
    Compute numerical approximation using rectangle or mid-point method in 
    an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
        a (int): left point of interval
        n (int): number of subintervals
        data (double): array that will hold values evaluated in function
    Returns:
        sum (double): pointer to result
    */
    double x=0.0;
    int idx;
    idx = blockIdx.x * blockDim.x + threadIdx.x;
    if(idx<=n-1){
        x=a+(idx+1/2.0)*h_hat;
        data[idx]=std::exp(-std::pow(x,2));
    }
    if(idx==0){
        *sum = thrust::reduce(thrust::device, data , data + n, (double)0, thrust::plus<double>());
    }
}

int main(int argc, char *argv[]){
    double suma=0.0;
    double *d_data;
    double *d_suma;
    double a=0.0, b=1.0;
    double h_hat;
    int n_threads_per_block=1024; 
    int n_bloques=2;
    double objetivo=0.7468241328124271;
    int n=n_bloques*n_threads_per_block;//número de subintervalos
    double time_spent;
    clock_t begin,end;
    cudaMalloc((void **)&d_data,sizeof(double)*n);
    cudaMalloc((void**)&d_suma,sizeof(double));
    h_hat=(b-a)/n;
    begin = clock();
    Rcf<<<n_bloques,n_threads_per_block>>>(d_data, a,h_hat,n,d_suma); 
    cudaDeviceSynchronize();
    end = clock();
    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    cudaMemcpy(&suma, d_suma, sizeof(double), cudaMemcpyDeviceToHost);
    suma=h_hat*suma;
    cudaFree(d_data) ;
    cudaFree(d_suma) ;
    printf("Integral de %f a %f = %1.15e\n", a,b,suma);
    printf("Error relativo de la solución: %1.15e\n", fabs(suma-objetivo)/fabs(objetivo));
    printf("Tiempo de cálculo en la gpu %.5f\n", time_spent);
    return 0;
}


Writing Rcf5.cu

In [46]:
%%bash
nvcc --compiler-options -Wall Rcf5.cu -o Rcf5.out

In [47]:
%%bash
./Rcf5.out


Integral de 0.000000 a 1.000000 = 7.468241401215338e-01
Error relativo de la solución: 9.786918140590463e-09
Tiempo de cálculo en la gpu 0.00014

In [48]:
%%bash
nvprof --normalized-time-unit s ./Rcf5.out


Integral de 0.000000 a 1.000000 = 7.468241401215338e-01
Error relativo de la solución: 9.786918140590463e-09
Tiempo de cálculo en la gpu 0.00014
==694== NVPROF is profiling process 694, command: ./Rcf5.out
==694== Profiling application: ./Rcf5.out
==694== Profiling result:
            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
                        %         s                   s         s         s
 GPU activities:    98.76  1.17e-04         1  1.17e-04  1.17e-04  1.17e-04  Rcf(double*, double, double, int, double*)
                     1.24  1.47e-06         1  1.47e-06  1.47e-06  1.47e-06  [CUDA memcpy DtoH]
      API calls:    99.36  0.093178         2  0.046589  4.18e-06  0.093174  cudaMalloc
                     0.31  2.88e-04        97  2.96e-06  2.85e-07  1.57e-04  cuDeviceGetAttribute
                     0.13  1.20e-04         1  1.20e-04  1.20e-04  1.20e-04  cudaDeviceSynchronize
                     0.08  7.32e-05         1  7.32e-05  7.32e-05  7.32e-05  cuDeviceTotalMem
                     0.06  5.18e-05         2  2.59e-05  5.26e-06  4.66e-05  cudaFree
                     0.03  2.69e-05         1  2.69e-05  2.69e-05  2.69e-05  cuDeviceGetName
                     0.02  1.98e-05         1  1.98e-05  1.98e-05  1.98e-05  cudaLaunchKernel
                     0.02  1.45e-05         1  1.45e-05  1.45e-05  1.45e-05  cudaMemcpy
                     0.00  2.12e-06         3  7.07e-07  3.85e-07  1.34e-06  cuDeviceGetCount
                     0.00  1.82e-06         1  1.82e-06  1.82e-06  1.82e-06  cuDeviceGetPCIBusId
                     0.00  9.84e-07         2  4.92e-07  3.02e-07  6.82e-07  cuDeviceGet
                     0.00  4.23e-07         1  4.23e-07  4.23e-07  4.23e-07  cuDeviceGetUuid

Para una visualización sobre la construcción del índice en el kernel utilizando blockDim.x*blockIdx.x + threadIdx.x ver An Even Easier Introduction to CUDA.

¿Más nodos?

Para este caso, incrementamos el número de bloques:


In [67]:
%%file Rcf6.cu
#include<stdio.h>
#include <thrust/reduce.h>
#include <thrust/execution_policy.h>

__global__ void Rcf(double *data, double a, double h_hat, int n, double *sum ) {
    /*
    Compute numerical approximation using rectangle or mid-point method in 
    an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
        a (int): left point of interval
        n (int): number of subintervals
        data (double): array that will hold values evaluated in function
    Returns:
        sum (double): pointer to result
    */
    double x=0.0;
    int idx;
    idx = blockIdx.x * blockDim.x + threadIdx.x;
    if(idx<=n-1){
        x=a+(idx+1/2.0)*h_hat;
        data[idx]=std::exp(-std::pow(x,2));
    }
    if(idx==0){
        *sum = thrust::reduce(thrust::device, data , data + n, (double)0, thrust::plus<double>());
    }
}

int main(int argc, char *argv[]){
    double suma=0.0;
    double *d_data;
    double *d_suma;
    double a=0.0, b=1.0;
    double h_hat;
    int n_threads_per_block=1024; 
    int n_bloques=0; //inicializamos en 0
    double objetivo=0.7468241328124271;
    int n=0; //inicializamos en 0
    double time_spent;
    clock_t begin,end;
    cudaDeviceProp properties;
    cudaGetDeviceProperties(&properties, 0);
    n_bloques = 256 * properties.multiProcessorCount; //elegimos un múltiplo del número de SM's
                                                    //y además que sea menor a 10^4 (para tener 10^6 aprox subintervalos)
    n = n_bloques*n_threads_per_block;//número de subintervalos
    cudaMalloc((void **)&d_data,sizeof(double)*n);
    cudaMalloc((void**)&d_suma,sizeof(double));
    h_hat=(b-a)/n;
    begin = clock();
    Rcf<<<n_bloques,n_threads_per_block>>>(d_data, a,h_hat,n,d_suma); 
    cudaDeviceSynchronize();
    end = clock();
    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    cudaMemcpy(&suma, d_suma, sizeof(double), cudaMemcpyDeviceToHost);
    suma=h_hat*suma;
    cudaFree(d_data) ;
    cudaFree(d_suma) ;
    printf("Integral de %f a %f = %1.15e\n", a,b,suma);
    printf("Error relativo de la solución: %1.15e\n", fabs(suma-objetivo)/fabs(objetivo));
    printf("Tiempo de cálculo en la gpu %.5f\n", time_spent);
    return 0;
}


Overwriting Rcf6.cu

In [68]:
%%bash
nvcc --compiler-options -Wall Rcf6.cu -o Rcf6.out

In [69]:
%%bash
./Rcf6.out


Integral de 0.000000 a 1.000000 = 7.468241328124761e-01
Error relativo de la solución: 6.555872157155699e-14
Tiempo de cálculo en la gpu 0.06257

In [70]:
%%bash
nvprof --normalized-time-unit s ./Rcf6.out


Integral de 0.000000 a 1.000000 = 7.468241328124761e-01
Error relativo de la solución: 6.555872157155699e-14
Tiempo de cálculo en la gpu 0.06288
==963== NVPROF is profiling process 963, command: ./Rcf6.out
==963== Profiling application: ./Rcf6.out
==963== Profiling result:
            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
                        %         s                   s         s         s
 GPU activities:   100.00  0.062851         1  0.062851  0.062851  0.062851  Rcf(double*, double, double, int, double*)
                     0.00  1.57e-06         1  1.57e-06  1.57e-06  1.57e-06  [CUDA memcpy DtoH]
      API calls:    58.26  0.089315         2  0.044657  5.76e-05  0.089257  cudaMalloc
                    41.00  0.062855         1  0.062855  0.062855  0.062855  cudaDeviceSynchronize
                     0.41  6.22e-04         2  3.11e-04  5.54e-05  5.67e-04  cudaFree
                     0.13  1.94e-04        97  2.00e-06  3.36e-07  7.47e-05  cuDeviceGetAttribute
                     0.10  1.58e-04         1  1.58e-04  1.58e-04  1.58e-04  cudaGetDeviceProperties
                     0.05  8.27e-05         1  8.27e-05  8.27e-05  8.27e-05  cuDeviceTotalMem
                     0.02  2.92e-05         1  2.92e-05  2.92e-05  2.92e-05  cuDeviceGetName
                     0.01  2.01e-05         1  2.01e-05  2.01e-05  2.01e-05  cudaLaunchKernel
                     0.01  1.77e-05         1  1.77e-05  1.77e-05  1.77e-05  cudaMemcpy
                     0.00  2.19e-06         1  2.19e-06  2.19e-06  2.19e-06  cuDeviceGetPCIBusId
                     0.00  1.99e-06         3  6.63e-07  3.77e-07  1.16e-06  cuDeviceGetCount
                     0.00  1.33e-06         2  6.66e-07  3.45e-07  9.87e-07  cuDeviceGet
                     0.00  5.68e-07         1  5.68e-07  5.68e-07  5.68e-07  cuDeviceGetUuid

$n=16777216 \approx 10^7$ subintervalos?


In [90]:
%%file Rcf7.cu
#include<stdio.h>
#include <thrust/reduce.h>
#include <thrust/execution_policy.h>

__global__ void Rcf(double *data, double a, double h_hat, int n, double *sum ) {
    /*
    Compute numerical approximation using rectangle or mid-point method in 
    an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
        a (int): left point of interval
        n (int): number of subintervals
        data (double): array that will hold values evaluated in function
    Returns:
        sum (double): pointer to result
    */
    double x=0.0;
    int idx;
    idx = blockIdx.x * blockDim.x + threadIdx.x;
    if(idx<=n-1){
        x=a+(idx+1/2.0)*h_hat;
        data[idx]=std::exp(-std::pow(x,2));
    }
    if(idx==0){
        *sum = thrust::reduce(thrust::device, data , data + n, (double)0, thrust::plus<double>());
    }
}

int main(int argc, char *argv[]){
    double suma=0.0;
    double *d_data;
    double *d_suma;
    double a=0.0, b=1.0;
    double h_hat;
    int n_threads_per_block=1024; 
    int n_bloques=0; //inicializamos en 0
    double objetivo=0.7468241328124271;
    int n=0; //inicializamos en 0
    double time_spent;
    clock_t begin,end;
    cudaDeviceProp properties;
    cudaGetDeviceProperties(&properties, 0);
    n_bloques = 4096 * properties.multiProcessorCount; //elegimos un múltiplo del número de SM's
                                                    //y además que sea menor a 10^8 (para tener 10^7 aprox subintervalos)
    n = n_bloques*n_threads_per_block;//número de subintervalos
    cudaMalloc((void **)&d_data,sizeof(double)*n);
    cudaMalloc((void**)&d_suma,sizeof(double));
    h_hat=(b-a)/n;
    begin = clock();
    Rcf<<<n_bloques,n_threads_per_block>>>(d_data, a,h_hat,n,d_suma); 
    cudaDeviceSynchronize();
    end = clock();
    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    cudaMemcpy(&suma, d_suma, sizeof(double), cudaMemcpyDeviceToHost);
    suma=h_hat*suma;
    cudaFree(d_data) ;
    cudaFree(d_suma) ;
    printf("Integral de %f a %f = %1.15e\n", a,b,suma);
    printf("Error relativo de la solución: %1.15e\n", fabs(suma-objetivo)/fabs(objetivo));
    printf("Tiempo de cálculo en la gpu %.5f\n", time_spent);
    return 0;
}


Overwriting Rcf7.cu

In [80]:
%%bash
nvcc --compiler-options -Wall Rcf7.cu -o Rcf7.out

In [81]:
%%bash
./Rcf7.out


Integral de 0.000000 a 1.000000 = 7.468241328123114e-01
Error relativo de la solución: 1.549029203572843e-13
Tiempo de cálculo en la gpu 0.94370

In [75]:
%%bash
nvprof --normalized-time-unit s ./Rcf7.out


Integral de 0.000000 a 1.000000 = 7.468241328123114e-01
Error relativo de la solución: 1.549029203572843e-13
Tiempo de cálculo en la gpu 0.94049
==1118== NVPROF is profiling process 1118, command: ./Rcf7.out
==1118== Profiling application: ./Rcf7.out
==1118== Profiling result:
            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
                        %         s                   s         s         s
 GPU activities:   100.00  0.939793         1  0.939793  0.939793  0.939793  Rcf(double*, double, double, int, double*)
                     0.00  1.44e-06         1  1.44e-06  1.44e-06  1.44e-06  [CUDA memcpy DtoH]
      API calls:    90.36  0.939797         1  0.939797  0.939797  0.939797  cudaDeviceSynchronize
                     8.80  0.091559         2  0.045780  6.07e-05  0.091499  cudaMalloc
                     0.76  7.94e-03         2  3.97e-03  1.58e-04  7.78e-03  cudaFree
                     0.03  3.02e-04        97  3.11e-06  3.18e-07  1.61e-04  cuDeviceGetAttribute
                     0.03  2.61e-04         1  2.61e-04  2.61e-04  2.61e-04  cudaGetDeviceProperties
                     0.01  8.03e-05         1  8.03e-05  8.03e-05  8.03e-05  cuDeviceTotalMem
                     0.00  2.99e-05         1  2.99e-05  2.99e-05  2.99e-05  cuDeviceGetName
                     0.00  2.02e-05         1  2.02e-05  2.02e-05  2.02e-05  cudaLaunchKernel
                     0.00  1.82e-05         1  1.82e-05  1.82e-05  1.82e-05  cudaMemcpy
                     0.00  2.19e-06         3  7.29e-07  4.02e-07  1.35e-06  cuDeviceGetCount
                     0.00  2.10e-06         1  2.10e-06  2.10e-06  2.10e-06  cuDeviceGetPCIBusId
                     0.00  1.33e-06         2  6.65e-07  3.42e-07  9.89e-07  cuDeviceGet
                     0.00  5.26e-07         1  5.26e-07  5.26e-07  5.26e-07  cuDeviceGetUuid

$n=117440512 \approx 10^8$ subintervalos?

Este caso se ejecutó en una GPU con las siguientes características:

----------------------
----device 0 ----
Device Name: Tesla K20Xm
Compute capability: 3.5
Clock rate: 732000
Unified memory: 1
 ---Memory Information for device 0 (results on bytes)---
Total global mem: 5977800704
Total constant Mem: 65536
Shared memory per thread block: 49152
Shared memory per SM: 49152
 ---MP Information for device 0 ---
SM count: 14
Threads in warp: 32
Max threads per SM: 2048
Max warps per SM: 64
Max threads per block: 1024
Max thread dimensions: (1024, 1024, 64)
Max grid dimensions: (2147483647, 65535, 65535)
----------------------
----device 1 ----
Device Name: Tesla K20Xm
Compute capability: 3.5
Clock rate: 732000
Unified memory: 1
 ---Memory Information for device 1 (results on bytes)---
Total global mem: 5977800704
Total constant Mem: 65536
Shared memory per thread block: 49152
Shared memory per SM: 49152
 ---MP Information for device 1 ---
SM count: 14
Threads in warp: 32
Max threads per SM: 2048
Max warps per SM: 64
Max threads per block: 1024
Max thread dimensions: (1024, 1024, 64)
Max grid dimensions: (2147483647, 65535, 65535)

In [1]:
%%file Rcf8.cu
#include<stdio.h>
#include <thrust/reduce.h>
#include <thrust/execution_policy.h>

__global__ void Rcf(double *data, double a, double h_hat, long int n, double *sum ) {
    /*
    Compute numerical approximation using rectangle or mid-point method in 
    an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
        a (int): left point of interval
        n (int): number of subintervals
        data (double): array that will hold values evaluated in function
    Returns:
        sum (double): pointer to result
    */
    double x=0.0;
    int idx;
    idx = blockIdx.x * blockDim.x + threadIdx.x;
    if(idx<=n-1){
        x=a+(idx+1/2.0)*h_hat;
        data[idx]=std::exp(-std::pow(x,2));
    }
    //esta alternativa no funciona... si el número de bloques rebasa el límite:
    //int num_threads=gridDim.x * blockDim.x;
    //int stride = num_threads;
    //for(i=idx; i<=n-1; i+=stride){
    //    if(idx<=n-1){
    //        x=a+(idx+1/2.0)*h_hat;
    //        data[idx]=std::exp(-std::pow(x,2));
    //    }
    //}
    
    if(idx==0){
        *sum = thrust::reduce(thrust::device, data , data + n, (double)0, thrust::plus<double>());
    }
}

int main(int argc, char *argv[]){
    double suma=0.0;
    double *d_data;
    double *d_suma;
    double a=0.0, b=1.0;
    double h_hat;
    int n_threads_per_block=1024; 
    long int n_bloques=0; //inicializamos en 0
    double objetivo=0.7468241328124271;
    long int n=0; //inicializamos en 0
    double time_spent;
    clock_t begin,end;
    cudaDeviceProp properties;    
    cudaGetDeviceProperties(&properties, 0);
    n_bloques = 8192 * properties.multiProcessorCount; //elegimos un múltiplo del número de SM's
                                                    //y además que sea menor a 10^9 (para tener 10^8 aprox subintervalos)
    n = n_bloques*n_threads_per_block;//número de subintervalos
    dim3 dimGrid(n_bloques,1,1);
    dim3 dimBlock(n_threads_per_block,1,1);
    cudaMalloc((void **)&d_data,sizeof(double)*n);
    cudaMalloc((void**)&d_suma,sizeof(double));
    h_hat=(b-a)/n;
    begin = clock();
    Rcf<<<dimGrid,dimBlock>>>(d_data, a,h_hat,n,d_suma); 
    cudaDeviceSynchronize();
    end = clock();
    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    cudaMemcpy(&suma, d_suma, sizeof(double), cudaMemcpyDeviceToHost);
    suma=h_hat*suma;
    cudaFree(d_data) ;
    cudaFree(d_suma) ;
    printf("Integral de %f a %f = %1.15e\n", a,b,suma);
    printf("Error relativo de la solución: %1.15e\n", fabs(suma-objetivo)/fabs(objetivo));
    printf("Tiempo de cálculo en la gpu %.5f\n", time_spent);
    return 0;
}


Writing Rcf8.cu

In [5]:
%%bash
nvcc --compiler-options -Wall Rcf8.cu -o Rcf8.out

In [6]:
%%bash
./Rcf8.out


Integral de 0.000000 a 1.000000 = 7.468241328125712e-01
Error relativo de la solución: 1.929596838999568e-13
Tiempo de cálculo en la gpu 12.95516

In [7]:
%%bash
nvprof --normalized-time-unit s ./Rcf8.out


Integral de 0.000000 a 1.000000 = 7.468241328125712e-01
Error relativo de la solución: 1.929596838999568e-13
Tiempo de cálculo en la gpu 13.00077
==173== NVPROF is profiling process 173, command: ./Rcf8.out
==173== Warning: Profiling results might be incorrect with current version of nvcc compiler used to compile cuda app. Compile with nvcc compiler 9.0 or later version to get correct profiling results. Ignore this warning if code is already compiled with the recommended nvcc version 
==173== Profiling application: ./Rcf8.out
==173== Profiling result:
            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
                        %         s                   s         s         s
 GPU activities:   100.00  12.99435         1  12.99435  12.99435  12.99435  Rcf(double*, double, double, long, double*)
                     0.00  2.46e-06         1  2.46e-06  2.46e-06  2.46e-06  [CUDA memcpy DtoH]
      API calls:    97.08  12.99439         1  12.99439  12.99439  12.99439  cudaDeviceSynchronize
                     2.21  0.296135         2  0.148068  1.54e-04  0.295981  cudaMalloc
                     0.69  0.092310         2  0.046155  1.33e-03  0.090979  cudaFree
                     0.00  5.54e-04       194  2.85e-06  2.41e-07  1.05e-04  cuDeviceGetAttribute
                     0.00  5.52e-04         2  2.76e-04  2.57e-04  2.95e-04  cuDeviceTotalMem
                     0.00  3.02e-04         1  3.02e-04  3.02e-04  3.02e-04  cudaLaunchKernel
                     0.00  2.43e-04         1  2.43e-04  2.43e-04  2.43e-04  cudaGetDeviceProperties
                     0.00  1.02e-04         2  5.10e-05  2.03e-05  8.17e-05  cuDeviceGetName
                     0.00  8.62e-05         1  8.62e-05  8.62e-05  8.62e-05  cudaMemcpy
                     0.00  6.53e-06         2  3.27e-06  2.73e-06  3.80e-06  cuDeviceGetPCIBusId
                     0.00  2.42e-06         4  6.06e-07  2.69e-07  1.19e-06  cuDeviceGet
                     0.00  1.94e-06         3  6.48e-07  3.04e-07  1.05e-06  cuDeviceGetCount
                     0.00  8.75e-07         2  4.37e-07  3.71e-07  5.04e-07  cuDeviceGetUuid

Obs: en la programación con CUDA-C es importante checar posibles errores de alojamiento de memoria. Una forma es con los tipos cudaError_t y cudaSuccess . Ver liga:


In [1]:
%%file Rcf9.cu
#include<stdio.h>
#include <thrust/reduce.h>
#include <thrust/execution_policy.h>

__global__ void Rcf(double *data, double a, double h_hat, long int n, double *sum ) {
    /*
    Compute numerical approximation using rectangle or mid-point method in 
    an interval.
    Nodes are generated via formula: x_i = a+(i+1/2)h_hat for i=0,1,...,n-1 and h_hat=(b-a)/n
    Args:
        a (int): left point of interval
        n (int): number of subintervals
        data (double): array that will hold values evaluated in function
    Returns:
        sum (double): pointer to result
    */
    double x=0.0;
    int idx;
    int num_threads=gridDim.x * blockDim.x;
    int stride = num_threads;
    int i;
    idx = blockIdx.x * blockDim.x + threadIdx.x;
    for(i=idx; i<=n-1; i+=stride){
        if(idx<=n-1){
            x=a+(idx+1/2.0)*h_hat;
            data[idx]=std::exp(-std::pow(x,2));
        }
}
    
    if(idx==0){
        *sum = thrust::reduce(thrust::device, data , data + n, (double)0, thrust::plus<double>());
    }
}

cudaError_t fun_checa_error(cudaError_t result) {
    if (result != cudaSuccess) {
        fprintf(stderr, "Error: %s\n", cudaGetErrorString(result));
    }
    return result;
}

int main(int argc, char *argv[]){
    double suma=0.0;
    double *d_data;
    double *d_suma;
    double a=0.0, b=1.0;
    double h_hat;
    int n_threads_per_block=1024; 
    long int n_bloques=0; //inicializamos en 0
    double objetivo=0.7468241328124271;
    long int n=0; //inicializamos en 0
    double time_spent;
    clock_t begin,end;
    cudaDeviceProp properties;    
    cudaGetDeviceProperties(&properties, 0);
    n_bloques = 32768 * properties.multiProcessorCount; //elegimos un múltiplo del número de SM's
                                                    //y además que sea menor a 10^9 (para tener 10^8 aprox subintervalos)
    n = n_bloques*n_threads_per_block;//número de subintervalos
    dim3 dimGrid(n_bloques,1,1);
    dim3 dimBlock(n_threads_per_block,1,1);
    fun_checa_error(cudaMalloc((void **)&d_data,sizeof(double)*n));
    fun_checa_error(cudaMalloc((void**)&d_suma,sizeof(double)));
    h_hat=(b-a)/n;
    begin = clock();
    Rcf<<<dimGrid,dimBlock>>>(d_data, a,h_hat,n,d_suma); 
    cudaDeviceSynchronize();
    end = clock();
    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    fun_checa_error(cudaMemcpy(&suma, d_suma, sizeof(double), cudaMemcpyDeviceToHost));
    suma=h_hat*suma;
    cudaFree(d_data);
    cudaFree(d_suma);
    printf("Integral de %f a %f = %1.15e\n", a,b,suma);
    printf("Error relativo de la solución: %1.15e\n", fabs(suma-objetivo)/fabs(objetivo));
    printf("Tiempo de cálculo en la gpu %.5f\n", time_spent);
    return 0;
}


Writing Rcf9.cu

In [2]:
%%bash
nvcc --compiler-options -Wall Rcf9.cu -o Rcf9.out

In [3]:
%%bash
./Rcf9.out


Integral de 0.000000 a 1.000000 = 0.000000000000000e+00
Error relativo de la solución: 1.000000000000000e+00
Tiempo de cálculo en la gpu 0.00048
Error: out of memory
Error: an illegal memory access was encountered

Preguntas de comprehensión:

1)¿Qué es y en qué consiste CUDA C?

2)¿Qué es un kernel?

3)¿Qué pieza de CUDA se encarga de asignar los bloques de cuda-threads a las SM’s?

4)¿Qué características (recursos compartidos, dimensiones, forma de agendar la ejecución en threads) tienen los bloques que se asignan a una SM al lanzarse y ejecutarse un kernel?

5)¿Qué es un warp?

6)Menciona los tipos de memorias que existen en las GPU’s.

7) Supón que tienes una tarjeta GT200 cuyas características son:

* Máximo número de threads que soporta una SM en un mismo instante en el tiempo: 1024
* Máximo número de threads en un bloque: 512
* Máximo número de bloques por SM: 8
* Número de SM’s que tiene esta GPU: 30

Responde:

a) ¿Cuál es la máxima cantidad de threads que puede soportar esta GPU en un mismo instante en el tiempo?
b) ¿Cuál es la máxima cantidad de warps por SM que puede soportar esta GPU en un mismo instante en el tiempo?
c) ¿Cuáles configuraciones de bloques y threads siguientes aprovechan la máxima cantidad de warps en una SM de esta GPU para un mismo instante en el tiempo?

    1.Una configuración del tipo: bloques de 64 threads y 16 bloques.
    2.Una configuración del tipo: bloques de 1024 threads y 1 bloque.
    3.Una configuración del tipo: bloques de 256 threads y 4 bloques.
    4.Una configuración del tipo: bloques de 512 threads y 8 bloques.

*Debes considerar las restricciones/características de la GPU dadas para responder pues algunas configuraciones infringen las mismas. No estamos considerando registers o shared memory.

Referencias

  1. N. Matloff, Parallel Computing for Data Science. With Examples in R, C++ and CUDA, 2014.

  2. CUDA

  3. 2.3.CUDA

Para más sobre Unified Memory revisar:

Es importante el manejo de errores por ejemplo en el alojamiento de memoria en la GPU. En este caso es útil revisar:

En stackoverflow encontramos a personas desarrolladoras de CUDA que resuelven preguntas muy útiles para continuar con el aprendizaje de CUDA C. Por ejemplo:

Otros sistemas de software para el Heterogeneous computing son: