parallel/reduction.cu

配套源码，运行方法见同目录 README。返回实验总览。

// Optional original CUDA example. Not compiled/run on the macOS CPU host.
// 256 threads per block; one input value per thread, padded with zero.
#include <cuda_runtime.h>
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <vector>
#define CUDA_CHECK(call) do { cudaError_t err = (call); if (err != cudaSuccess) { \
    std::fprintf(stderr, "%s:%d: %s\n", __FILE__, __LINE__, cudaGetErrorString(err)); std::exit(1); } } while (0)

__global__ void block_sum(const float *input, float *partials, int n) {
    __shared__ float scratch[256];
    unsigned lane = threadIdx.x;
    unsigned i = blockIdx.x * blockDim.x + lane;
    scratch[lane] = i < static_cast<unsigned>(n) ? input[i] : 0.0f;
    __syncthreads();
    for (unsigned stride = blockDim.x / 2; stride; stride /= 2) {
        if (lane < stride) scratch[lane] += scratch[lane + stride];
        __syncthreads(); // ALL threads reach each barrier, including inactive lanes.
    }
    if (lane == 0) partials[blockIdx.x] = scratch[0];
}

int main() {
    constexpr int n = 100003; // Tail not a multiple of the block size.
    constexpr int blocks = (n + 255) / 256;
    std::vector<float> input(n), partials(blocks);
    double expected = 0;
    for (int i = 0; i < n; ++i) { input[i] = float(i % 17) / 16.0f; expected += input[i]; }
    float *device_input = nullptr, *device_partials = nullptr;
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void **>(&device_input), input.size() * sizeof(float)));
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void **>(&device_partials), partials.size() * sizeof(float)));
    CUDA_CHECK(cudaMemcpy(device_input, input.data(), input.size() * sizeof(float), cudaMemcpyHostToDevice));
    block_sum<<<blocks, 256>>>(device_input, device_partials, n);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());
    CUDA_CHECK(cudaMemcpy(partials.data(), device_partials, partials.size() * sizeof(float), cudaMemcpyDeviceToHost));
    double result = 0;
    for (float part : partials) result += part;
    CUDA_CHECK(cudaFree(device_input));
    CUDA_CHECK(cudaFree(device_partials));
    double error = std::fabs(result - expected);
    if (error > 1e-5 * std::max(1.0, std::fabs(expected))) { std::fprintf(stderr, "reduction mismatch\n"); return 1; }
    std::printf("reduction: result=%.9f expected=%.9f error=%.9g\n", result, expected, error);
}