GPU编程 CUDA C++ 数组归约的示例和解释

数组归约是指将一个数组中的所有元素经过某种操作后，得到一个最终结果的过程。例如，将一个数组中的所有元素相加，就是一种数组归约操作。在CUDA中，可以使用reduce函数来实现数组归约。

示例代码如下：

#include <stdio.h>
#define N 1024

__global__ void reduce(int *g_idata, int *g_odata) {
    extern __shared__ int sdata[];

    // 每个线程加载一个元素到共享内存
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
    sdata[tid] = g_idata[i];
    __syncthreads();

    // 归约操作
    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }

    // 将归约结果存储到全局内存中
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}

int main(void) {
    int *a, *d_a, *d_b;
    int size = N * sizeof(int);

    // 分配内存空间
    a = (int *)malloc(size);
    cudaMalloc((void **)&d_a, size);
    cudaMalloc((void **)&d_b, size);

    // 初始化数组
    for (int i = 0; i < N; i++) {
        a[i] = i;
    }

    // 将数组复制到设备上
    cudaMemcpy(d_a, a, size, cudaMemcpyHostToDevice);

    // 归约操作
    int block_size = 512;
    reduce<<<(N + block_size - 1) / block_size, block_size, block_size * sizeof(int)>>>(d_a, d_b);

    // 将结果从设备上复制回主机内存
    int result;
    cudaMemcpy(&result, d_b, sizeof(int), cudaMemcpyDeviceToHost);

    printf("sum: %d\n", result);

    // 释放内存空间
    free(a);
    cudaFree(d_a);
    cudaFree(d_b);

    return 0;
}

在上面的示例代码中，首先定义了一个大小为N的整型数组a，然后将该数组复制到设备上。接着定义了一个reduce函数，该函数使用共享内存实现了数组归约操作。最后，在主函数中调用reduce函数进行归约操作，并将结果从设备上复制回主机内存。最终，输出结果即为数组中所有元素的和。

---

或者：

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>

__global__ void reduce(int *input, int *output, int size) {
    extern __shared__ int sdata[];
    int tid = threadIdx.x;
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    sdata[tid] = i < size ? input[i] : 0;
    __syncthreads();
    for (int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    if (tid == 0) {
        output[blockIdx.x] = sdata[0];
    }
}

int main() {
    const int size = 100000;
    const int threads_per_block = 256;
    const int blocks = (size + threads_per_block - 1) / threads_per_block;
    int *input = (int*)malloc(size * sizeof(int));
    for (int i = 0; i < size; i++) {
        input[i] = i;
    }
    int *d_input, *d_output;
    cudaMalloc(&d_input, size * sizeof(int));
    cudaMalloc(&d_output, blocks * sizeof(int));
    cudaMemcpy(d_input, input, size * sizeof(int), cudaMemcpyHostToDevice);
    reduce<<<blocks, threads_per_block, threads_per_block * sizeof(int)>>>(d_input, d_output, size);
    int *output = (int*)malloc(blocks * sizeof(int));
    cudaMemcpy(output, d_output, blocks * sizeof(int), cudaMemcpyDeviceToHost);
    int sum = 0;
    for (int i = 0; i < blocks; i++) {
        sum += output[i];
    }
    printf("sum = %d\n", sum);
    free(input);
    free(output);
    cudaFree(d_input);
    cudaFree(d_output);
    return 0;
}

首先，我们定义了一个名为“reduce”的CUDA内核函数。该函数接受三个参数：输入数组、输出数组和数组大小。内核中的第一行代码使用“extern shared”关键字定义了一个共享内存数组，“sdata”。这个数组的大小等于每个线程块中的线程数（即“threads_per_block”）乘以每个整数的大小（即“sizeof(int)”）。

接下来，我们获取当前线程的ID（即“tid”）和输入数组的索引（即“i”）。我们使用输入索引来访问输入数组的元素，并将它们存储在共享内存数组中。如果输入索引超出了数组大小，则将共享内存中的元素设置为零。

然后，我们调用“__syncthreads()”函数来同步所有线程的执行。这是因为我们需要确保所有线程都已将其输入值存储在共享内存中，才能开始归约操作。

接下来，我们使用一个循环来执行归约操作。在每次循环迭代中，我们将“s”设置为当前线程块的线程数的一半。然后，如果当前线程ID小于“s”，则将当前线程的共享内存值与其相邻的线程的共享内存值相加。这个过程会重复执行，直到“s”为零为止。

最后，如果当前线程的ID为零，则将归约结果存储在输出数组中。这个过程对于每个线程块都会执行一次，因此输出数组的大小等于线程块的数量（即“blocks”）。

在主函数中，我们首先定义了输入数组的大小（即“size”）、每个线程块中的线程数（即“threads_per_block”）和线程块的数量（即“blocks”）。我们还使用标准库函数“malloc”分配了输入和输出数组的内存，并将输入数组初始化为顺序整数。

接下来，我们使用CUDA函数“cudaMalloc”分配了设备内存，并使用“cudaMemcpy”函数将输入数组从主机内存复制到设备内存。

然后，我们调用“reduce”内核函数，并传递输入数组、输出数组和数组大小作为参数。我们还传递了共享内存大小作为第三个参数，这个大小等于每个线程块中的线程数乘以每个整数的大小。

最后，我们使用“cudaMemcpy”函数将输出数组从设备内存复制到主机内存，并计算输出数组中的所有元素的和。我们释放了所有内存，然后返回零。

在此示例中，我们使用了一些简单的技巧来优化内核函数的性能，例如使用共享内存来提高内存访问效率，并使用循环执行归约操作。这些技巧可用于各种不同的CUDA应用程序中，以提高其性能和效率。