TechThings

Posted on Edited on

本篇介绍 cuda samples 中的 mergeSort.

大体上来讲, mergeSort 分为两个阶段.

  1. 对含有 SHARED_SIZE_LIMIT (即 1024) 个元素的数组进行排序.
  2. 合并多个有序数组.

其中第一个阶段调用一次函数 mergeSortShared 结束. 而第二个阶段需要循环调用三个函数: generateSampleRanks, mergeRanksAndIndices 和 mergeElementaryIntervals. 我将分别讲述这两个阶段.

第一阶段: 段内排序

函数 mergeSortShared 比较简单, 就是对完整的数据执行 mergeSortSharedKernel 函数. 其函数声明为

1
2
3
4
template <uint sortDir>
__global__ void mergeSortSharedKernel(uint *d_DstKey, uint *d_DstVal,
                                      uint *d_SrcKey, uint *d_SrcVal,
                                      uint arrayLength);

该 kernel 函数对一个长为 arrayLength 的数组排序, 在代码中 arrayLength 接收的参数为宏变量 SHARED_SIZE_LIMIT (即 1024). 调用该函数的 grid 和 block 大小分别为 N / arrayLength 和 arrayLength / 2. 每个 block 排序一段长为 arrayLength 的数组.

在函数 mergeSortSharedKernel 中排序通过一个循环进行. 该循环在长度为 stride 的数组已经排序基础上, 将 2 个长为 stride 的数组进行合并. 因此 stride 的大小从 1 开始每次增大 1 倍直到 arrayLength / 2. 在循环开始之前, 长度 stride 的数组已经是排序了的状态.

循环中数量为 stride 的 thread 为一组, 排序相邻的两个长为 stride 的数组, 记前一个数组为 A, 后一个数组为 B. 因此需要计算每一组 thread 内部的 thread ID 以及这组 thread 对应数据的起始位置.

1
2
3
4
5
6
// 计算组内 thread ID, [0, thread)
uint lPos = thread & (stride - 1);

// 计算该组 thread 对应数据的起始位置
// 每组 thread 处理连续的两个 thread 数组
uint *baseKey = s_key + 2 * (threadIdx.x - lPos);

相对应的, 每个 thread 在两个数组中都有一个对应的数据, 分别为 keyA 和 keyB.

1
2
uint keyA = baseKey[lPos + 0];
uint keyB = baseKey[lPos + stride];

接下来找到 keyA 和 keyB 在两个数组合并后的位置. 因为数组 A 和数组 B 都是排序的, 因此在数组 A 中有 lPos 个数比 keyA 小, 只需要找到在数组 B 中有多少个数比 keyA 小, 两者相加就得到了 keyA 合并后的位置. 这个寻找的过程可以通过二分法查找. keyB 同理也能找到合并后的位置. 但是这里有个问题就是某个 thread 对应的 keyA 和 keyB 相同, 那么这样计算得到的最后的位置也是相同的, 产生了冲突.

代码中解决这个可能的冲突的方法是 keyA 在数组 B 中寻找小于 keyA 的数, keyB 在数组 A 中寻找小于等于 keyB 的数. 这样做的好处是保持了排序的稳定性. 数组 A 的元素总是排在数组 B 的相同元素之前.

1
2
uint posA = binarySearchExclusive<sortDir>(keyA, baseKey + stride, stride, stride) + lPos;
uint posB = binarySearchInclusive<sortDir>(KeyB, baseKey + 0, stride, stride) + lPos;

函数 binarySearchExclusive 和 binarySearchInclusive 有着相似的结构, 只有内部的微弱区别. 通过二分法查找 val 在数组中的位置, 即返回数组中小于 (binarySearchExclusive) 或小于等于 (binarySearchInclusive) val 的元素的数量.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
// val 待查找的元素
// data 数组起始位置
// L 数组长度
// stride >= L, 并行是 2 的整数幂
template <uint sortDir>
static inline __device__ uint binarySearchExclusive(uint val, uint *data,
                                                    uint L, uint stride) {
  if (L == 0) {
    return 0;
  }

  uint pos = 0;

  for (; stride > 0; stride >>= 1) {
    uint newPos = umin(pos + stride, L);

	// binarySearchInclusive 这里的 < 变成 <=
    if ((sortDir && (data[newPos - 1] < val)) ||
        (!sortDir && (data[newPos - 1] > val))) {
      pos = newPos;
    }
  }

  return pos;
}

merge sort 的第一个阶段就完成了. 该阶段将每 SHARED_SIZE_LIMIT (即 1024) 个元素进行排序, 得到了 N / SHARED_SIZE_LIMIT 个有序数组.

第二阶段: 合并有序段

第二个阶段循环合并连续 2 个长为 stride 的有序数组. 因此 stride 的大小从 SHARED_SIZE_LIMIT (即 1024) 开始每次增大 1 倍直到 N. 这个阶段主要分成三个部分: 生成排序, 合并排序和合并. 也就是调用三个函数 generateSampleRanks, mergeRanksAndIndices 和 mergeElementaryIntervals.

思路解析

要合并 2 个连续的有序数组, 当然可以像第一阶段一样启动足够多的线程查找每个元素合并后的位置. 但是这样的方法在后期的时延会增加很多, 因为 stride 长度倍增. 这里采用的方法是先分组再查找. 要理解这部分的代码, 需要理解数组 d_RanksA, d_RanksB, d_LimitsA 和 d_LimitsB. 这个会在后面详细讲解.

先把 stride 分为大小为 SAMPLE_STRIDE 的子数组. d_RanksA 表示 2 个 stride 中所有的子数组的起始元素在 stride A 中的位置, 即 A 中小于 (或小于等于) 某个起始元素的元素数. 这一部分对应的就是函数 generateSampleRanks 的部分.

接着把 d_RanksA 排序, 将排序后的值填入对应的 d_LimitsA 中. 这一部分对应的就是函数 mergeRanksAndIndices 的部分.

经过前两步, 就把 A 和 B 每个都分为了 2×stride/SAMPLE_STRIDE2 \times stride / SAMPLE\_STRIDE 个范围. 之后合并对应的范围就得到 2 个 stride 合并后的结果. 这一部分对应的就是函数 mergeElementaryIntervals.

代码解析

generateSamplesRanks

这一步生成 d_RanksA 和 d_RanksB. 通过 kernel 函数 generateSampleRanksKernel 实现的. 无论 stride 大小, 该 kernel 函数的总线程数都不变. 就是 N/(2×SAMPLE_STRIDE)N / (2 \times SAMPLE\_STRIDE).

1
2
3
4
5
6
// 计算总线程数量
uint lastSegmentElements = N % (2 * stride);
uint threadCount =
	(lastSegmentElements > stride)
	    ? (N + 2 * stride - lastSegmentElements) / (2 * SAMPLE_STRIDE)
        : (N - lastSegmentElements) / (2 * SAMPLE_STRIDE);

调用函数 generateSampleRanksKernel 的 grid 大小为 threadCount / 256, block 大小为 256. 在该 kernel 函数中, 把所有的线程分组. 每 stride / SAMPLE_STRIDE 个线程组成一组. 因此, 每个线程的全局 ID 和组内 ID 可以通过如下方式计算.

1
2
3
4
5
6
// 线程的全局 ID
uint pos = blockIdx.x * blockDim.x + threadIdx.x;

// stride / SAMPLE_STRIDE 个线程为一组
// 线程的组内 ID
const uint i = pos & ((stride / SAMPLE_STRIDE) - 1);

每个 thread 对应两个长为 SAMPLE_STRIDE 的数组以及 d_RanksA 和 d_RanksB 中的两个元素, 就可以计算每个线程组的数据偏移量.

1
2
3
4
5
// 计算每个线程组的数据偏移量
const uint segmentBase = (pos - i) * (2 * SAMPLE_STRIDE);
d_SrcKey += segmentBase;
d_RanksA += segmentBase / SAMPLE_STRIDE;
d_RanksB += segmentBase / SAMPLE_STRIDE;

每个线程组对应连续 2 个长为 stride 的数组, 前一个数组称为数组 A, 后一个称为 B. 每个线程在 A 和 B 中都有对应的一个长为 SAMPLE_STRIDE 的子数组. 就可以计算得到 A 和 B 的长度以及子数组的个数.

1
2
3
4
5
6
7
// A 和 B 的长度
const uint segmentElementsA = stride;
const uint segmentElementsB = umin(stride, N - segmentBase - stride);

// 子数组的个数
const uint segmentSamplesA = getSampleCount(segmentElementsA);
const uint segmentSamplesB = getSampleCount(segmentElementsB);

最后求在 A 和 B 两个数组中, 小于 (或者小于等于) 每个子数组起始元素的元素数. 这也是通过 binarySearchExclusive 和 binarySearchInclusive 两个 kernel 函数完成.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
if (i < segmentSamplesA) {
	d_RanksA[i] = i * SAMPLE_STRIDE;
    d_RanksB[i] = binarySearchExclusive<sortDir>(
        d_SrcKey[i * SAMPLE_STRIDE], d_SrcKey + stride, segmentElementsB,
        nextPowerOfTwo(segmentElementsB));
}

if (i < segmentSamplesB) {
    d_RanksB[(stride / SAMPLE_STRIDE) + i] = i * SAMPLE_STRIDE;
    d_RanksA[(stride / SAMPLE_STRIDE) + i] = binarySearchInclusive<sortDir>(
        d_SrcKey[stride + i * SAMPLE_STRIDE], d_SrcKey + 0, segmentElementsA,
        nextPowerOfTwo(segmentElementsA));
}

mergeRanksAndIndices

这一步把 d_RanksA 和 d_RanksB 中的内容排序后赋值给 d_LimitsA 和 d_LimitsB.

对 A 和 B 分别调用 kernel 函数 mergeRanksAndIndicesKernel, 通过 d_RanksA 生成 d_LimitsA. 该 kernel 函数有大小 <<<threadCount / 256, 256>>>. 其中 threadCount = N / (2 * SAMPLE_STRIDE).

同样 stride / SAMPLE_STRIDE 个 thread 为一组. 计算组内 ID 和数组偏移同之前一样.

1
2
3
4
5
6
7
8
9
10
// thread ID
uint pos = blockIdx.x * blockDim.x + threadIdx.x;

// 组内 ID
const uint i = pos & ((stride / SAMPLE_STRIDE) - 1);

// 数组偏移
const uint segmentBase = (pos - i) * (2 * SAMPLE_STRIDE);
d_Ranks += (pos - i) * 2;
d_Limits += (pos - i) * 2;

需要将 d_RanksA 排序. 对 A 中的每个 d_RanksA 查找其在 B 中的位置.

1
2
3
4
5
6
7
8
if (i < segmentSamplesA) {
    uint dstPos = binarySearchExclusive<1U>(
                      d_Ranks[i], d_Ranks + segmentSamplesA, segmentSamplesB,
                      nextPowerOfTwo(segmentSamplesB)) +
                  i;
    d_Limits[dstPos] = d_Ranks[i];
}

同样, 对 B 中的每个 d_RanksA 查找其在 A 中的位置.

1
2
3
4
5
6
7
if (i < segmentSamplesA) {
    uint dstPos = binarySearchExclusive<1U>(
                      d_Ranks[i], d_Ranks + segmentSamplesA, segmentSamplesB,
                      nextPowerOfTwo(segmentSamplesB)) +
                  i;
    d_Limits[dstPos] = d_Ranks[i];
}

经过排序, 实际上是把一个 stride 分为了 2 * stride / SAMPLE_STRIDE 个部分. 划分的标准是每个 SAMPLE_STRIDE 的起始元素在 A 或 B 中的位置.

mergeElementaryIntervals

这一步合并 2 个 stride. 调用 kernel 函数 mergeElementaryIntervalsKernel, 其大小为 <<<N / SAMPLE_STRIDE, SAMPLE_STRIDE>>>. 之前已经把每个 stride 分为了 2 * SAMPLE_STRIDE 个部分. 每个 block 负责 2 个 stride 中对应的一部分, 将他们合并.

首先是把所有的 block 分组, 每组有 2 * stride / SAMPLE_STRIDE 个 block. 再计算组内 ID 和数据偏移.

1
2
3
4
5
6
7
8
9
// 组内 ID
const uint intervalI = blockIdx.x & ((2 * stride) / SAMPLE_STRIDE - 1);

// 数据偏移
const uint segmentBase = (blockIdx.x - intervalI) * SAMPLE_STRIDE;
d_SrcKey += segmentBase;
d_SrcVal += segmentBase;
d_DstKey += segmentBase;
d_DstVal += segmentBase;

还需要计算合并后的起始位置, 该 block 的合并长度等数据.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
__shared__ uint startSrcA, startSrcB, lenSrcA, lenSrcB, startDstA, startDstB;

// 本 block 读取数据的起始位置
startSrcA = d_LimitsA[blockIdx.x];
startSrcB = d_LimitsB[blockIdx.x];

// 本 block 读取数据的终止位置
// 下一个 block 读取数据的开始位置或者全 stride 长度
uint endSrcA = (intervalI + 1 < segmentSamples) ? d_LimitsA[blockIdx.x + 1]
                                                    : segmentElementsA;
uint endSrcB = (intervalI + 1 < segmentSamples) ? d_LimitsB[blockIdx.x + 1]
                                                    : segmentElementsB;

// 本 block 读取的数据长度
lenSrcA = endSrcA - startSrcA;
lenSrcB = endSrcB - startSrcB;
// 本 block 合并完需要写入的起始位置
// A 的起始位置就是之前 block 的合并的长度和
// B 的起始位置就是 A 的起始位置 + 长度
startDstA = startSrcA + startSrcB;
startDstB = startDstA + lenSrcA;

把 global 内存中的数据加载到 shared 内存.

1
2
3
4
5
6
7
8
9
10
11
if (threadIdx.x < lenSrcA) {
    s_key[threadIdx.x + 0] = d_SrcKey[0 + startSrcA + threadIdx.x];
    s_val[threadIdx.x + 0] = d_SrcVal[0 + startSrcA + threadIdx.x];
}

if (threadIdx.x < lenSrcB) {
    s_key[threadIdx.x + SAMPLE_STRIDE] =
        d_SrcKey[stride + startSrcB + threadIdx.x];
    s_val[threadIdx.x + SAMPLE_STRIDE] =
        d_SrcVal[stride + startSrcB + threadIdx.x];
}

调用 merge 函数合并 s_key 中的数据. merge 函数比较简单, 找到每个元素合并后的位置.

最后把数据存到 global 内存的对应位置.

1
2
3
4
5
6
7
8
9
10
if (threadIdx.x < lenSrcA) {
    d_DstKey[startDstA + threadIdx.x] = s_key[threadIdx.x];
    d_DstVal[startDstA + threadIdx.x] = s_val[threadIdx.x];
}

if (threadIdx.x < lenSrcB) {
    d_DstKey[startDstB + threadIdx.x] = s_key[lenSrcA + threadIdx.x];
    d_DstVal[startDstB + threadIdx.x] = s_val[lenSrcA + threadIdx.x];
}

Posted on

我在照着STL 源码刨析写 STL 的时候出现了一个编译错误. 当时的情况是我要用一个派生类的指针给一个基类指针的引用赋值, 类似于:

1
2
3
4
5
class Base {};
class Derived : public Base {};

Derived *d = new Derived();
Base *&b = d;

但是编译器报错, 错误信息类似于:

1
error: cannot bind non-const lvalue reference of type ‘Base*&’ to an rvalue of type ‘Base*’

在 C++ 中, 派生类的指针可以直接给基类指针赋值以实现多态的特性. 但是却不可以给基类指针的引用赋值. 这个报错当时令我很困惑. 经过查阅并分析报错信息, 该语句的执行过程应该是首先创建了一个类型为 Base* 的临时值, 用 d 来初始化, 再把这个临时值赋予变量 b. 因为这个临时值是右值不能给左值赋值, 所有会有这个报错.

这个问题可以用两种方式解决.

第一种方式是通过常量引用可以指向右值实现的. 就是把变量 b 声明为常量指针, 即 Base * const &b = d, 这样就解决了这个问题.

但是在我当时的代码中不能通过这个方法解决. 因为报错的函数是红黑树的一个迭代函数, 有可能改变 d 的值. 这样的解决方法会复制 d 的值, 使得函数内改变的局部值无法作用到外部的 d 上.

第二种方式是先把 d 做一个强制类型转换再赋值给 b, 即 Base *&b = (Base*&)d. 这样就可以在函数中通过修改 b 的值修改 d.

下面是一个示例代码. 验证两种方法.

1
2
3
4
5
6
7
8
9
10
class Base {};
class Derived : public Base {};

Derived *pd = new Derived();
Base * const &pb1 = pd;   // 方法 1
Base *&pb2 = (Base*&)pd;  // 方法 2

cout << "addr of pd:\t" << &pd << endl;
cout << "addr of pb1:\t" << &pb1 << endl;
cout << "addr of pb2:\t" << &pb2 << endl;

编译没有报错, 运行得到下面的结果:

1
2
3
add of pd:	0x7ffcdec09978
add of pb1:	0x7ffcdec09980
add of pb2:	0x7ffcdec09978

说明方法 1 创建了一个新的指针, 而方法 2 绑定在原来的指针上.

Posted on Edited on

这个作业有三个部分, 都是 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 就行.

1
2
3
4
5
6
7
__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 实现有十几倍的提升.

1
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 上的实现分为三个部分:

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

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

首先是三个辅助函数.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
__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 函数.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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);
    }
}

运行测试程序查看得分.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 测试程序
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 加速该函数.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
__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;
}

运行测试程序查看得分.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 测试程序
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 的性能勉强达到要求.

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
// 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, 可以直接调用.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
// 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;
}

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
// 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];

__device__ __inline__ void
shadePixel(int circleIndex, float2 pixelCenter, float4 p_rad, float4* imagePtr) {

    float diffX = p_rad.x - pixelCenter.x;
    float diffY = p_rad.y - pixelCenter.y;
    float pixelDist = diffX * diffX + diffY * diffY;

    float rad = p_rad.w;
    float maxDist = rad * rad;

    if (pixelDist > maxDist)
        return;

    float3 rgb;
    float alpha;

    if (cuConstRendererParams.sceneName == SNOWFLAKES || cuConstRendererParams.sceneName == SNOWFLAKES_SINGLE_FRAME) {

        const float kCircleMaxAlpha = .5f;
        const float falloffScale = 4.f;

        float normPixelDist = sqrt(pixelDist) / rad;
        rgb = lookupColor(normPixelDist);

        float maxAlpha = .6f + .4f * (1.f-p_rad.z);
        maxAlpha = kCircleMaxAlpha * fmaxf(fminf(maxAlpha, 1.f), 0.f); // kCircleMaxAlpha * clamped value
        alpha = maxAlpha * exp(-1.f * falloffScale * normPixelDist * normPixelDist);

    } else {
        rgb = smem_colors[circleIndex];
        alpha = .5f;
    }

    float oneMinusAlpha = 1.f - alpha;

    // BEGIN SHOULD-BE-ATOMIC REGION

    float4 existingColor = *imagePtr;
    float4 newColor;
    newColor.x = alpha * rgb.x + oneMinusAlpha * existingColor.x;
    newColor.y = alpha * rgb.y + oneMinusAlpha * existingColor.y;
    newColor.z = alpha * rgb.z + oneMinusAlpha * existingColor.z;
    newColor.w = alpha + existingColor.w;

    *imagePtr = newColor;

    // END SHOULD-BE-ATOMIC REGION
}

__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);
}

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 测试程序
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           |
-------------------------------------------------------------------------
0%