十亿行数据挑战：CUDA申请出战（从17分钟到17秒）创作不易，方便的话点点关注，谢谢文章结尾有最新热度的文章，感兴趣

创作不易，方便的话点点关注，谢谢

文章结尾有最新热度的文章，感兴趣的可以去看看。

本文是经过严格查阅相关权威文献和资料，形成的专业的可靠的内容。全文数据都有据可依，可回溯。特别申明：数据和资料已获得授权。本文内容，不涉及任何偏颇观点，用中立态度客观事实描述事情本身

文章有点长(4500字阅读时长：10分)，期望您能坚持看完，并有所收获。

导读

在我学习CUDA的过程中，我决定用它来应对十亿行数据挑战。这个挑战本身很简单，但用CUDA来实现并没有这么容易。在这里，我将分享我的解决方案，它在V100显卡上运行耗时16.8秒。当然，这绝不是最快的解决方案，但它是首个此类（不使用cudf，仅手写内核）方案。我向其他CUDA爱好者发起挑战，看谁能让它运行得更快。

纯C++基准实现

那我们就先用纯C++代码做个基准实现吧。我的CUDA代码应该比这个更快才对。

实现方法很直接：逐行读取文件，解析出城市名和温度值，然后将它们累积到一个标准模板库（STL）的map中。

while (getline(file, line)){
istringstream iss(line);
    string station;
float temp;
getline(iss, station,';');
    iss >> temp;

auto it = stationStats.find(station);
if(it == stationStats.end()){
        stationStats[station]={temp, temp, temp,1};
}else{
Stat& s = it->second;
        s.min =min(s.min, temp);
        s.max =max(s.max, temp);
        s.sum += temp;
        s.count++;
}
}

ofstream measurements("measurements.out");
for(auto& pair : stationStats){
constStat& s = pair.second;
float mean = s.sum / s.count;
    measurements << pair.first <<"="<< s.min <<"/";
    measurements << fixed <<setprecision(1)<< mean <<"/";
    measurements << s.max << endl;
}

这段代码运行需要16.5分钟。让我们用CUDA来改进它吧。

工作划分方法

十亿个线程？CUDA以及其他并行编程应用程序接口（API）的优势在于，你可以在多个进程间并行化工作负载。对于CUDA来说，它采用的是单指令多线程（SIMT）模型——单个指令可以并行地在多个线程中执行。

很好，那我们就使用十亿个线程来同时处理十亿行数据吧！

遗憾的是，我们不能直接启动十亿个线程。我们首先需要为每个线程准备好要处理的每行数据缓冲区。然而，准备这十亿个行缓冲区需要逐行读取整个文件（除非这些行已经被放在十亿个文件中了，但那样的话这就变成十亿个文件挑战了）。

设置这些缓冲区所涉及的工作量本质上会复制基准实现的工作量，使得这种方法适得其反。

使用字节偏移量

解决方案是准备文件偏移量，而不是行缓冲区。这些偏移量是通过迭代获取的，按照期望的划分大小（=文件总大小 / 期望划分的份数）逐步遍历整个文件缓冲区，并标记换行符的位置：

long long split_size = size / num_parts;
longlong offset =0;
std::vector<Part> parts;
while(offset < size){
longlong seek_offset = std::max(offset + split_size - MAX_CITY_BYTE,0LL);
if(seek_offset > size){
        parts.back().length += size - offset;
break;
}
    file.seekg(seek_offset, std::ios::beg);
char buf[MAX_CITY_BYTE];
    file.read(buf, MAX_CITY_BYTE);

// 在附近查找换行符。
// 它将是当前偏移量和下一个偏移量之间的边界。
    std::streamsize n = file.gcount();
    std::streamsize newline =-1;
for(int i = n -1; i >=0;--i){
if(buf[i]=='\n'){
            newline = i;
break;
}
}
int remaining = n - newline -1;
longlong next_offset = seek_offset + n - remaining;
    parts.push_back({offset, next_offset - offset});
    offset = next_offset;
}

这种方式比读取整个文件要快得多，因为我们处理的是整型字节值。例如，假设我们要将14GB的输入文件划分给两个线程。上述while循环会迭代两次（偏移量从0GB到7GB再到14GB）。对比基于行的处理方式，我们需要迭代5亿（= 10亿 / 2）次才能将5亿行数据加载到第一个分区中。

实际上，我们需要的线程不止两个。但十亿个线程又太多了。不是因为我们的GPU无法处理十亿个线程，而是因为查找十亿个偏移量会成为整个运行时间的瓶颈。我们甚至还没开始启动CUDA内核呢。我们需要尽可能减少准备时间。

在我的解决方案中，我创建了100万个分区，这在整个16.8秒的运行时间里占用了2.6秒。（相比之下，仅创建1亿个分区就需要超过3分钟。）

CUDA内核

其余的时间都花在了CUDA内核上（终于到这一步了！）。其背后的思路很简单。每个线程索引到文件缓冲区的不同部分，解析出城市名和温度值，并更新最小值、最大值和平均值统计信息。

然而，由于以下原因（按棘手程度排序），实际实现起来并不简单：

1. CUDA的AtomicMin（原子最小值操作）和AtomicMax（原子最大值操作）只适用于整型值。我们需要自己编写能处理浮点值的变体。
2. CUDA中没有std::string。是时候编写我们自己的atof（字符串转浮点数函数）、strcmp（字符串比较函数）、getline（读取一行函数）了。准备好处理空终止符'\0'吧。
3. 也没有std::map。我们该如何将城市字符串传递给数组索引查找表并传入CUDA内核呢？

让我们逐个来看这些问题。

用于浮点数的`AtomicMin`和`AtomicMax`

CUDA中的任何原子操作都可以用atomicCAS（原子比较并交换操作）来编写。让我们参照官方编程指南中的示例，为浮点值编写AtomicMin和AtomicMax函数。以下是AtomicMin函数的代码：

__device__ static float atomicMin(float* address, float val) {
int* address_as_i =(int*) address;
int old =*address_as_i, assumed;
do{
        assumed = old;
        old =::atomicCAS(address_as_i, assumed,
// 使用`fmaxf`来实现原子最大值操作。
        __float_as_int(::fminf(val, __int_as_float(assumed))));
}while(assumed!= old);
return __int_as_float(old);
}

现在我们就可以原子性地更新温度的最小、最大值（浮点型）了。

C字符串

虽然我们没有std::string，但我们有char*。字符串本质上就是一个8位字符组成的数组，原始文件缓冲区（例如“Hamburg;12.0\nBulawayo;8.9\nPalembang;38.8...”）也没什么不同。

每个线程遍历这个字符数组，偏移量和长度各不相同（在我们进行工作划分步骤时已计算好）：

char city[MAX_CITY_BYTE];
char floatstr[5];// 最长的温度浮点数字符串是 -99.9，即5个字节

for(int i =0; i < parts[bx].length; i++){// bx是全局线程索引
char c = buffer[parts[bx].offset - buffer_offset + i];
if(parsing_city){// 城市字符
if(c ==';'){
            city[index]='\0';
            index =0;
            parsing_city =false;
}else{
            city[index]= c;
            index++;
}
}else{// 浮点数字符
if(c =='\n'){// 到达行尾
            floatstr[index]='\0';

int stat_index =get_index(cities, city, n_city);
float temp =cuda_atof(floatstr);

// CUDA内核的核心部分。
// （原子性地）更新温度统计信息。
// 本质上和简单的C++版本思路相同。
atomicMin(&stats[stat_index].min, temp);
atomicMax(&stats[stat_index].max, temp);
atomicAdd(&stats[stat_index].sum, temp);
atomicAdd(&stats[stat_index].count,1);

// 为下一次读取行做重置
            parsing_city =true;
            index =0;
            floatstr[0]='\0';
            city[0]='\0';
}else{
            floatstr[index]= c;
            index++;
}
}
}

这样写不太美观，但由于我们没有getline这样的便利函数，所以这是必要的。解析完每一行后，我们现在就得到了一对字符串：一个城市名字符串char city[]和一个温度字符串char floatstr[]。（后者需要从字符串转换为浮点数，由于CUDA中没有atof函数，我又自己编写了一个）

城市字符串到索引查找

GPU哈希表？我们如何为每个城市存储温度统计信息呢？在C++中，我们依靠哈希表——以城市字符串作为键，温度统计信息（浮点数）作为值。但在CUDA中，我们没有这么方便的std::map。

好吧，那我们自己写一个。能有多难呢？结果发现，难到几乎不可能，因为我的城市字符串作为键有100个字节长。

虽然我在网上找到了一些实现方式，但这些方法由于原子操作受限于有限的位数（即使在CPU上也是如此），所以它们仅限于32位的键。

需要明确的是，处理哈希表冲突并不一定需要原子操作，但当冲突可能在多个线程间发生时，也就是在并行环境下有并发插入哈希表的情况时，就需要原子操作了。

先斩后奏（变通做法）

所以在这一点上，我稍微变通了一下原始挑战的规则。我假设输入文件会附带一份所有城市的列表。这个列表就是data/weather_stations.csv文件，它实际上是用于生成官方挑战中十亿行数据的文件。

变通解决方案：排序后的城市列表 + 二分查找知道了所有可能的城市列表后，我就避免使用哈希表了。我只需对所有城市的列表进行排序，然后将其传递给CUDA内核。我使用排序后的索引作为查找表。例如，假设排序后的城市列表是["A","B","C"]。给定城市“B”，查找索引就是它的位置，即1。

排序很重要，因为我们可以通过二分查找在对数时间内找到索引。虽然这比哈希表的常数时间查找要慢，但比线性搜索4万多个城市条目要快得多。

// 从没想过我会在CUDA中写二分查找，但现在就是写了。
// 感谢力扣（LeetCode）！
__device__ int get_index(char* cities, char* city_target, int n_city) {
int left =0;
int right = n_city -1;
while(left <= right){
int mid = left +(right - left)/2;
constchar* city_query = cities + mid * MAX_CITY_BYTE;

int cmp =cuda_strcmp(city_query, city_target);
if(cmp ==0)
return mid;
elseif(cmp <0)
            left = mid +1;
else
            right = mid -1;
}
return-1;
}

现在我们终于完成了！内核的核心部分归结为以下代码：

int stat_index = get_index(cities, city, n_city);
float temp = cuda_atof(floatstr);

atomicMin(&stats[stat_index].min, temp);
atomicMax(&stats[stat_index].max, temp);
atomicAdd(&stats[stat_index].sum, temp);
atomicAdd(&stats[stat_index].count, 1);

做这一切仅仅是为了执行这4个原子操作。

性能分析

在V100显卡上，CUDA解决方案运行耗时16.8秒。与我们的C++基准实现的16.5分钟相比，性能提升了60倍。以下是重现此结果的脚本。

有趣的是，这个内核在T4显卡上运行大约要慢1.5倍，所以我使用ncu工具在这两种设备上进行了性能分析。有一个引起我注意的地方是控制分歧（control divergence）方面的差异。

（V100）

~ ncu --section SourceCounters fast data/measurements.txt 1000000600000
...
Section:SourceCounters
--------------------------------------------------
MetricNameMetricUnitMetricValue
--------------------------------------------------
BranchInstructionsRatio%0.18
BranchInstructions              inst 53,736,337,800
BranchEfficiency%92.73
Avg.DivergentBranches10,041,990.22
------------------------------------ --------------

（T4）

~ ncu --section SourceCounters fast data/measurements.txt 1000000600000
...
Section:SourceCounters
--------------------------------------------------
MetricNameMetricUnitMetricValue
--------------------------------------------------
BranchInstructionsRatio%0.17
BranchInstructions              inst 53,870,435,921
BranchEfficiency%92.73
Avg.DivergentBranches20,156,806.42
------------------------------------ --------------

我原本预计会有很多控制分歧（因为我在内核中到处都使用了if语句和循环），但没想到在T4上分歧会更严重。与V100相比，它的平均分歧分支数是V100的两倍。这会不会是内核在T4上运行更慢的原因呢？而且，如果有这么多分歧分支，为什么分支效率还这么高呢？

显然，我还在学习ncu性能分析器的使用。欢迎大家在这方面给予我指导。

可能的优化——使用共享内存私有化

在我写完这篇博客时，我意识到我的Stat结构体不需要保存城市字符数组。这样的话，每个结构体将是16字节（最小值、最大值、总和、计数这几个成员），整个统计信息数组将是16N字节，其中N是唯一城市的数量。

在这里，N = 41000（基于data/weather_stations.csv文件），所以整个数组将是66KB。这个大小应该足够小，可以放入共享内存中（对于Volta架构，每个流式多处理器（SM）有96KB的共享内存）。然后，每个线程块可以更新统计信息的私有版本，以减少原子操作的争用。这应该能得到一个更快的解决方案。