Este artículo aborda la implementación y optimización de tres problemas clásicos de computación paralela utilizando CUDA. El entorno de desarrollo utilizado consiste en una GPU RTX 3090 con el toolkit CUDA 12.2. Los ejercicios cubren operaciones vectoriales básicas, algoritmos de recorrido (scan) y un problema de renderizado gráfico que requiere una gestión cuidadosa de la memoria y la concurrencia.

Operación SAXPY en GPU

El primer ejercicio consiste en implementar la operación SAXPY (Single-Precision A·X Plus Y). Dados dos vectores X e Y, y una constante scale, se calcula el vector resultado Z tal que Z[i] = scale * X[i] + Y[i]. La implementación directa en CUDA asigna un hilo por elemento, realizando lecturas y escrituras directas en la memoria global.

A continuación se muestra el código del kernel optimizado para esta tarea:

__global__ void saxpy_cuda(int numElements, float scale, float* vecA, float* vecB, float* resultVec) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if (tid < numElements) {
        resultVec[tid] = scale * vecA[tid] + vecB[tid];
    }
}

Esta implementación aprovecha la alta anchura de banda de la memoria global de la GPU para lograr un rendimiento significativamente superior a una implementación en CPU, alcanzando tasas de transferencia efectivas de aproximadamente 484 GB/s en el hardware de prueba.

Suma de Prefijos (Paralel Prefix-Sum)

El segundo bloque se centra en el algoritmo de "scan" o suma de prefijos exclusiva. Para un vector de entrada X, se genera un vector Y donde Y[i] = sum(X[0]...X[i-1]). La implementación se basa en el algoritmo de Blelloch, que consta de dos fases: upsweep (reducción) y downsweep (distribución).

Implementación del Scan Exclusivo

Dado que el algoritmo estándar requiere vectores de longitud potencia de dos, se debe rellenar (padding) el arreglo original si es necesairo. Las funciones auxiliares accumulate_up y propagate_down realizan las fases del algoritmo.

__global__ void accumulate_up(int* data, int len, int stride) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    int jump = stride * 2;
    int idx = tid * jump;
    
    if (idx < len) {
        data[idx + jump - 1] += data[idx + stride - 1];
    }
}

__global__ void propagate_down(int* data, int len, int stride) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    int jump = stride * 2;
    int idx = tid * jump;

    if (idx < len) {
        int temp = data[idx + stride - 1];
        data[idx + stride - 1] = data[idx + jump - 1];
        data[idx + jump - 1] += temp;
    }
}

La función controladora coordina estas fases:

void exclusive_scan(int* d_input, int length, int* d_output) {
    int paddedLen = nextPow2(length);
    int blockSize = 256;
    int gridSize = (paddedLen + blockSize - 1) / blockSize;

    // Fase de reducción (upsweep)
    for (int s = 1; s < paddedLen; s *= 2) {
        accumulate_up<<<gridsize blocksize="">>>(d_output, paddedLen, s);
    }

    // Reiniciar el último elemento
    cudaMemset(&d_output[paddedLen - 1], 0, sizeof(int));

    // Fase de distribución (downsweep)
    for (int s = paddedLen / 2; s >= 1; s /= 2) {
        propagate_down<<<gridsize blocksize="">>>(d_output, paddedLen, s);
    }
}
</gridsize></gridsize>

Detección de Elementos Repetidos (Find Repeats)

Utilizando el scan exclusivo, se puede implementar una función para encontrar elementos adyacentes repetidos en un arreglo. El proceso implica generar un arreglo de banderas (flags), aplicar el scan exclusivo sobre este arreglo para obtener índices compactados y finalmente dispersar los índices de los elementos repetidos al arreglo de salida.

__global__ void detect_duplicates(int* d_input, int len, int* d_flags) {
    extern __shared__ int s_buffer[];
    int gid = blockIdx.x * blockDim.x + threadIdx.x;

    s_buffer[threadIdx.x] = d_input[gid];
    // Cargar el elemento del siguiente bloque para manejo de fronteras
    if (threadIdx.x == blockDim.x - 1) {
        s_buffer[blockDim.x] = d_input[gid + 1];
    }
    __syncthreads();

    if (gid < len) {
        d_flags[gid] = (s_buffer[threadIdx.x] == s_buffer[threadIdx.x + 1]);
    }
}

Renderizador de Círculos

El problema final consiste en renderizar círculos con transparencia, donde el orden de renderizado es crítico para la corrección visual. La optimización de este kernel evoluciona a través de tres etapas principales.

Estrategia Inicial y Problemas de Rendimiento

La primera aproximación asigna un hilo por píxel. Cada hilo itera sobre la lista completa de círculos para determinar si afecta su píxel. Esto genera un bajo randimiento debido a que muchos círculos no intersectan el bloque de píxeles procesado por un bloque de hilos, resultando en cálculos redundantes y acceso ineficiente a memoria.

Optimización mediante Descarte y Scan en Memoria Compartida

Para mejorar la eficiencia, se implementa una fase de "culling" o descarte. Cada bloque de hilos determina qué círculos intersectan su región de la imagen. Se utilizan operaciones de scan en memoria compartida para compactar la lista de círculos activos para ese bloque.

El siguiente paso optimiza aún más el rendimiento mediante tres técnicas clave:

1. Uso de la función circleInBoxConservative para pruebas de intersección eficientes.
2. Agrupación de círculos activos en un búfer de memoria compartida hasta llenarlo antes de procesarlos, reduciendo la sobrecarga de sincronización.
3. Empaquetamiento de datos (posición y radio) en estructuras float4 para optimizar las transacciones de memoria.

La implementación final del kernel demuestra un rendimiento notablemente superior:

const int BUF_SIZE = SCAN_BLOCK_DIM;
__shared__ float4 s_pos_radius[BUF_SIZE];
__shared__ float3 s_color[BUF_SIZE];
__shared__ uint s_scan_temp[BUF_SIZE * 2];
__shared__ uint s_flags[BUF_SIZE];
__shared__ uint s_indices[BUF_SIZE];

__global__ void render_circles_cuda() {
    float4 pixelColor;
    int tid = threadIdx.y * blockDim.x + threadIdx.x;
    int pixelX = blockIdx.x * blockDim.x + threadIdx.x;
    int pixelY = blockIdx.y * blockDim.y + threadIdx.y;

    // Coordenadas del bloque
    int bx1 = blockIdx.x * blockDim.x, bx2 = bx1 + blockDim.x;
    int by1 = blockIdx.y * blockDim.y, by2 = by1 + blockDim.y;

    int imgW = cuConstRendererParams.imageWidth;
    int imgH = cuConstRendererParams.imageHeight;
    float invW = 1.0f / imgW, invH = 1.0f / imgH;

    // Cargar color actual a registro
    pixelColor = __ldcs((float4*)cuConstRendererParams.imageData + pixelY * imgW + pixelX);

    int activeCount = 0;

    for (int circleIdx = 0; circleIdx < cuConstRendererParams.numCircles; circleIdx += BUF_SIZE) {
        int batch = min(BUF_SIZE, cuConstRendererParams.numCircles - circleIdx);
        int isRelevant = 0;

        float3 pos = *((float3*)cuConstRendererParams.position + circleIdx + tid);
        float radius = *(cuConstRendererParams.radius + circleIdx + tid);
        float3 col = *((float3*)cuConstRendererParams.color + circleIdx + tid);

        if (tid < batch) {
            isRelevant = circleInBoxConservative(pos.x, pos.y, radius, bx1 * invW, bx2 * invW, by2 * invH, by1 * invH);
        }

        if (tid < BUF_SIZE) {
            s_flags[tid] = isRelevant;
            s_indices[tid] = isRelevant;
        }
        __syncthreads();

        sharedMemExclusiveScan(tid, s_flags, s_indices, s_scan_temp, nextPow2(batch));

        int numActive = s_indices[nextPow2(batch) - 1] + s_flags[nextPow2(batch) - 1];

        // Si el búfer está lleno, procesar los círculos almacenados
        if (activeCount + numActive > BUF_SIZE) {
            for (int i = 0; i < activeCount; ++i) {
                float4 pr = s_pos_radius[i];
                float2 normPos = make_float2(invW * (pixelX + 0.5f), invH * (pixelY + 0.5f));
                shadePixel(i, normPos, pr, &pixelColor);
            }
            activeCount = 0;
        }

        __syncthreads();
        if (s_flags[tid] == 1 && tid < batch) {
            int writeIdx = s_indices[tid] + activeCount;
            s_pos_radius[writeIdx] = make_float4(pos.x, pos.y, pos.z, radius);
            s_color[writeIdx] = col;
        }
        activeCount += numActive;
    }

    // Procesar restantes
    if (activeCount > 0) {
        __syncthreads();
        for (int i = 0; i < activeCount; ++i) {
            float4 pr = s_pos_radius[i];
            float2 normPos = make_float2(invW * (pixelX + 0.5f), invH * (pixelY + 0.5f));
            shadePixel(i, normPos, pr, &pixelColor);
        }
    }

    __stwt((float4*)cuConstRendererParams.imageData + pixelY * imgW + pixelX, pixelColor);
}

Esta versión final logra cumplir con los requisitos de rendimiento en escenas complejas como rand100k y snowsingle, optimizando tanto el acceso a memoria como la lógica de cómputo dentro del kernel.