并行算法的实现与性能分析

### 并行算法的实现与性能分析 #### 1. 并行性的重要性 从程序员的角度来看，如果当今的计算机硬件是 100 GHz 的单核 CPU，而不是 3 GHz 的多核 CPU，那会非常方便，我们就无需关注并行性了。然而，计算机硬件正朝着多核 CPU 的方向发展，为了充分利用硬件性能，程序员必须使用高效的并行模式。 #### 2. 并行算法 并行编程是指利用多核硬件的编程方式。若硬件无法提供并行的优势，对算法进行并行化就毫无意义。与顺序算法相比，并行算法在算法层面上速度较慢，但其优势在于能够将算法分布到多个处理单元上执行。 衡量算法并行扩展性的简单方法是： - A：算法在单个 CPU 核心上顺序执行所需的时间。 - B：算法并行执行所需的时间乘以核心数。 若 A 和 B 相等，则算法的并行化效果完美；B 相对于 A 越大，算法的并行化效果越差。算法的并行化效果取决于每个元素的处理独立性。例如，`std::transform()` 很容易并行化，因为每个元素的处理与其他元素完全独立。理论上，对于 n 个核心，其执行速度将是顺序执行的 n 倍。但实际上，创建线程、上下文切换等诸多参数会限制并行执行。 由于并行算法的计算成本通常高于顺序算法，在某些情况下，即使顺序算法速度较慢，我们也可能会选择它。例如，当我们优化的目标是低能耗而非低计算时间时。 #### 3. 实现并行 `std::transform()` 虽然从算法角度看，`std::transform()` 易于实现，但实际上，实现一个基本的并行版本比乍看之下要复杂得多。一个简单的并行实现思路如下： 1. 将元素划分为与计算机核心数对应的块。 2. 并行地在单独的任务中执行每个块。 3. 等待所有任务完成。 以下是一个简单的实现代码： ```cpp template <typename SrcIt, typename DstIt, typename Func> auto par_transform_naive(SrcIt first, SrcIt last, DstIt dst, Func f) { auto n = static_cast<size_t>(std::distance(first, last)); auto num_tasks = std::max(std::thread::hardware_concurrency(), 1); auto chunk_sz = std::max(n / num_tasks, 1); auto futures = std::vector<std::future<void>>{}; futures.reserve(num_tasks); for (size_t task_idx = 0; task_idx < num_tasks; ++task_idx) { auto start_idx = chunk_sz * task_idx; auto stop_idx = std::min(chunk_sz * (task_idx + 1), n); auto fut = std::async([first, dst, start_idx, stop_idx, &f](){ std::transform(first+start_idx, first+stop_idx, dst+start_idx, f); }); futures.emplace_back(std::move(fut)); } for (auto& fut : futures) { fut.wait();} } ``` #### 4. 性能评估 我们通过两种场景来评估简单实现的性能： 1. 使用昂贵的函数 `heavy_f` 处理 32 个元素。 2. 使用廉价的函数 `light_f` 处理 1 亿个元素。 ```cpp // Low number of elements - heavy transform function auto heavy_f = [](float v) { auto sum = v; for (size_t i = 0; i < 100'000'000; ++i) { sum += (i*i*i*sum); } return sum; }; auto measure_heavy() { auto n = 32; auto src = std::vector<float>(n); auto dst = std::vector<float>(n); std::transform(src.begin(), src.end(), dst.begin(), heavy_f); par_transform_naive(src.begin(), src.end(), dst.begin(), heavy_f); } // High number of elements - light transform function auto light_f = [](float v) { auto sum = v; for (size_t i = 0; i < 10; ++i) { sum += (i*i*i*sum); } return sum; }; auto measure_light() { auto n = 100'000'000; auto src = std::vector<float>(n); auto dst = std::vector<float>(n); std::transform(src.begin(), src.end(), dst.begin(), light_f); par_transform_naive(src.begin(), src.end(), dst.begin(), light_f); } ``` 性能结果如下表所示： | 算法 | 时间（越低越好） | 加速比（越高越好） | | --- | --- | --- | | `std::transform()`（处理 32 个元素，`heavy_f`） | 2913914 微秒 | 1.00 x | | `par_transform_naive()`（处理 32 个元素，`heavy_f`） | 364596 微秒 | 7.99 x | | `std::transform()`（处理 1 亿个元素，`light_f`） | 407859 微秒 | 1.00 x | | `par_transform_naive()`（处理 1 亿个元素，`light_f`） | 88887 微秒 | 4.60 x | 在八核心 CPU 上，处理 32 个元素时并行化效果近乎完美，速度提升约 8 倍。而处理 1 亿个元素时，并行版本速度约为顺序版本的 5 倍，这是因为并行化会受到内存带宽等诸多外部参数的影响。 #### 5. 简单实现的缺点 简单实现仅在我们是唯一使用硬件的应用程序，且每个块的计算成本相同时效果较好。但实际情况并非如此，我们需要一个通用的并行实现。当每个块的计算成本不相等时，实现会受到耗时最长的块的限制；当应用程序或操作系统有其他进程需要处理时，操作无法并行处理所有块。将操作拆分为更小的块可以使并行化适应当前条件，避免单个任务拖慢整个操作。 #### 6. 分治法 为了实现通用的并行转换，我们采用分治法，将范围递归地划分为更小的范围。具体步骤如下： 1. 将输入范围划分为两个范围。若输入范围小于指定阈值，则处理该范围；否则，将范围拆分为两部分： - 一部分分支到另一个任务中递归处理。 - 一部分在调用线程中递归处理。 以下是分治法的实现代码： ```cpp template <typename SrcIt, typename DstIt, typename Func> auto par_transform(SrcIt first,SrcIt last,DstIt dst,Func f,size_t chunk_sz) { const auto n = static_cast<size_t>(std::distance(first, last)); if (n <= chunk_sz) { std::transform(first, last, dst, f); return; } const auto src_middle = std::next(first, n/2); auto future = std::async([=, &func]{ par_transform(first, src_middle, dst, f, chunk_sz); }); const auto dst_middle = std::next(dst, n/2); par_transform(src_middle, last, dst_middle, f, chunk_sz); future.wait(); } ``` #### 7. 性能评估 我们创建一个根据输入值耗时不同的 `transform_func`，并使用 `std::iota()` 生成递增的值范围。通过不同的块大小进行评估，结果如下表所示： | 函数 | 块大小 | 任务数 | 微秒 | 加速比 | | --- | --- | --- | --- | --- | | `std::transform()` | 10000000 (= n) | 1 | 844629 | 1.00x | | `par_transform_naive()` | 1250000 (= n / 8) | 8 | 222933 | 3.79x | | `par_transform()` | 1000000 | 10 | 210942 | 4.00x | | `par_transform()` | 100000 | 100 | 148066 | 5.70x | | `par_transform()` | 10000 | 1000 | 144189 | 5.86x | | `par_transform()` | 1000 | 10000 | 152123 | 5.55x | | `par_transform()` | 100 | 100000 | 208969 | 4.04x | | `par_transform()` | 10 | 1000000 | 1536680 | 0.55x | 从表中可以看出，在这种情况下，块大小约为 10000 个元素时性能最佳。块太大，性能会受到处理最终块所需时间的限制；块太小，创建和调用任务的开销会比计算本身更大。 若将 `transform_func` 改为固定计算成本的函数，再次执行代码，结果如下表所示： | 函数 | 块大小 | 任务数 | 微秒 | | --- | --- | --- | --- | | `std::transform()` | 10000000 (= n) | 1 | 815498 | | `par_transform_naive()` | 1250000 (= n / 8) | 8 | 129403 | | `par_transform_chunks()` | 1000000 | 10 | 184041 | | `par_transform_chunks()` | 100000 | 100 | 132248 | | `par_transform_chunks()` | 10000 | 1000 | 131812 | | `par_transform_chunks()` | 1000 | 10000 | 141705 | | `par_transform_chunks()` | 100 | 100000 | 179279 | | `par_transform_chunks()` | 10 | 1000000 | 1542512 | 最快的版本是 `par_transform_naive()`，但差距不大。这表明，通用实现使用大量任务比根据 CPU 核心数确定任务数更为明智，应让异步任务的调度器进行调度。 #### 8. 实现并行 `std::count_if` 我们可以使用相同的分治法来实现并行版本的 `std::count_if()`，不同之处在于需要累加返回值： ```cpp template <typename It, typename Pred> auto par_count_if(It first, It last, Pred pred, size_t chunk_sz) { auto n = static_cast<size_t>(std::distance(first, last)); if (n <= chunk_sz) return std::count_if(first, last, pred); auto middle = std::next(first, n/2); auto future = std::async([=, &pred]{ return par_count_if(first, middle, pred, chunk_ ```