CMU15418-Assignment2-解析

这个作业有三个部分, 都是 CUDA 编程. 前两个比较简单, 最后一个比较难.

本文的运行环境:

RTX 3090
CUDA 12.2

Part 1: SAXPY

用 CUDA 实现一个在 GPU 上运行的 SAXPY 程序. 输入两个数组 X, Y 以及一个常数 scale, 输出一个数组 Z, 并且 Z[i] = scale * X[i] + Y[i].

直接访问和存储 global memory 就行.

__global__ void
saxpy_kernel(int N, float alpha, float* x, float* y, float* result) {
    int index = blockIdx.x * blockDim.x + threadIdx.x;

    if (index < N)
       result[index] = alpha * x[index] + y[index];
}

运行结果如下. 相比于 Assignment 1 中的 CPU 实现有十几倍的提升.

Overall time: 0.461 ms        [484.865 GB/s]

part 2: Parallel Prefix-Sum

这个 part 又分成两部分. 一部分是 Exclusive Prefix Sum 的 CUDA 实现, 另一部分是在此基础上实现的 Find Repeats

Exclusive Prefix Sum

这部分需要实现一个求数组前缀和的程序, 该前缀和不包括当前元素. 即输入数组 X, 输出数组 Y, 其中 Y[i] = sum(X[j] for j in [0, i-1]). 作业中给了 CPU 上的并行实现, 这里只要把它移植到 CUDA 上就可以. CPU 上的实现分为三个部分:

upsweep.
设置结果数组的最后一位为 0.
downsweep.

作业中提供的 CPU 算法只适用于长度为 2 的整数幂的数组. 当输入数组长度不满足要求时, 在 GPU 上申请足够大的数组 (大于等于数组长度, 并且长度是 2 的整数幂), 把数组复制到这段数组上, 末尾置 0. 在这段数组上执行 exclusive scan, 把结果中的对应部分复制回 CPU.

首先是三个辅助函数.

__global__ void upsweep(int* device_result, int N, int twod) {
    size_t id = blockDim.x * blockIdx.x + threadIdx.x;
    size_t twod1 = twod * 2;
    size_t i = id * twod1;
    if (i < N)
        device_result[i + twod1 - 1] += device_result[i + twod - 1];
}

__global__ void downsweep(int* device_result, int N, int twod) {
    size_t id = blockDim.x * blockIdx.x + threadIdx.x;
    size_t twod1 = twod * 2;
    size_t i = id * twod1;
    if (i < N) {
        int t = device_result[i + twod - 1];
        device_result[i + twod - 1] = device_result[i + twod1 - 1];
        device_result[i + twod1 - 1] += t;
    }
}

__global__ void set(int* device_result, int N, int idx, int value) {
    if (idx < N) device_result[idx] = value;
}

接着实现 exclusive scan 函数.

void exclusive_scan(int* device_start, int length, int* device_result)
{
    length = nextPow2(length);
    const int blockDim = 256;
    const int gridDim  = (length + blockDim - 1) / blockDim;

    for (int twod = 1; twod < length; twod *= 2) {
        upsweep<<<gridDim, blockDim>>>(device_result, length, twod);
    }

    set<<<1, 1>>>(device_result, length, length - 1, 0);

    for (int twod = length / 2; twod >= 1; twod /= 2) {
        downsweep<<<gridDim, blockDim>>>(device_result, length, twod);
    }
}

运行测试程序查看得分.

# 测试程序
make check_scan

# 测试结果
-------------------------------------------------------------------------
| Element Count   | Fast Time       | Your Time       | Score           |
-------------------------------------------------------------------------
| 10000           | 0.190           | 0.092           | 1.25            |
| 100000          | 0.325           | 0.127           | 1.25            |
| 1000000         | 0.481           | 0.293           | 1.25            |
| 2000000         | 0.725           | 0.522           | 1.25            |
-------------------------------------------------------------------------
|                                   | Total score:    | 5/5             |
-------------------------------------------------------------------------

Find Repeats

Find Repeats 是查找数组中的相邻重复元素, 并把重复位置的序号按照先后顺序放到另一个数组中. 输入数组 X, 输出数组 Y, 对于 y = Y[i] 有 X[y] == X[y + 1].

可以通过 Exclusive Prefix Sum 实现 Find Repeats. 首先在数组 X 中查找重复元素得到重复数组 flag, 若 X[i] == X[i+1], 则 flag[i] = 1. 接着通过 Exclusive Prefix Sum 求 flag 的前缀和得到位置数组 idx. 若 flag[i] == 1, 则 idx[i] 表示了该元素在输出数组中的位置. 最后把结果输出到数组 Y, 若 flag[i] == 1, 则 Y[idx[i]] = i.

代码中 findAdjacentSame 实现查找重复元素和函数. 这里通过 shared memory 加速该函数.

__global__ void findAdjacentSame(int* device_input, int length, int* device_output) {
    extern __shared__ int smem[];

    size_t global_idx = blockDim.x * blockIdx.x + threadIdx.x;

    smem[threadIdx.x] = device_input[global_idx];
    if (threadIdx.x == blockDim.x - 1)
        smem[blockDim.x] = device_input[global_idx + 1];
    __syncthreads();

    device_output[global_idx] = (smem[threadIdx.x] == smem[threadIdx.x + 1]);
}

__global__ void setRepeats(int* device_idxs, int* device_same, int length, int* device_output) {
    size_t id = blockDim.x * blockIdx.x + threadIdx.x;
    if (device_same[id] == 1) {
        device_output[device_idxs[id]] = id;
    }
}

int find_repeats(int *device_input, int length, int *device_output) {
    // device_same 就是上述的数组 flag
    // device_idxs 就是上述的数组 idx
    int *device_same;
    int *device_idxs;
    cudaMalloc(&device_same, nextPow2(length) * sizeof(int));
    cudaMalloc(&device_idxs, nextPow2(length) * sizeof(int));

    size_t blockDim = 256;
    size_t gridDim = (length + blockDim - 1) / blockDim;

    // 查找重复元素
    findAdjacentSame<<<gridDim, blockDim, (blockDim + 1) * sizeof(int)>>>(device_input, length, device_same);

    cudaMemcpy(device_idxs, device_same, nextPow2(length) * sizeof(int), cudaMemcpyDeviceToDevice);
    // 对 flag 求前缀和
    exclusive_scan(device_input, length, device_idxs);

    // 输出结果
    setRepeats<<<gridDim, blockDim>>>(device_idxs, device_same, length, device_output);

    int count;
    cudaMemcpy(&count, device_idxs + length - 1, sizeof(int), cudaMemcpyDeviceToHost);

    cudaFree(device_same);
    cudaFree(device_idxs);

    return count;
}

运行测试程序查看得分.

# 测试程序
make check_find_repeats

# 测试结果
-------------------------------------------------------------------------
| Element Count   | Fast Time       | Your Time       | Score           |
-------------------------------------------------------------------------
| 10000           | 0.227           | 0.087           | 1.25            |
| 100000          | 0.384           | 0.142           | 1.25            |
| 1000000         | 0.815           | 0.615           | 1.25            |
| 2000000         | 1.243           | 0.942           | 1.25            |
-------------------------------------------------------------------------
|                                   | Total score:    | 5/5             |
-------------------------------------------------------------------------

Part 3: A Simple Circle Renderer

这一部分要实现一个渲染器, 作业代码中有实现可用的 RefRenderer, 为渲染器的具体功能提供参考. 这个渲染器主要用来渲染圆, 由于不同的圆之间会有重叠, 不同的计算顺序会导致重叠部分不同. CUDA 上的渲染器首先应该保证功能正确, 即不同的圆在渲染的时候存在先后关系, 必须按照给定的顺序渲染重叠部分. 这部分主要就是修改文件 cudaRenderer.cu 中的函数 kernelRenderCircles. 我也是经过了多次优化, 才使得 cudaRenderer 的性能勉强达到要求.

最开始的版本性能很低, 仅能正确渲染图像. 我考虑每个线程处理一个像素, 按照顺序渲染每个圆. 每个线程运行类似的代码, 处理不同的数据, 任务之间没有相互关系. 这比较符合这门课的主题. 就有了下面的代码.

// shared memory
const int smem_size = SCAN_BLOCK_DIM;
__shared__ float smem_ps[smem_size * 3];
__shared__ float smem_rads[smem_size];

__global__ void kernelRenderCircles() {

    int pixelX = blockIdx.x * blockDim.x + threadIdx.x;
    int pixelY = blockIdx.y * blockDim.y + threadIdx.y;

    for (int index = 0; index < cuConstRendererParams.numCircles; ++index) {
        int index3 = 3 * index;

        // 这里所有线程从 global memory 中获取相同的数据, 带宽利用率极低, 因此可以考虑使用 shared memory 加速
        float3 p = *(float3*)(&cuConstRendererParams.position[index3]);
        float  rad = cuConstRendererParams.radius[index];

        short imageWidth = cuConstRendererParams.imageWidth;
        short imageHeight = cuConstRendererParams.imageHeight;
        short minX = static_cast<short>(imageWidth * (p.x - rad));
        short maxX = static_cast<short>(imageWidth * (p.x + rad)) + 1;
        short minY = static_cast<short>(imageHeight * (p.y - rad));
        short maxY = static_cast<short>(imageHeight * (p.y + rad)) + 1;

        short screenMinX = (minX > 0) ? ((minX < imageWidth) ? minX : imageWidth) : 0;
        short screenMaxX = (maxX > 0) ? ((maxX < imageWidth) ? maxX : imageWidth) : 0;
        short screenMinY = (minY > 0) ? ((minY < imageHeight) ? minY : imageHeight) : 0;
        short screenMaxY = (maxY > 0) ? ((maxY < imageHeight) ? maxY : imageHeight) : 0;

        float invWidth = 1.f / imageWidth;
        float invHeight = 1.f / imageHeight;

        if (pixelX < screenMinX || pixelX >= screenMaxX) continue;
        if (pixelY < screenMinY || pixelY >= screenMaxY) continue;

        float4* imgPtr = (float4*)(&cuConstRendererParams.imageData[4 * (pixelY * imageWidth + pixelX)]);
        float2 pixelCenterNorm = make_float2(invWidth * (static_cast<float>(pixelX) + 0.5f),
                                             invHeight * (static_cast<float>(pixelY) + 0.5f));

        // imgPtr 位于 global memory 中, shadePixel 修改其内容使得性能低下
        shadePixel(index, pixelCenterNorm, p, imgPtr);
    }
}

static inline
size_t didUp(size_t a, size_t b) {
    return (a + b - 1) / b;
}

void
CudaRenderer::render() {

    // 256 threads per block is a healthy number
    dim3 blockDim(16, 16);
    dim3 gridDim(didUp(image->width, blockDim.x), didUp(image->height, blockDim.y));

    kernelRenderCircles<<<gridDim, blockDim>>>();
    cudaDeviceSynchronize();
}

上面这个版本的代码在第 10 行和 11 行处获取 global memory 的内容, 但是一个 warp 内的所有线程访问相同的地址. 我认为这会导致带宽利用率低, 因此想通过 shared memory 改进这一块性能. 另一方面, 函数 shadePixel 中会修改 imgPtr 的内容, 会使得线程反复访问 global memory, 降低了函数的性能. 因为每个线程处理一个像素, 所以可以通过把 imgPtr 对应的数据读取到寄存器避免反复访问 global memory.

但是实际修改了之后发现性能提升并不明显. 经过分析, 程序运行 rgb 渲染的时候性能还算可以, 但是运行 rand10k 和 rand100k 的时候性能极低. 这样的性能降低肯定是因为中间处理各个像素导致的, 而不是两端的 global memory 访存. 这是因为 rgb, rand10k 和 rand100k 两端的 global memory 访存相同. 因此进一步的优化需要减少 for 循环中处理像素的时间.

一个 block 有 256 个线程, 渲染完整中一块大小为 16 * 16 部分的图像. 这是比较小的, 因此大部分的圆实际上都不会覆盖这一块图像. 也就说对于一个 block 来讲, 大部分圆是不需要渲染的. 但是在 for 循环内部, 256 个线程会同时计算同一个圆是否需要渲染, 导致了相同的计算重复 256 次, 降低了程序性能.

为了改进程序性能, 应该每个线程读取一个圆, 判断这个圆是否与该 block 对应的图像相交, 如果相交就留待计算. 具体的步骤是设置一个数组 need 表示某个圆是否与这块图像相交. need[i] = 1 表示第 i 个圆与图像相交. 接着再通过 part 2 中的 Exclusive Prefix Sum 对数组 need 求和得到数组 idx. 当 need[i] 为 1 时, idx[i] 就表示需要处理的圆的序号. 这把这些序号都整理到结果数组 Y 中, Y 中的每一个值都代表了该 block 需要渲染的圆的序号. 接着就可以完成实际渲染的任务了. 这样就有了下面这个版本的代码. part 2 是在 global memory 上求 Exclusive Prefix Sum, 这里需要在 shared memory 上求. 在文件 exclusiveScan.cu_inl 给出了求 shared memory 的例程 sharedMemExclusiveScan, 可以直接调用.

// shared memory
const int smem_size = SCAN_BLOCK_DIM;
__shared__ float smem_ps[smem_size * 3];
__shared__ float smem_colors[smem_size * 3];
__shared__ float smem_rads[smem_size];
__shared__ uint smem_idxs[smem_size * 2];
__shared__ uint smem_needs[smem_size];
__shared__ uint smem_sums[smem_size];

__global__ void kernelRenderCircles() {
    float4 img;

    const int threadId = threadIdx.y * blockDim.x + threadIdx.x;

    const unsigned int globalIdX = blockIdx.x * blockDim.x + threadIdx.x;
    const unsigned int globalIdY = blockIdx.y * blockDim.y + threadIdx.y;

    const unsigned int x1 = blockIdx.x * blockDim.x;
    const unsigned int x2 = x1 + blockDim.x;
    const unsigned int y1 = blockIdx.y * blockDim.y;
    const unsigned int y2 = y1 + blockDim.y;

    short imageWidth = cuConstRendererParams.imageWidth;
    short imageHeight = cuConstRendererParams.imageHeight;

    // 读取 imageData 的数据到寄存器
    img = *((float4*)(cuConstRendererParams.imageData) + globalIdY * imageWidth + globalIdX);

    for (int index = 0; index < cuConstRendererParams.numCircles; index += smem_size) {
        int index3 = 3 * index;
        int size = min(smem_size, cuConstRendererParams.numCircles - index);

        // 先把 global memory 中的数据加载到 shared memory
        __syncthreads();
        if (threadId < size) {
            smem_ps[threadId] = __ldca(&cuConstRendererParams.position[index3 + threadId]);
            smem_ps[threadId + size] = __ldca(&cuConstRendererParams.position[index3 + threadId + size]);
            smem_ps[threadId + size * 2] = __ldca(&cuConstRendererParams.position[index3 + threadId + size * 2]);
            smem_rads[threadId] = __ldca(&cuConstRendererParams.radius[index + threadId]);
            smem_colors[threadId] = __ldca(&cuConstRendererParams.color[index3 + threadId]);
            smem_colors[threadId + size] = __ldca(&cuConstRendererParams.color[index3 + threadId + size]);
            smem_colors[threadId + size * 2] = __ldca(&cuConstRendererParams.color[index3 + threadId + size * 2]);
        }
        __syncthreads();

        int need_cal = 0;

        // 判断哪些圆与本 block 对应的图像相交
        if (threadId < size) {

            float3 p = ((float3*)smem_ps)[threadId];
            float rad = smem_rads[threadId];

            short minX = static_cast<short>(imageWidth * (p.x - rad));
            short maxX = static_cast<short>(imageWidth * (p.x + rad)) + 1;
            short minY = static_cast<short>(imageHeight * (p.y - rad));
            short maxY = static_cast<short>(imageHeight * (p.y + rad)) + 1;

            short screenMinX = (minX > 0) ? ((minX < imageWidth) ? minX : imageWidth) : 0;
            short screenMaxX = (maxX > 0) ? ((maxX < imageWidth) ? maxX : imageWidth) : 0;
            short screenMinY = (minY > 0) ? ((minY < imageHeight) ? minY : imageHeight) : 0;
            short screenMaxY = (maxY > 0) ? ((maxY < imageHeight) ? maxY : imageHeight) : 0;

            need_cal = max(screenMaxX, x2) - min(screenMinX, x1) < x2 - x1 + screenMaxX - screenMinX &&
                       max(screenMaxY, y2) - min(screenMinY, y1) < y2 - y1 + screenMaxY - screenMinY;
        }

        if (threadId < smem_size) {
            smem_sums[threadId] = need_cal;
            smem_needs[threadId] = need_cal;
            smem_idxs[threadId] = 0;
        }

        __syncthreads();
        size = nextPow2(size);

        // 通过 Exclusive Prefix Sum 求的相交圆的序号
        sharedMemExclusiveScan(threadId, smem_needs, smem_sums, smem_idxs, size);

        if (smem_needs[threadId] == 1 && threadId < size) {
            smem_idxs[smem_sums[threadId]] = threadId;
        }
        __syncthreads();

        size = smem_sums[size - 1] + smem_needs[size - 1];

        // 渲染相交圆
        for (int i = 0; i < size; ++i) {
            float3 p = ((float3*)smem_ps)[smem_idxs[i]];
            float rad = smem_rads[smem_idxs[i]];

            short minX = static_cast<short>(imageWidth * (p.x - rad));
            short maxX = static_cast<short>(imageWidth * (p.x + rad)) + 1;
            short minY = static_cast<short>(imageHeight * (p.y - rad));
            short maxY = static_cast<short>(imageHeight * (p.y + rad)) + 1;

            short screenMinX = (minX > 0) ? ((minX < imageWidth) ? minX : imageWidth) : 0;
            short screenMaxX = (maxX > 0) ? ((maxX < imageWidth) ? maxX : imageWidth) : 0;
            short screenMinY = (minY > 0) ? ((minY < imageHeight) ? minY : imageHeight) : 0;
            short screenMaxY = (maxY > 0) ? ((maxY < imageHeight) ? maxY : imageHeight) : 0;

            float invWidth = 1.f / imageWidth;
            float invHeight = 1.f / imageHeight;

            if (globalIdX < screenMinX || globalIdX >= screenMaxX) continue;
            if (globalIdY < screenMinY || globalIdY >= screenMaxY) continue;
            float2 pixelCenterNorm = make_float2(invWidth * (static_cast<float>(globalIdX) + 0.5f),
                                                 invHeight * (static_cast<float>(globalIdY) + 0.5f));
            shadePixel(smem_idxs[i], pixelCenterNorm, p, &img);
        }

	    }


    *((float4*)(cuConstRendererParams.imageData) + globalIdY * imageWidth + globalIdX) = img;
}

上述代码的性能得到了极大提高, 但是距离要求还有一定的差距. 看了知乎上的一篇解析后, 发现还有三处可以改进性能:

判断圆与图像相交可以使用文件 circleBoxTest.cu_inl 中的函数 circleInBoxConservative.
等到 smem_idxs 满了一并处理. 当前是把 256 个圆中需要渲染的部分抓紧处理, 每次可能只能处理 10 个圆左右. 因此可以每次把相交圆的 position, rad 和 color 存到 shared 中知道缓冲区满了再处理.
把 position 和 rad 合并为一个 float4.

再次修改就有了下面的代码.

// shared memory
const int smem_size = SCAN_BLOCK_DIM;
__shared__ float4 smem_p_rad[smem_size];
__shared__ float3 smem_colors[smem_size];
__shared__ uint smem_scan_scratch[smem_size * 2];
__shared__ uint smem_needs[smem_size];
__shared__ uint smem_sums[smem_size];

__global__ void kernelRenderCircles() {
    float4 img;

    const int threadId = threadIdx.y * blockDim.x + threadIdx.x;

    const unsigned int globalIdX = blockIdx.x * blockDim.x + threadIdx.x;
    const unsigned int globalIdY = blockIdx.y * blockDim.y + threadIdx.y;

    const unsigned int x1 = blockIdx.x * blockDim.x;
    const unsigned int x2 = x1 + blockDim.x;
    const unsigned int y1 = blockIdx.y * blockDim.y;
    const unsigned int y2 = y1 + blockDim.y;

    short imageWidth = cuConstRendererParams.imageWidth;
    short imageHeight = cuConstRendererParams.imageHeight;
    float invWidth = 1.f / imageWidth;
    float invHeight = 1.f / imageHeight;

    int offset = 0;

	// 读取 imageData 中的数据到寄存器
    img = __ldcs((float4*)(cuConstRendererParams.imageData) + globalIdY * imageWidth + globalIdX);

    for (int index = 0; index < cuConstRendererParams.numCircles; index += smem_size) {
        int size = min(smem_size, cuConstRendererParams.numCircles - index);
        int need_cal = 0;

        float3 p = *((float3*)cuConstRendererParams.position + index + threadId);
        float3 color = *((float3*)cuConstRendererParams.color + index + threadId);
        float rad = *(cuConstRendererParams.radius + index + threadId);

		// 判断圆与图像块是否相交
        if (threadId < size) {
            need_cal = circleInBoxConservative(p.x, p.y, rad, x1 * invWidth, x2 * invWidth, y2 * invHeight, y1 * invHeight);
        }

        if (threadId < smem_size) {
            smem_sums[threadId] = need_cal;
            smem_needs[threadId] = need_cal;
        }

        __syncthreads();
        size = nextPow2(size);

        sharedMemExclusiveScan(threadId, smem_needs, smem_sums, smem_scan_scratch, size);
        __syncthreads();
        size = smem_sums[size - 1] + smem_needs[size - 1];
		// 缓冲区满了就处理缓冲区的数据
        if (offset + size > smem_size) {
            for (int i = 0; i < offset; ++i) {
                float4 p_rad = smem_p_rad[i];
    
                float2 pixelCenterNorm = make_float2(invWidth * (static_cast<float>(globalIdX) + 0.5f),
                                                     invHeight * (static_cast<float>(globalIdY) + 0.5f));
                shadePixel(i, pixelCenterNorm, p_rad, &img);
            }

            offset = 0;
        }

		// 缓冲区没满就把 position, rad 和 color 存到缓冲区
        __syncthreads();
        if (smem_needs[threadId] == 1 && threadId < smem_size) {
            smem_p_rad[smem_sums[threadId] + offset] = make_float4(p.x, p.y, p.z, rad);
            smem_colors[smem_sums[threadId] + offset] = color;
        }
        offset += size;
    }

	// 处理缓冲区剩下的内容
    if (offset) {
        __syncthreads();
        for (int i = 0; i < offset; ++i) {
            float4 p_rad = smem_p_rad[i];

            float2 pixelCenterNorm = make_float2(invWidth * (static_cast<float>(globalIdX) + 0.5f),
                                                 invHeight * (static_cast<float>(globalIdY) + 0.5f));
            shadePixel(i, pixelCenterNorm, p_rad, &img);
        }
    }

    __stwt((float4*)(cuConstRendererParams.imageData) + globalIdY * imageWidth + globalIdX, img);
}

这样程序的性能勉强达到要求. 经过测试, 有一定几率拿到满分.

# 测试程序
make check

# 测试结果
-------------------------------------------------------------------------
| Scene Name      | Fast Time (Tf)  | Your Time (T)   | Score           |
-------------------------------------------------------------------------
| rgb             | 0.0769          | 0.0920          | 12              |
| rand10k         | 1.1239          | 1.0955          | 12              |
| rand100k        | 10.7049         | 10.3191         | 12              |
| pattern         | 0.1531          | 0.1693          | 12              |
| snowsingle      | 7.1459          | 8.4992          | 11              |
| biglittle       | 6.7287          | 3.8204          | 12              |
-------------------------------------------------------------------------
|                                   | Total score:    | 71/72           |
-------------------------------------------------------------------------