Una GPU está compuesta por múltiples Streaming Multiprocessors (SM). Cada SM agrupa un conjunto de Streaming Processors (SP) junto con recursos como registros, memoria compartida (Shared Memory), caché L1, planificadores y unidades de carga/almacenamiento. Desde el punto de vista de la ejecución, los SP dentro de un SM se organizan en grupos llamados Warp, compartiendo un mismo controlador de memoria y de instrucciones.

La arquitectura se denomina SIMT (Single Instruction, Multiple Threads): en cada ciclo se emite una misma instrucción para todos los hilos de un warp. Un warp contiene 32 hilos. Por esta razón conviene que el número total de hilos lanzados sea múltiplo de 32; de lo contrario, el hardware genera una máscara de ejecución y los hilos sobrantes quedan inactivos, desperdiciando capacidad de cómputo.

Consulta de capacidades del dispositivo

La estructrua cudaDeviceProp permite conocer los límites físicos de la GPU. El siguiente programa muestra cómo recuperar esta información:

#include <iostream>
#include <cuda_runtime.h>

int main() {
    int gpu_id = 0;
    cudaDeviceProp prop;
    cudaGetDeviceProperties(&prop, gpu_id);

    std::cout << "GPU seleccionada: " << prop.name << std::endl;
    std::cout << "Multiprocesadores: " << prop.multiProcessorCount << std::endl;
    std::cout << "Bloques por SM: " << prop.maxBlocksPerMultiProcessor << std::endl;
    std::cout << "Hilos por SM: " << prop.maxThreadsPerMultiProcessor << std::endl;
    std::cout << "Hilos por bloque: " << prop.maxThreadsPerBlock << std::endl;
    std::cout << "Registros por bloque: " << prop.regsPerBlock << std::endl;
    std::cout << "Memoria compartida por bloque: "
              << prop.sharedMemPerBlock / 1024.0f << " KB" << std::endl;
    std::cout << "Warps por SM: "
              << prop.maxThreadsPerMultiProcessor / 32 << std::endl;
    std::cout << "Grid máximo: "
              << prop.maxGridSize[0] << " x "
              << prop.maxGridSize[1] << " x "
              << prop.maxGridSize[2] << std::endl;

    return 0;
}

Para compilar con CMake se puede emplear la sintaxis moderna para CUDA:

cmake_minimum_required(VERSION 3.18)
project(ejemplo CUDA)

find_package(CUDAToolkit REQUIRED)

add_executable(inspeccion inspeccion.cu)
target_link_libraries(inspeccion CUDA::cudart)

Modelo de programación: grids, bloques e hilos

En CUDA, el trabajo se organiza en una jerarquía de tres niveles:

• Grid: contiene todos los bloques que ejecutan un mismo kernel.
• Block: conjunto de hilos que pueden cooperar mediante memoria compartida y barreras.
• Thread: unidad mínima de ejecución.

Cada hilo dispone de registros y memoria local privada. Los hilos de un mismo bloque comparten memoria compartida. Todos los hilos del grid comparten la memoria global, constante y de texturas. Cada grid posee sus propias instancias de estas memorias globales.

Un SM ejecuta los warps de un bloque hasta completarlos antes de tomar otro bloque. El planificador puede alternar entre warps listos para ocultar la latencia de memoria y maximizar el uso de la unidad funcional.

Sintaxis de lenzamiento y configuración

Un kernel se invoca con la notación <<<Dg, Db, Ns, S>>>:

• Dg (dim3): dimensiones del grid. Hasta CUDA 10.x los límites en x e y eran 65535; arquitetcuras recientes permiten valores mucho mayores.
• Db (dim3): dimensiones del bloque. El producto total de hilos está limitado por la capacidad de cómputo (por ejemplo, 1024 en arquitecturas recientes).
• Ns (size_t): bytes adicionales de memoria compartida asignados dinámicamente por bloque.
• S (cudaStream_t): flujo de ejecución; por defecto el flujo 0.

Estrategias de indexado de hilos

Combinando dimensiones 1D, 2D y 3D para grids y bloques se obtienen 15 formas distintas de calcular un índice lineal. A continuación se ilustran varias de ellas con kernels que operan sobre vectores:

#include <cuda_runtime.h>
#include <device_launch_parameters.h>

// Bloque 1D, hilo 1D
__global__ void idx_b1d_t1d(int *r, const int *p, const int *q) {
    int gid = threadIdx.x + blockIdx.x * blockDim.x;
    r[gid] = p[gid] + q[gid];
}

// Bloque 1D, hilo 2D
__global__ void idx_b1d_t2d(int *r, const int *p, const int *q) {
    int tid = threadIdx.x + threadIdx.y * blockDim.x;
    int gid = tid + blockIdx.x * blockDim.x * blockDim.y;
    r[gid] = p[gid] + q[gid];
}

// Bloque 1D, hilo 3D
__global__ void idx_b1d_t3d(int *r, const int *p, const int *q) {
    int tid = threadIdx.x
            + threadIdx.y * blockDim.x
            + threadIdx.z * blockDim.x * blockDim.y;
    int gid = tid + blockIdx.x * blockDim.x * blockDim.y * blockDim.z;
    r[gid] = p[gid] + q[gid];
}

// Bloque 2D, hilo 1D
__global__ void idx_b2d_t1d(int *r, const int *p, const int *q) {
    int bid = blockIdx.x + blockIdx.y * gridDim.x;
    int gid = threadIdx.x + bid * blockDim.x;
    r[gid] = p[gid] + q[gid];
}

// Bloque 3D, hilo 1D
__global__ void idx_b3d_t1d(int *r, const int *p, const int *q) {
    int bid = blockIdx.x
            + blockIdx.y * gridDim.x
            + blockIdx.z * gridDim.x * gridDim.y;
    int gid = threadIdx.x + bid * blockDim.x;
    r[gid] = p[gid] + q[gid];
}

// Bloque 2D, hilo 2D
__global__ void idx_b2d_t2d(int *r, const int *p, const int *q) {
    int tid = threadIdx.x + threadIdx.y * blockDim.x;
    int bid = blockIdx.x + blockIdx.y * gridDim.x;
    int gid = tid + bid * blockDim.x * blockDim.y;
    r[gid] = p[gid] + q[gid];
}

// Bloque 2D, hilo 3D
__global__ void idx_b2d_t3d(int *r, const int *p, const int *q) {
    int tid = threadIdx.x
            + threadIdx.y * blockDim.x
            + threadIdx.z * blockDim.x * blockDim.y;
    int bid = blockIdx.x + blockIdx.y * gridDim.x;
    int gid = tid + bid * blockDim.x * blockDim.y * blockDim.z;
    r[gid] = p[gid] + q[gid];
}

// Bloque 3D, hilo 2D
__global__ void idx_b3d_t2d(int *r, const int *p, const int *q) {
    int tid = threadIdx.x + threadIdx.y * blockDim.x;
    int bid = blockIdx.x
            + blockIdx.y * gridDim.x
            + blockIdx.z * gridDim.x * gridDim.y;
    int gid = tid + bid * blockDim.x * blockDim.y;
    r[gid] = p[gid] + q[gid];
}

// Bloque 3D, hilo 3D
__global__ void idx_b3d_t3d(int *r, const int *p, const int *q) {
    int tid = threadIdx.x
            + threadIdx.y * blockDim.x
            + threadIdx.z * blockDim.x * blockDim.y;
    int bid = blockIdx.x
            + blockIdx.y * gridDim.x
            + blockIdx.z * gridDim.x * gridDim.y;
    int gid = tid + bid * blockDim.x * blockDim.y * blockDim.z;
    r[gid] = p[gid] + q[gid];
}

Configuración de grid y bloque para alta ocupación

La ocupación mide la proporción de hilos activos respecto al máximo soportado por un SM. Para ocultar la latencia de memoria es conveniente mantener la mayor ocupación posible.

Sea Tmax el máximo de hilos por SM y Bmax el máximo de bloques por SM. Para alcanzar el 100 % de ocupación se necesita al menos Tmax / Bmax hilos por bloque. Además, el tamaño de bloque debe ser divisor de Tmax, ya que los bloques se asignan atómicamente a un SM.

En arquitecturas comunes:

• Volta, Turing y Ampere anteriores a Ampere GA10x: 2048 / 32 = 64 hilos mínimos.
• Ampere GA10x (RTX 3090, por ejemplo): 1536 / 16 = 96 hilos mínimos.

El máximo común divisor que cumple los requisitos anteriores para estas familias es 512, cuyos divisores mayores o iguales a 96 son 128, 256 y 512. Por tanto, valores robustos para el tamaño de bloque son 128, 256 o 512. El grid debe ser lo suficientemente grande para saturar todos los SM; por ejemplo:

dim3 bloque(256);
dim3 grid((n + bloque.x - 1) / bloque.x);

Ejemplo práctico: multiplicación de matrices

El siguiente kernel realiza el producto de dos matrices cuadradas. Cada hilo calcula un elemento de la matriz resultado:

#include <iostream>
#include <cuda_runtime.h>

constexpr int N = 20;

__global__ void producto_matrices(const float *M, const float *K, float *R, int n) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int fila = blockIdx.y * blockDim.y + threadIdx.y;

    if (col >= n || fila >= n) return;

    float acum = 0.0f;
    for (int k = 0; k < n; ++k) {
        float a = M[fila * n + k];
        float b = K[k * n + col];
        acum += a * b;
    }
    R[fila * n + col] = acum;
}

int main() {
    const int elementos = N * N;
    const size_t bytes = elementos * sizeof(float);

    float *M_host = new float[elementos];
    float *K_host = new float[elementos];
    float *R_host = new float[elementos];

    for (int i = 0; i < elementos; ++i) {
        M_host[i] = 1.0f;
        K_host[i] = 2.0f;
    }

    float *M_dev, *K_dev, *R_dev;
    cudaMalloc((void**)&M_dev, bytes);
    cudaMalloc((void**)&K_dev, bytes);
    cudaMalloc((void**)&R_dev, bytes);

    cudaMemcpy(M_dev, M_host, bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(K_dev, K_host, bytes, cudaMemcpyHostToDevice);

    dim3 bs(16, 16);
    dim3 gs((N + bs.x - 1) / bs.x, (N + bs.y - 1) / bs.y);

    producto_matrices<<<gs, bs>>>(M_dev, K_dev, R_dev, N);
    cudaDeviceSynchronize();

    cudaMemcpy(R_host, R_dev, bytes, cudaMemcpyDeviceToHost);

    for (int i = 0; i < elementos; ++i) {
        std::cout << R_host[i] << "\n";
    }

    cudaFree(M_dev);
    cudaFree(K_dev);
    cudaFree(R_dev);

    delete[] M_host;
    delete[] K_host;
    delete[] R_host;

    return 0;
}

En este ejemplo se usa memoria device explícita con transferencias controladas por el programador. También es posible utilizar cudaMallocManaged para simplificar la gestión de memoria unificada, aunque el control manual suele ofrecer mejor rendimiento.