数组规约加法（任意长度）

▶ 采用两种方法计算任意长度的数组规约加法。第一种是将原数组收缩到一个固定长度的共享内存数组中，再利用二分规约计算求和，这里测试了固定长度分别取2，4，8……1024的情形；第二种方法是将原数组分段，每段使用各自的共享内存进行二分规约，结果汇总后再次分段，规约，直至数组长度减小到1，输出结果。

● 代码

#include <stdio.h>
#include <stdlib.h>
#include <windows.h>
#include <time.h>
#include <math.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#define ARRAY_SIZE      (1024 * 7 + 119)
#define WIDTH           512
#define SEED            1   //(unsigned int)clock()
#define CEIL(x,y)       (int)(( x - 1 ) /  y + 1)    
#define MIN(x,y)        ((x)<(y)?(x):(y))
//#define RESIZE(x)       MIN(WIDTH, 1 << (int)ceil(log(x) / log(2)))

typedef int format;    // int or float

void checkCudaError(cudaError input)
{
    if (input != cudaSuccess)
    {
        printf("
	find a cudaError!");
        exit(1);
    }
    return;
}

// 不超过1024个元素的，任意个元素的规约加法。调用时指定与原数组长度相同大小的共享内存数组
__global__ void add_1(const format *in, format *out)
{
    extern __shared__ format sdata[];
    sdata[threadIdx.x] = in[threadIdx.x];

    for (int s = blockDim.x; s > 0; s = (s - 1) / 2 + 1)
    {        
        sdata[threadIdx.x] += (threadIdx.x >= s / 2) ? 0 : sdata[threadIdx.x + (s - 1) / 2 + 1];// 考虑加数是否为奇数个，后半线程都加 0
        __syncthreads();
    }

    if (threadIdx.x == 0)
        *out = sdata[0];
    return;
}

// 将原数组收缩到blockDim.x的尺度上进行二分规约。调用时指定与线程数相同大小的共享内存数组
__global__ void add_2(const format *in, format *out, int size)// size为原数组大小，可以远大于blockDim.x
{
    extern __shared__ format sdata[];
    sdata[threadIdx.x] = (threadIdx.x < size) ? in[threadIdx.x] : 0;

    // 收尾操作，在另外两个规约加法函数中没有的
    for (int id = threadIdx.x + blockDim.x; id < size; id += blockDim.x)
        sdata[threadIdx.x] += in[id];
    __syncthreads();

    for (int jump = blockDim.x / 2; jump > 0; jump /= 2)
    {
        sdata[threadIdx.x] += (threadIdx.x >= jump) ? 0 : sdata[threadIdx.x + jump];
        __syncthreads();
    }
    __syncthreads();

    if (threadIdx.x == 0)
        *out = sdata[0];
    return;
}
// 限制元素个数为2的整数次幂的规约加法。调用时指定共享内存数组大小为 1 << ((log(size) - 1) / 2 + 1) 
__global__ void add_simple(const format *in, format *out, int size)
{
    extern __shared__ format sdata[];
    sdata[threadIdx.x] = (blockIdx.x * blockDim.x + threadIdx.x < size) ? in[blockIdx.x * blockDim.x + threadIdx.x] : 0;

    for (int jump = blockDim.x / 2; jump > 0; jump /= 2)
    {
        sdata[threadIdx.x] += (threadIdx.x < jump) ? sdata[threadIdx.x + jump] : 0;
        __syncthreads();
    }
    __syncthreads();

    if (threadIdx.x == 0)
        out[blockIdx.x] = sdata[0];
    return;
}

int main()
{
    int i, j, leftLength;
    format h_in[ARRAY_SIZE], h_out[10 + 1], cpu_sum;
    format *d_in, *d_out, *d_temp1, *d_temp2;
    cudaEvent_t start, stop;
    float elapsedTime[10 + 1];
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    checkCudaError(cudaMalloc((void **)&d_in, sizeof(format) * ARRAY_SIZE));
    checkCudaError(cudaMalloc((void **)&d_out, sizeof(format)));
    checkCudaError(cudaMalloc((void **)&d_temp1, sizeof(format) * CEIL(ARRAY_SIZE, WIDTH)));
    checkCudaError(cudaMalloc((void **)&d_temp2, sizeof(format) * CEIL(CEIL(ARRAY_SIZE, WIDTH), WIDTH)));

    srand(SEED);
    for (i = 0, cpu_sum = 0; i < ARRAY_SIZE; i++)
    {
        h_in[i] = rand() - RAND_MAX / 2;
        cpu_sum += h_in[i];
    }
    cudaMemcpy(d_in, h_in, sizeof(format) * ARRAY_SIZE, cudaMemcpyHostToDevice);

    // 方法一
    for (i = 0; i < 10; i++)
    {
        cudaEventRecord(start, 0);
        add_2 << < 1, 2 << i, sizeof(format) * (2 << i) >> > (d_in, d_out, ARRAY_SIZE);
        cudaMemcpy(&h_out[i], d_out, sizeof(format), cudaMemcpyDeviceToHost);
        cudaDeviceSynchronize();
        cudaEventRecord(stop, 0);
        cudaEventSynchronize(stop);
        cudaEventElapsedTime(&elapsedTime[i], start, stop);
    }
    // 方法二
    cudaEventRecord(start, 0);
    for (i = 0; i < CEIL(log(ARRAY_SIZE), log(WIDTH)); i++)// 一次循环完成一次迭代，共需要 CEIL(log(ARRAY_SIZE), log(WIDTH))次迭代
    {
        if (i == 0)// 首次迭代，将d_in算到d_temp1中去
        {
            add_simple << <CEIL(ARRAY_SIZE, WIDTH), WIDTH, sizeof(format) * WIDTH >> > (d_in, d_temp1, ARRAY_SIZE);
            cudaDeviceSynchronize();
            leftLength = ARRAY_SIZE / WIDTH + (bool)(ARRAY_SIZE % WIDTH);
        }
        else if (i % 2)// 非首次的偶数次，把d_temp1算到d_temp2中去
        {
            add_simple << <CEIL(leftLength, WIDTH), WIDTH, sizeof(format) * WIDTH >> > (d_temp1, d_temp2, WIDTH);
            cudaDeviceSynchronize();
            leftLength = leftLength / WIDTH + (bool)(leftLength % WIDTH);
        }
        else// 非首次的奇数次，把d_temp2算到d_temp1中去
        {
            add_simple << <CEIL(leftLength, WIDTH), WIDTH, sizeof(format) * WIDTH >> > (d_temp2, d_temp1, WIDTH);
            cudaDeviceSynchronize();
            leftLength = leftLength / WIDTH + (bool)(leftLength % WIDTH);
        }
    }
    if (i % 2)//说明经历了i次规约，i奇结果在d_temp1中，i偶结果在d_temp2中
        cudaMemcpy(&h_out[10], d_temp1, sizeof(format), cudaMemcpyDeviceToHost);
    else
        cudaMemcpy(&h_out[10], d_temp2, sizeof(format), cudaMemcpyDeviceToHost);
    cudaDeviceSynchronize();
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime[10], start, stop);

    printf("
	CPU:			%f", (float)cpu_sum);
    for (i = 0; i < 10; i++)
        printf("
	GPU, %4d threads:	%d	Time:	%6.2f ms", 2 << i, (int)h_out[i], elapsedTime[i]);
    printf("
	GPU,   new method:	%d	Time:	%6.2f ms", (int)h_out[10], elapsedTime[10]);

    cudaFree(d_in);
    cudaFree(d_out);
    cudaFree(d_temp1);
    cudaFree(d_temp2);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);

    getchar();
    return 0;
}

● 输出结果，可见随着一个线程块内的线程数的增大，每次全局内存读取的带宽增大，计算所需迭代次数减小，但共享内存结果的汇总耗时增加，总时间先减小后增大，在blockDim.x在256和512左右接近最优。采用第二种方法在CPU中逻辑较为复杂，没有去的较好的改进。