OpenCL 归约 1

▶ 照着书上的代码，写了几个一步归约的计算，只计算一步，将原数组归约到不超过 1024 个工作项

 ● 代码

// kernel.cl
__kernel void reduce01(__global uint* input, __global uint* output, __local uint* sdata)
{
    const unsigned int tid = get_local_id(0), blockSize = get_local_size(0);
    unsigned int s;

    sdata[tid] = input[get_global_id(0)];
    barrier(CLK_LOCAL_MEM_FENCE);

    // 三种写法，用一种就够
    // 1、模法，问题： % 运算很慢
    for (s = 1; s < blockSize; s <<= 1)
    {
        if (tid % (2 * s) == 0)                 
            sdata[tid] += sdata[tid + s];
        barrier(CLK_LOCAL_MEM_FENCE);
    }
    // 2、间隔缩短法，问题：首次迭代只用一半的工作项，之后每次迭代活跃的工作项持续减少
    for (s = blockSize / 2; s > 0; s >>= 1)
    {
        if (tid < s) 
            sdata[tid] += sdata[tid + s];
        barrier(CLK_LOCAL_MEM_FENCE);
    }
    // 3、间隔增长法，问题：当间隔等于某几个数的时候会产生
    unsigned int index;
    for (s = 1; s < blockSize; s <<= 1)
    {
        if ((index = 2 * s * tid) < blockSize)
            sdata[index] += sdata[index + s];
        barrier(CLK_LOCAL_MEM_FENCE);
    }

    if (tid == 0)
        output[get_group_id(0)] = sdata[0];
}

__kernel void reduce02(__global uint* input, __global uint* output, __local uint* sdata)
{
    const unsigned int tid = get_local_id(0), bid = get_group_id(0), blockSize = get_local_size(0);
    const unsigned int index = bid * (blockSize * 2) + tid;
    unsigned int s;
    
    sdata[tid] = input[index] + input[index + blockSize];// 读入局部内存时就进行一次归约
    barrier(CLK_LOCAL_MEM_FENCE);

    // 两种写法，用一种就够
    // 1、不手动展开循环，仍然有工作项浪费的问题
    for (s = blockSize / 2; s > 0; s >>= 1)
    {
        if (tid < s)
            sdata[tid] += sdata[tid + s];
        barrier(CLK_LOCAL_MEM_FENCE);
    }
    // 2、手动展开最后的循环
    for (s = blockSize / 2; s > 32; s >>= 1)// BUG：如果从 64 开始手工归约，在这一行有且仅有一个工作项会算出 640 = 512 + 128 来，其他行却没问题
    {
        if (tid < s)
            sdata[tid] += sdata[tid + s];
        barrier(CLK_LOCAL_MEM_FENCE);
    }
    if (tid < 32)                           // 手动展开最后的归约，注意同步，书中源代码中没有同步，计算结果是错的
    {
        if (blockSize >= 64)
            sdata[tid] += sdata[tid + 32];
        barrier(CLK_LOCAL_MEM_FENCE);
        if (blockSize >= 32)
            sdata[tid] += sdata[tid + 16];
        barrier(CLK_LOCAL_MEM_FENCE);
        if (blockSize >= 16)
            sdata[tid] += sdata[tid + 8];
        barrier(CLK_LOCAL_MEM_FENCE);
        if (blockSize >= 8)
            sdata[tid] += sdata[tid + 4];
        barrier(CLK_LOCAL_MEM_FENCE);
        if (blockSize >= 4)
            sdata[tid] += sdata[tid + 2];
        barrier(CLK_LOCAL_MEM_FENCE);
        if (blockSize >= 2)
            sdata[tid] += sdata[tid + 1];
        barrier(CLK_LOCAL_MEM_FENCE);
    }

    if (tid == 0)
        output[bid] = sdata[0];
}

// main.c
#include <stdio.h>
#include <stdlib.h>
#include <cl.h>

#define BLOCK_SIZE  256                 // 工作组内最大工作项数为 1024
#define DATA_SIZE   (BLOCK_SIZE * 1024) // 一维最大工作组数为1024

const char *sourceText = "D:/Code/OpenCL/OpenCLProjectTemp/OpenCLProjectTemp/kernel.cl";

int readText(const char* kernelPath, char **pcode)// 读取文本文件放入 pcode，返回字符串长度
{
    FILE *fp;
    int size;
    //printf("<readText> File: %s
", kernelPath);
    fopen_s(&fp, kernelPath, "rb");
    if (!fp)
    {
        printf("Open kernel file failed
");
        getchar();
        exit(-1);
    }
    if (fseek(fp, 0, SEEK_END) != 0)
    {
        printf("Seek end of file failed
");
        getchar();
        exit(-1);
    }
    if ((size = ftell(fp)) < 0)
    {
        printf("Get file position failed
");
        getchar();
        exit(-1);
    }
    rewind(fp);
    if ((*pcode = (char *)malloc(size + 1)) == NULL)
    {
        printf("Allocate space failed
");
        getchar();
        exit(-1);
    }
    fread(*pcode, 1, size, fp);
    (*pcode)[size] = '';
    fclose(fp);
    return size + 1;
}

int main()
{    
    cl_int status;
    cl_uint nPlatform;
    clGetPlatformIDs(0, NULL, &nPlatform);
    cl_platform_id *listPlatform = (cl_platform_id*)malloc(nPlatform * sizeof(cl_platform_id));
    clGetPlatformIDs(nPlatform, listPlatform, NULL);
    cl_uint nDevice;
    clGetDeviceIDs(listPlatform[0], CL_DEVICE_TYPE_ALL, 0, NULL, &nDevice);
    cl_device_id *listDevice = (cl_device_id*)malloc(nDevice * sizeof(cl_device_id));
    clGetDeviceIDs(listPlatform[0], CL_DEVICE_TYPE_ALL, nDevice, listDevice, NULL);
    cl_context context = clCreateContext(NULL, nDevice, listDevice, NULL, NULL, &status);
    cl_command_queue queue = clCreateCommandQueue(context, listDevice[0], CL_QUEUE_PROFILING_ENABLE, &status);

    //const unsigned int nGroup = DATA_SIZE / BLOCK_SIZE;     // reduce01 使用
    const unsigned int nGroup = DATA_SIZE / BLOCK_SIZE / 2; // reduce02 使用
    int *hostA = (cl_int*)malloc(sizeof(cl_int) * DATA_SIZE);
    int *hostB = (cl_int*)malloc(sizeof(cl_int) * nGroup);
    int i;
    unsigned long refSum;
    srand(97);
    for (i = 0, refSum = 0L; i < DATA_SIZE; refSum += (hostA[i++] = 1));// rand()));
    memset(hostB, 0, sizeof(int) * nGroup);
    cl_mem deviceA = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(cl_int) * DATA_SIZE, hostA, &status);
    cl_mem deviceB = clCreateBuffer(context, CL_MEM_WRITE_ONLY, sizeof(cl_int) * nGroup, NULL, &status);
       
    char *code;
    size_t codeLength = readText(sourceText, &code);
    cl_program program = clCreateProgramWithSource(context, 1, (const char**)&code, &codeLength, &status);
    status = clBuildProgram(program, nDevice, listDevice, NULL, NULL, NULL);
    if (status)
    {
        char info[10000];
        clGetProgramBuildInfo(program, listDevice[0], CL_PROGRAM_BUILD_LOG, 10000, info, NULL);
        printf("
%s
", info);
    }
    //cl_kernel kernel = clCreateKernel(program, "reduce01", &status);
    cl_kernel kernel = clCreateKernel(program, "reduce02", &status);   
    
    clSetKernelArg(kernel, 0, sizeof(cl_mem), (void*)&deviceA);
    clSetKernelArg(kernel, 1, sizeof(cl_mem), (void*)&deviceB);
    clSetKernelArg(kernel, 2, BLOCK_SIZE * sizeof(cl_int), NULL);

    size_t globalSize = DATA_SIZE, localSize = BLOCK_SIZE;
    cl_event ev;
    //cl_ulong startTime, endTime;
    status = clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &globalSize, &localSize, 0, NULL, &ev);
    clFinish(queue);
    //clGetEventProfilingInfo(ev, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &startTime, NULL);  // 不启用计时，因为一趟归约时间太短
    //clGetEventProfilingInfo(ev, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &endTime, NULL);
    //printf("Time:%lu.%lu
", (endTime - startTime) / 1000000000, (endTime - startTime) % 1000000000);
    
    clEnqueueReadBuffer(queue, deviceB, CL_TRUE, 0, sizeof(cl_int) * nGroup, hostB, 0, NULL, NULL);
    for (i = 0; i < nGroup; refSum -= hostB[i++]);
    printf("Result %s.
", (refSum == 0) ? "correct" : "error");

    free(hostA);
    free(hostB);
    free(code);
    free(listPlatform);
    free(listDevice);
    clReleaseContext(context);
    clReleaseCommandQueue(queue);
    clReleaseProgram(program);
    clReleaseKernel(kernel);
    clReleaseEvent(ev);
    clReleaseMemObject(deviceA);
    clReleaseMemObject(deviceB);
    getchar();
    return 0;
}

● 输出结果

Result correct.