STANDARD PARALLELISM

Bryce Adelstein Lelbach
HPC Programming Models Architect
Standard C++ Library Evolution Chair, US Programming Languages Chair
@blelbach

目录

• 约定 (Conventions)

  • 命名空间别名
  • 类模板参数推导 (CTAD)

• C++ 可移植性

  • 我们为何需要“加速匝道” (We Need On-Ramps)
  • 什么使得 C++ 可移植？ (What Makes C++ Portable?)
  • 新的可移植性契约 (The New Portability Contract)
  • C++ 标准是描述性的，而非规定性的。
  • 实现自由 (Implementation Freedom)

• C++ 标准行为与并行性

  • 标准并行性位于库中
  • C++ 标准并行性的支柱

• C++ 标准算法

  • 串行算法 (C++98)
  • 并行算法 (C++17)

• C++ 标准库执行策略 (Execution Policy)

• 单词计数示例（使用 std::transform_reduce 和并行策略）
• C++20 Ranges 介绍

• 并行算法与数据结构示例

  • std::for_each 和 Ranges 概念图 (std::vector)
  • std::for_each 和 Ranges 概念图 (std::iota)
  • 使用 Ranges 进行矩阵操作
  • 并行 for_each 阻塞行为
  • 变换与并行 for_each 操作
  • 过滤与归约操作

• C++ 并行应用案例

  • LULESH
  • STLBM
  • M-AIA

• 异步模型: Senders & Receivers

  • 当前 C++ 现状与解决方案
  • C++ 并发编程示例
  • 核心概念定义
  • 执行上下文与调度器
  • 发送器函数属性
  • 发送器链式操作示例
  • 发送器适配器 (Sender Adaptors)
  • 发送器工厂 (Sender Factories)
  • Sender 消费者
  • Sender/Receiver 链式操作数据流示例
  • inclusive_scan_async 适配器实现
  • 麦克斯韦方程组 (Maxwell's Equations) 求解
  • 麦克斯韦方程组性能加速
  • Palabos 碳捕集模拟

• C++ 标准算法的演进

• 多维数据抽象与 std::mdspan

  • 多维数据抽象的缺失与解决方案
  • std::mdspan 介绍
  • std::extents 概述
  • std::mdspan 模板定义及示例

  • 多维数组布局与 std::mdspan

    • std::mdspan 基础用法
    • 行主序 (Row-Major AKA Right)
    • 行主序与列主序对比
    • 用户自定义步长 (User-Defined Strides)
    • 布局 (Layouts) 的能力

  • std::mdspan 与现有矩阵库的兼容性

  • 自定义矩阵类与 std::mdspan 的集成
  • 使用 std::mdspan 进行多维操作和并行计算
  • submdspan 功能
  • std::mdspan 和 std::submdspan 示例

• 线性代数并行算法示例

约定 (Conventions)

命名空间别名

幻灯片首先定义了一系列命名空间别名，以简化代码书写：
* namespace stdv = std::views;
* namespace stdr = std::ranges;
* namespace ex = std::execution;
* namespace this_thread = std::this_thread;
这些约定旨在提高C++并行编程代码的可读性和简洁性。
（Page 3）

类模板参数推导 (Class Template Argument Deduction - CTAD)

幻灯片还提到了CTAD，它允许编译器在创建类模板实例时推断模板参数，从而减少冗余代码。例如：
* std::tuple t{3.14, 42}; → std::tuple<double, int>
* std::array a{0, 1, 1, 0}; → std::array<int, 4>
（Page 4）

C++ 可移植性

我们为何需要“加速匝道” (We Need On-Ramps)

C++标准并行性为代码从非并行执行过渡到高性能并行执行提供了关键的“加速匝道”。这体现在以下层次：
* Code that doesn't run in parallel (不并行运行的代码)：基础层。
* C++ Standard Parallelism (C++标准并行性)：作为从非并行代码到高性能并行代码的入口，其性能比非并行代码快10倍。
* Platform Specific Extensions (平台特定扩展)：在标准并行性之上，提供更高的性能，达到“光速车道”。

（Page 5）

什么使得C++可移植？ (What Makes C++ Portable?)

传统上，C++的可移植性由一系列语言特性决定，这些特性在20世纪非常重要。然而，随着技术发展，这些特性在21世纪的重要性发生了变化。

以下表格对比了这些特性在20世纪和21世纪的重要性：

特性	Important in the 20th century? (20世纪是否重要？)	Important in the 21st century? (21世纪是否重要？)
Non-8-bit char	✅	❌
Noncommittal sizeof	✅	❌
Non-2's comp int	✅	❌
Non-IEEE float	✅	❌
Non-endian pointers	✅	❌
Aligned addressing	✅	✅
Segmented memory	✅	✅

（Page 9）

新的可移植性契约 (The New Portability Contract)

为了适应现代计算环境，C++的可移植性契约已扩展，包含以下核心要素：

• 内存模型 (Memory Model)
• 执行模型 (Execution Model)
• 前向进展保证 (Forward Progress Guarantees)
• 并发原语 (Concurrency Primitives)
• 并行原语 (Parallelism Primitives)
• 异步模型 (Asynchrony Model)

这些要素共同定义了C++在多核、异构和分布式系统上实现高性能和可移植并行程序的基础。
（Page 11）

C++标准是描述性的，而非规定性的。 (The C++ Standard is descriptive, not prescriptive.)

C++标准为不同平台提供了足够的通用性，以确保可移植性和一致性。同时，它也给予了每个平台足够的自由度，以选择最适合自身的设计和实现。

这意味着：
* 标准指定了足够的内容以实现跨平台的可移植性和一致性。
* 标准给予了每个平台足够的自由度来选择正确的设计。

（Page 13）

实现自由 (Implementation Freedom)

C++标准在可移植性和实现自由之间取得了平衡。它为开发者提供了跨平台编写可移植代码的基础，同时允许不同的平台根据其特定架构和优化需求进行定制化的实现。

幻灯片指出，实现定义和未定义行为通常是一种特性，而非缺陷。C++ 标准更像是“指导原则”而非严格规则，这赋予了实现者一定的自由度。

• 实现定义行为与未定义行为: 这些在C++中通常被视为一种特性，而非缺陷，为编译器和库实现者提供了优化和适应特定平台的机会。
• C++ 标准的性质: C++ 标准更像是“指导原则”，而不是严格的强制规则，这意味着存在一定程度的实现自由。

（Page 15, 17）

C++ 标准行为与并行性

标准并行性位于库中 (Standard Parallelism is in the Library)

C++ 的标准并行性功能主要通过标准库提供。

• C++ 标准并行性: 标准库提供了实现并行化的功能。
  (images/page-0018.jpg "Page 18")
  (images/page-0019.jpg "Page 19")

C++ 标准并行性的支柱 (Pillars of C++ Standard Parallelism)

C++ 标准并行化基于三个主要支柱：

1. 分发到供应商优化并行库的常见算法 (Common Algorithms that Dispatch to Vendor-Optimized Parallel Libraries)
   这些是标准库中已经提供的算法，它们可以利用底层供应商优化的并行库来实现高效的并行执行。这使得开发者能够利用底层硬件的并行能力，而无需关心具体实现细节。
   

2. 编写可在任何地方运行的并行算法的工具 (Tools to Write Your Own Parallel Algorithms that Run Anywhere)
   C++ 提供了工具和机制，允许开发者编写自己的并行算法，并且这些算法可以在不同的硬件和系统上运行。Senders & Receivers 机制提供了一套工具，让开发者能够编写自己的并行算法，并且这些算法可以在任何地方运行（跨不同的执行环境和设备）。示例代码展示了如何使用 sender 和 bulk 操作来构建并行计算任务。
   

3. 将并行调用组合成任务图的机制 (Mechanisms for Composing Parallel Invocations into Task Graphs)
   C++ 标准并行性支持将多个并行操作组合成任务图，以实现更复杂的并行工作流和依赖管理。Senders & Receivers 还提供了将独立的并行调用组合成复杂任务图的机制。这使得更高级的并行程序设计和优化成为可能，例如任务依赖管理和调度。
   
   

C++ 标准算法

串行算法 (C++98)

C++98 提供了多种串行算法，例如：

• std::vector<T> x{...};
• std::for_each(begin(x), end(x), f);
• std::for_each(begin(x), end(x), g);
• std::for_each(begin(x), end(x), h);
  (images/page-0023.jpg "Page 23")

更多串行算法示例：

• std::vector<double> x{...}, y{...};
• double dot_product = std::transform_reduce(begin(x), end(x), begin(y));
• std::span<std::string_view> s{...};
• std::sort(begin(s), end(s));
• std::unordered_map<std::string_view, int> db{...};
• std::vector<std::pair<std::string_view, int>> m{...};
• std::copy_if(begin(db), end(db), begin(m), [](auto e) { return e.second > 0; });
  (images/page-0024.jpg "Page 24")

C++ 标准算法概览（部分列表）：

并行算法 (C++17)

C++17 引入了执行策略（Execution Policy），使得许多标准算法可以并行执行。通过添加执行策略，如 ex::par unseq，可以指示算法以并行且无序的方式执行。

串行与并行 std::for_each 对比：

• 串行 (C++98):

std::vector<T> x{...};
std::for_each(begin(x), end(x), f);
std::for_each(begin(x), end(x), g);
std::for_each(begin(x), end(x), h);

• 并行 (C++17):

std::vector<T> x{...};
std::for_each(ex::par unseq, begin(x), end(x), f);
std::for_each(ex::par unseq, begin(x), end(x), g);
std::for_each(ex::par unseq, begin(x), end(x), h);

更多并行算法示例：

• double dot_product = std::transform_reduce(ex::par unseq, begin(x), end(x), begin(y));
• std::sort(ex::par unseq, begin(s), end(s));
• std::copy_if(ex::par unseq, begin(db), end(db), begin(m), [](auto e) { return e.second > 0; });
  (images/page-0027.jpg "Page 27")

C++ 标准库执行策略 (Execution Policy)

执行策略定义了算法操作的发生方式和顺序。

下表概述了 C++ 标准库中不同的执行策略，包括操作的发生位置和操作的序列性。

执行策略	Operations occur ...	Operations are ...	操作发生位置	操作序列性
std::execution::seq	In the calling thread	Indeterminately sequenced	在调用线程中发生	序列不确定（Indeterminately sequenced）
std::execution::unseq	In the calling thread	Unsequenced	在调用线程中发生	非序列化（Unsequenced）
std::execution::par	(Not shown)	(Not shown)	可能在多个线程中发生	在每个线程中序列不确定
std::execution::par_unseq	(Not shown)	(Not shown)	可能在多个线程中发生	非序列化（Unsequenced）

(images/page-0028.jpg "Page 28")
(images/page-0029.jpg "Page 29")
(images/page-0030.jpg "Page 30")

单词计数示例（使用 std::transform_reduce 和并行策略）

以下系列幻灯片逐步展示了如何使用 std::transform_reduce 结合 ex::par_unseq 策略来实现一个并行单词计数函数。

word_count 函数签名和示例字符串

一个简单的 word_count 函数签名，用于计算 std::string_view 中的单词数量，并定义了一个名为 frost 的 std::string_view 作为示例输入。

std::size_t word_count(std::string_view s) {
  // ...
}

std::string_view frost = "Whose woods these are I think I know.
"
                         "His house is in the village though;
"
                         "He will not see me stopping here   
"
                         "To watch his woods fill up with snow.
";

word_count(frost);

（Page 33）

引入 std::transform_reduce

word_count 函数的实现开始使用 std::transform_reduce 并指定 ex::par_unseq 并行执行策略。

std::size_t word_count(std::string_view s) {
  if (s.empty()) return 0;
  return std::transform_reduce(ex::par_unseq, ...);
}

std::string_view frost = "Whose woods these are I think I know.
"
                         "His house is in the village though;
"
                         "He will not see me stopping here   
"
                         "To watch his woods fill up with snow.
";

word_count(frost);

（Page 34）

std::transform_reduce 的迭代器范围

指定了 std::transform_reduce 的迭代器范围：从 begin(s) 到 end(s) - 1，以及第二个序列的起始迭代器 begin(s) + 1。

std::size_t word_count(std::string_view s) {
  if (s.empty()) return 0;
  return std::transform_reduce(ex::par_unseq,
                               begin(s), end(s) - 1, begin(s) + 1,
                               ...);
}

std::string_view frost = "Whose woods these are I think I know.
"
                         "His house is in the village though;
"
                         "He will not see me stopping here   
"
                         "To watch his woods fill up with snow.
";

word_count(frost);

（Page 35-36）

转换函数 Lambda

引入了用于 transform_reduce 的转换函数 lambda [](char l, char r) { return std::isspace(l) && !std::isspace(r); }，它在遇到从空白字符到非空白字符的转换时返回 true，这可以用于识别单词的开头。

std::size_t word_count(std::string_view s) {
  if (s.empty()) return 0;
  return std::transform_reduce(ex::par_unseq,
                               begin(s), end(s) - 1, begin(s) + 1,
                               ...,
                               [](char l, char r) { return std::isspace(l) && !std::isspace(r); }
                               );
}

std::string_view frost = "Whose woods these are I think I know.
"
                         "His house is in the village though;
"
                         "He will not see me stopping here   
"
                         "To watch his woods fill up with snow.
";

word_count(frost);

（Page 37）

初始结果的二进制表示

展示了在 transform_reduce 执行后的一个中间结果，可能是转换函数输出的布尔值序列（1代表单词开头，0代表其他）。

（Page 38）

修正初始值

修正了 transform_reduce 的初始值，以正确处理字符串开头是单词的情况：std::size_t(!std::isspace(s.front())) ? 1 : 0,。

std::size_t word_count(std::string_view s) {
  if (s.empty()) return 0;
  return std::transform_reduce(ex::par_unseq,
                               begin(s), end(s) - 1, begin(s) + 1,
                               std::size_t(!std::isspace(s.front())) ? 1 : 0,
                               ...,
                               [](char l, char r) { return std::isspace(l) && !std::isspace(r); }
                               );
}

std::size_t result = 00000100000100001010000010100000
                     100010000100100100010000000100000000
                     100100001000100100000000100000000000
                     10010000100010000010000100100000;

（Page 39）

归约操作 std::plus() 和最终计数

添加了 std::plus() 作为归约操作，将所有识别出的单词开头（即转换函数返回 true 的位置）累加起来，从而得到最终的单词计数。

std::size_t word_count(std::string_view s) {
  if (s.empty()) return 0;
  return std::transform_reduce(ex::par_unseq,
                               begin(s), end(s) - 1, begin(s) + 1,
                               std::size_t(!std::isspace(s.front())) ? 1 : 0,
                               std::plus(),
                               [](char l, char r) { return std::isspace(l) && !std::isspace(r); }
                               );
}

std::size_t result = 1 + 1 + 1 + 1 + 1+1 + 1+1 +
                     1 + 1 + 1+1 +1 + 1 + 1 + 1 +
                     1+1 + 1 + 1+1 + 1 +
                     1+1 + 1 + 1 + 1 + 1+1 + 1 ;

（Page 40-41）

C++20 Ranges 介绍

C++20 标准库引入了范围（ranges）的概念，它与迭代器不同，具有可组合性（composable）且可以是惰性的（lazy）。
（Page 42）

并行算法与数据结构示例

std::for_each 和 Ranges 概念图（std::vector）

展示了如何使用 std::for_each 和 ex::par_unseq 策略来处理 std::vector。图示说明了输入元素如何被独立的 lambda 操作并行处理。

std::vector x{...};

std::for_each(
    ex::par_unseq,
    begin(x), end(x),
    [...] (auto& obj) { ... });

（Page 43）

std::for_each 和 Ranges 概念图（std::iota）

展示了如何使用 std::iota 创建一个序列，并用 std::for_each 结合 ex::par_unseq 策略对其进行并行处理。

auto v = std::iota(1, N);

std::for_each(
    ex::par_unseq,
    begin(v), end(v),
    [...] (auto idx) { ... });

（Page 44）

使用 Ranges 进行矩阵操作

此代码片段展示了如何结合 std::span、std::cartesian_product 和 std::iota 来高效地进行矩阵样式的操作，并使用 std::for_each 和 ex::par_unseq 策略并行处理这些操作。

std::span A{input, N * M};
std::span B{output, M * N};

auto v = std::cartesian_product(
    std::iota(0, N),
    std::iota(0, M));

std::for_each(ex::par_unseq,
              begin(v), end(v),
              [=] (auto idx) {
                auto [i, j] = idx;
                B[i + j * N] = A[i * M + j];
              });

本页展示了使用 std::span 和 stdv::cartesian_product 的C++代码片段，并用图示说明了输入数据和操作流程。
- 代码功能: 定义了输入和输出的 std::span，并使用 stdv::cartesian_product 生成一个包含 (0,0) 到 (N-1, M-1) 笛卡尔积的 v。
- 操作图示: 显示了一个二维的 Input 网格，其中每个元素都表示为 (row, column)。下方是针对这些输入元素执行的 Operations，每个操作用一个圆圈内的 T 表示，暗示对每个元素或组合进行处理。

（Page 45, 46）

并行 for_each 阻塞行为

本页展示了对 std::vector 进行两次独立 std::for_each 操作的代码，均使用了 ex::par_unseq 并行策略。
- 代码功能: 对向量 x，首先对每个元素应用函数 f，然后再次对每个元素应用函数 g。ex::par_unseq 提示这些操作是并行且无序的。
- 操作图示: 显示了 Input 向量 x[0], x[1], x[2]...。Operations 部分清晰地展示了两层独立的并行操作：第一层是所有 f 操作，第二层是所有 g 操作。这暗示了两个独立的并行执行阶段。

此页与Page 47的代码相同，但进一步强调了两次 std::for_each 调用之间的阻塞性质。
- 代码功能: 同Page 47。
- 操作图示: 除了Page 47的并行操作图示外，右侧添加了“Blocked”标签和表示阻塞的波浪线，明确指出第一次 std::for_each 完成后，第二次 std::for_each 才能开始，尽管每个 for_each 内部的操作是并行的。

（Page 47, 48）

变换与并行 for_each 操作

本页展示了使用 stdv::transform 和 std::for_each 的C++代码，实现了一种链式或流水线式的并行操作。
- 代码功能: 对向量 x，首先使用 stdv::transform 将函数 f 应用于每个元素，并将结果存储在 v 中。然后，对 v 中的每个元素应用函数 g。
- 操作图示: Input 向量 x[0], x[1], x[2]...。Operations 部分展示了一个“流水线”结构：f 应用于 x[i] 的结果直接输入到 g 应用于该结果。这意味着对于每个独立的输入元素，f 和 g 是按序执行的，但不同元素之间的 f->g 链可以并行执行。

（Page 49）

过滤与归约操作

本页展示了使用 stdv::filter 和 std::reduce 的C++代码，实现了并行过滤和归约操作。
- 代码功能: 对向量 x，首先使用 stdv::filter 根据一个Lambda表达式（[](auto e) { return e > 0; }）筛选出所有大于0的元素，并将结果存储在 v 中。然后，使用 std::reduce 对 v 中的元素执行归约操作（例如求和）。
- 操作图示: Input 向量 x[0], x[1], x[2]...。Operations 部分展示了两个阶段：
1. 过滤阶段: 每个元素都经过一个谓词 P，不满足条件的元素被“移除”（图示中没有后续操作），满足条件的元素则进入下一阶段。
2. 归约阶段: 满足条件的元素（图示中为 x[0] 和 x[2]）通过 + 操作进行归约。

（Page 50）

C++ 并行应用案例

LULESH

LULESH是一个用于非结构化网格上的拉格朗日显式冲击流体动力学的小型应用程序。它旨在强调向量化、并行开销和节点内并行性。该应用大约有9000行C++代码。它存在多种并行版本，包括MPI、OpenMP、OpenACC、CUDA、RAJA、Kokkos和标准C++。

• 性能加速比:

  • OpenMP C++ GCC (64c EPYC 7742): 1.00
  • OpenMP C++ NVC++ (64c EPYC 7742): 1.03
  • Standard C++ GCC + TBB (64c EPYC 7742): 1.53
  • Standard C++ NVC++ (64c EPYC 7742): 2.08
  • Standard C++ NVC++ A100: 13.57
  • 所有使用"Standard C++ Code"的实现均大幅超越了基线OpenMP C++ GCC版本，特别是在A100 GPU上实现了13.57倍的加速。

• 项目链接: https://github.com/LLNL/LULESH
  
  （Page 51）

STLBM

STLBM是一个用于多目标（包括多核CPU和GPU）并行Lattice-Boltzmann模拟的框架。它使用C++标准并行性实现，不依赖语言扩展、外部库、供应商特定代码注释或预编译步骤。
- 碰撞模型加速比:
- 2s 20c Xeon 6148: 1.00
- A100: 12.30
- 使用相同的标准C++代码，在A100上实现了12.30倍的加速。

• 引用: "我们很高兴地发现NVIDIA 'stdpar' 对C++标准并行算法的实现……结果产生了最先进的性能，具有高度的教学意义，并在跨平台CPU/GPU编程领域引入了范式转变。" — Jonas Latt教授，日内瓦大学。
• 项目链接: https://gitlab.com/unigehpfs/stlbm
  
  （Page 52）

M-AIA

M-AIA是一个用于航空航天流体和噪声模拟的软件包。求解器包括有限体积、Navier-Stokes和Lattice-Boltzmann方法。它正从OpenMP切换到C++标准并行性。
- 加速比:
- OpenMP C++ GCC (64c EPYC 7742): 1
- Standard C++ NVC++ (64c EPYC 7742): 0.9
- Standard C++ NVC++ A100: 5
- 在A100 GPU上实现了5倍的加速。

（Page 53）

异步模型: Senders & Receivers

C++17中引入的C++并行算法非常出色，但这仅仅是故事的开始。
（Page 54）

本页展示了对 std::vector<std::string_view> 执行 std::sort 和 std::unique 操作的代码，均使用了 ex::par_unseq 并行策略。

• 代码功能: 首先并行地对字符串视图向量 s 进行排序，然后并行地移除重复项，保持唯一性。
  （Page 55）

当前 C++ 现状与解决方案

当前C++的现状：
- 没有标准的异步模型。
- 没有标准的方式来表达事物应该在哪里执行。
解决方案即将到来：Senders & Receivers。
（Page 58-59）

C++ 并发编程示例

本页展示了一个使用 ex::scheduler、ex::sender 和 ex::then 实现异步任务调度的C++代码示例。
- 代码功能:
1. ex::scheduler auto sch = thread_pool.scheduler();: 从线程池获取一个调度器。
2. ex::sender auto begin = ex::schedule(sch);: 创建一个发送器，表示在调度器上调度一个初始任务。
3. ex::sender auto hi = ex::then(begin, [] { return 13; });: 创建一个发送器 hi，它在 begin 完成后执行一个Lambda函数并返回 13。
4. ex::sender auto add = ex::then(hi, [](int a) { return a + 42; });: 创建一个发送器 add，它在 hi 完成后执行一个Lambda函数，将 hi 的结果 a 加上 42。
5. auto [i] = this_thread::sync_wait(add).value();: 同步等待 add 发送器完成，并获取其结果。

• 这个示例展示了如何链式地创建和执行异步任务，任务的调度和执行通过发送器和接收器模型进行管理。
  （Page 60）

在这些幻灯片中，展示了一个C++代码片段，演示了基于ex::scheduler和ex::sender的异步编程模型。

ex::scheduler auto sch = thread_pool.scheduler();

ex::sender auto begin = ex::schedule(sch);
ex::sender auto hi    = ex::then(begin, [] { return 13; });
ex::sender auto add   = ex::then(hi, [] (int a) { return a + 42; });

auto [i] = this_thread::sync_wait(add).value();

该示例展示了如何：
* 从一个线程池获取一个调度器（sch）。
* 使用调度器ex::schedule(sch)创建一个初始的发送者（begin）。
* 使用ex::then将任务链接起来：
* hi发送者在前一个任务begin完成后执行一个lambda表达式，返回13。
* add发送者在前一个任务hi完成后执行一个lambda表达式，接收hi的返回值a并计算a + 42。

• this_thread::sync_wait(add).value()阻塞当前线程，等待add发送者的结果，并获取其值。
  （Page 61-65）

核心概念定义

• 调度器（Schedulers）：是执行上下文的句柄。
• 发送者（Senders）：表示异步工作。
• 接收者（Receivers）：处理异步信号。
  （Page 66-68）

执行上下文与调度器

调度器是执行上下文的句柄。
（Page 69）

执行上下文示例：
* CPU线程池（CPU Thread Pool）：

(Page 70)
* GPU流（GPU Stream）：

(Page 71)
* 当前线程（Current Thread）：

(Page 72)

调度器与执行上下文的关联：
不同的调度器可以对应不同的执行上下文，例如CPU线程池、GPU流和当前线程。

(Page 73)

调度器之间的交互：
调度器可以在不同的执行上下文之间或在同一个执行上下文内部进行工作调度和通信。

(Page 74)

调度器支持复杂的任务流和依赖管理，允许工作在不同的上下文之间无缝地流转。

(Page 75)

• 调度器生成发送器。
• 发送器代表异步工作。
• 发送器构成任务图的节点。

• 

• 发送器是惰性的（lazy）。
• 当发送器的工作完成时，它会向与其连接的接收器发送信号。

• 

• 接收器（Receiver）是用于接收异步工作结果的组件。

• 

• 接收器可以通过以下三种方式接收信号：

  • set_value(values...)：接收成功的值。
  • set_error(error)：接收发生的错误。
  • set_done()：接收完成信号。

• 

• 异步操作的生命周期包括：

  1. 调度器（Scheduler）通过 schedule 操作生成发送器（Sender）。
  2. 发送器与接收器（Receiver）通过 connect 建立连接。
  3. 连接后生成一个操作状态（Operation State）。
  4. 操作状态通过 start 启动。
  5. 经过一段时间，操作完成。
  6. 操作完成后，通知接收器结果。

• 

• 示例函数签名：sender auto f(sender auto p, ...);，表明发送器可以在函数之间传递和作为返回值。
  （Page 76-89）

发送器函数属性

发送器函数（sender auto f(sender auto p, ...);）具有以下属性：

• 接受一个或多个发送器作为参数。
• 返回一个发送器。
• 可管道化（类比*nix shell）：snd | f | g 等同于 g(f(snd))。
  （Page 90-93）

发送器链式操作示例

以下代码示例演示了如何使用发送器进行链式异步操作，包括数据传输、排序、去重和打印。

示例1：嵌套式链式调用

std::vector<std::string_view> v{...};

ex::sender auto s = for_each_async(
    ex::transfer(
        unique_async(
            sort_async(
                ex::transfer_just(gpu_stream_scheduler{}, v)
            )
        ),
        thread_pool.scheduler()
    ),
    [] (std::string_view e)
    { std::print(file, "{}
", e); }
);

this_thread::sync wait(s);

（Page 94）

示例2：分解为中间发送器

std::vector<std::string_view> v{...};

ex::sender auto s0 = ex::transfer_just(gpu_stream_scheduler{}, v);
ex::sender auto s1 = sort_async(s0);
ex::sender auto s2 = unique_async(s1);
ex::sender auto s3 = ex::transfer(s2, thread_pool.scheduler());
ex::sender auto s4 = for_each_async(s3,
    [] (std::string_view e)
    { std::print(file, "{}
", e); }
);

this_thread::sync wait(s);

（Page 95）

示例3：使用管道操作符

std::vector<std::string_view> v{...};

ex::sender auto s = ex::transfer_just(gpu_stream_scheduler{}, v)
    | sort_async
    | unique_async
    | ex::transfer(thread_pool.scheduler())
    | for_each_async([] (std::string_view e)
        { std::print(file, "{}
", e); });

this_thread::sync wait(s);

（Page 96）

发送器适配器 (Sender Adaptors)

下表总结了常用的发送器适配器及其返回发送器的语义。

Sender Adaptor	Semantics Of Returned Sender
then(sender auto last, invocable auto f)	用 last 发送的值调用 f。
bulk(sender auto last, shape auto n, invocable auto body)	对 n 中的每个索引，用 last 发送的值调用 body。
transfer(sender auto last, scheduler auto sch)	为下一个发送器切换到 sch。
split(sender auto last)	可以连接到多个接收器。
when_all(sender auto... inputs)	将多个发送器组合成一个聚合。
ensure_started(sender auto last)	连接并启动 last。

（Page 97-101）

发送器工厂 (Sender Factories)

下表总结了发送器工厂及其返回发送器的语义。

Sender Factories	Semantics Of Returned Sender
schedule(scheduler auto sch)	在 sch 上完成。
just(T&&... ts)	发送值 ts。

（Page 102-104）

Sender 消费者

• sync_wait(sender auto snd): 阻塞直到 snd 完成，并返回或抛出其发送的任何值。
  • 返回：values-sent-by-sender
    （Page 105-108）

Sender/Receiver 链式操作数据流示例

幻灯片通过一个示例展示了 before | then(f) | after; 这样的链式操作，并逐步揭示其实现细节和数据流。

• 概念链式操作

before | then(f) | after;

（Page 109）

• Sender 的声明与构成

sender auto before_snd = ...;
sender auto then_f_snd = then_sender(before_snd, f);
sender auto after_snd = after_sender(then_f_snd);

上图直观地展示了 Sender 的嵌套结构：`after` 包含 `then(f)`，而 `then(f)` 又包含 `before`。 （Page 110-111）

• 连接操作的顺序

return connect(after_snd, ...);
  return connect(then_f_snd, after_rcv);
    return connect(before_snd, then_f_rcv);

连接操作从最外层的 \1 开始，递归地向内连接。 （Page 112）

• Receiver 的反向构成
  
  与 Sender 的构成相反，Receiver 的嵌套顺序是 before 包含 then(f)，而 then(f) 又包含 after。这反映了数据流的反向路径。
  （Page 113）
• 值传递过程

set_value(before_rcv, ...);
  set_value(then_f_rcv, before_val);
    set_value(after_rcv, f(before_val));

值从最内层的 \1 开始传递，经过 \1 的处理（应用函数 \1），最终到达 \1。 （Page 114）

inclusive_scan_async 适配器实现

本节展示了一个 inclusive_scan_async 的 sender adaptor 的逐步实现。

• 基本结构
  一个 inline constexpr sender adaptor，它是一个 lambda 表达式，返回一个 ex::sender auto。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (...) -> ex::sender auto {
}

（Page 115）

• 添加参数
  为 inclusive_scan_async 引入参数：last (ex::sender), init (初始值), tile_count (分块计数)。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
}

（Page 116）

• then 操作的引入
  在 inclusive_scan_async 的 lambda 体中，使用 return last | ex::then(...) 结构，定义了一个接受 std::random_access_range 类型 input 的 lambda 函数。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
    | ex::then([=] (std::random_access_range auto input) {
    })
}

（Page 117）

• then 内部逻辑
  在 ex::then 内部，创建了一个 std::vector partials 用于存储部分结果，并设置其第一个元素为 init。然后通过 send_values 返回 input 和 partials。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
    | ex::then([=] (std::random_access_range auto input) {
        std::vector<std::range_value_t<decltype(input)>> partials(tile_count + 1);
        partials[0] = init;
        return send_values(input, std::move(partials));
      })
}

（Page 118, 119）

• bulk 操作的引入
  在 ex::then 之后，链式添加了 ex::bulk 操作。它接受 tile_count 和一个 lambda 函数，该 lambda 函数接收 i (索引), input 和 partials 作为参数。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
    | ex::then([=] (std::random_access_range auto input) {
        std::vector<std::range_value_t<decltype(input)>> partials(tile_count + 1);
        partials[0] = init;
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
      [=] (std::size_t i, auto input, auto partials) {
      })
}

（Page 120）

异步包含扫描适配器代码

本页展示了 inclusive_scan_async 函数的初始骨架，它使用 ex::sender 适配器。该函数定义了两个主要阶段：
- ex::then: 用于初始化操作，将 partials 数组的第一个元素 partials[0] 设置为 init，并发送值。
- ex::bulk: 用于并行处理数据块。它计算每个块（tile）的大小、起始和结束位置。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
    | ex::then([=] (std::random_access_range auto input) {
        std::vector<std::range_value_t<decltype(input)>> partials(tile_count + 1);
        partials[0] = init;
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
        [=] (std::size_t i, auto input, auto partials) {
          auto tile_size = (input.size() + tile_count - 1) / tile_count;
          auto start     = i * tile_size;
          auto end       = std::min(input.size(), (i + 1) * tile_size);
        })
  ;
};

（Page 121）

本页通过视觉方式展示了并行 inclusive_scan 的分块处理（tiling）概念。一个输入序列 [a, b, c, d, e, f, g, h, i] 被划分为三个块。每个块独立地执行其内部的 inclusive_scan 操作，产生部分结果：
- 块1：[a, ab, abc]
- 块2：[d, de, def]
- 块3：[g, gh, ghi]

（Page 122）

本页明确指出每个块内的部分 inclusive_scan 是通过 std::inclusive_scan 完成的。这进一步确认了每个数据块都是独立进行扫描。

（Page 123）

本页更新了 inclusive_scan_async 代码，在 ex::bulk 阶段为每个块执行 std::inclusive_scan。
每个块（由 start 和 end 定义）的 inclusive_scan 结果会存储回 input 范围内的相应位置。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
    | ex::then([=] (std::random_access_range auto input) {
        std::vector<std::range_value_t<decltype(input)>> partials(tile_count + 1);
        partials[0] = init;
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
        [=] (std::size_t i, auto input, auto partials) {
          auto tile_size = (input.size() + tile_count - 1) / tile_count;
          auto start     = i * tile_size;
          auto end       = std::min(input.size(), (i + 1) * tile_size);
          std::inclusive_scan(begin(input) + start, begin(input) + end, begin(input) + start);
        })
  ;
};

（Page 124）

本页展示了在并行 inclusive_scan 中，后续块如何依赖前一个块的最终累积和。
- 第一个块的结果 abc 需要用于更新第二个块。
- 第二个块的结果 def 需要用于更新第三个块。
这强调了并行扫描后进行修正的重要性，因为各个块是独立计算的。

（Page 125）

本页对 ex::bulk 阶段的代码进行了微调。现在，它获取每个块的 inclusive_scan 结果的最后一个元素（即该块的累积总和）。*--std::inclusive_scan(...) 表示获取 inclusive_scan 结果范围的最后一个元素。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
    | ex::then([=] (std::random_access_range auto input) {
        std::vector<std::range_value_t<decltype(input)>> partials(tile_count + 1);
        partials[0] = init;
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
        [=] (std::size_t i, auto input, auto partials) {
          auto tile_size = (input.size() + tile_count - 1) / tile_count;
          auto start     = i * tile_size;
          auto end       = std::min(input.size(), (i + 1) * tile_size);
          *--std::inclusive_scan(begin(input) + start, begin(input) + end, begin(input) + start);
        })
  ;
};

（Page 126）

本页进一步修改了代码，将每个块的最终累积和存储到 partials 数组中。partials[i + 1] 接收了第 i 个块的累积和。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
    | ex::then([=] (std::random_access_range auto input) {
        std::vector<std::range_value_t<decltype(input)>> partials(tile_count + 1);
        partials[0] = init;
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
        [=] (std::size_t i, auto input, auto partials) {
          auto tile_size = (input.size() + tile_count - 1) / tile_count;
          auto start     = i * tile_size;
          auto end       = std::min(input.size(), (i + 1) * tile_size);
          partials[i + 1] = *--std::inclusive_scan(begin(input) + start, begin(input) + end, begin(input) + start);
        })
  ;
};

（Page 127）

本页可视化了 partials 数组的形成过程。
- abc 是第一个块的最终累积和。
- def 是第二个块的最终累积和。
- ghi 是第三个块的最终累积和。
这些值被收集到一个名为 partials 的数组中，其内容为 [abc, def, ghi]。

（Page 128）

本页在 partials = [abc, def, ghi] 的下方显示 std::inclusive_scan。这表明接下来将对 partials 数组本身执行一次 inclusive_scan 操作，以计算全局的块级前缀和。

（Page 129）

本页添加了一个新的 ex::then 阶段，它对 partials 数组执行 std::inclusive_scan。这确保 partials 数组中的值变为全局正确的块结束前缀和。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
    | ex::then([=] (std::random_access_range auto input) {
        std::vector<std::range_value_t<decltype(input)>> partials(tile_count + 1);
        partials[0] = init;
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
        [=] (std::size_t i, auto input, auto partials) {
          auto tile_size = (input.size() + tile_count - 1) / tile_count;
          auto start     = i * tile_size;
          auto end       = std::min(input.size(), (i + 1) * tile_size);
          partials[i + 1] = *--std::inclusive_scan(begin(input) + start, begin(input) + end, begin(input) + start);
        })
    | ex::then([] (auto input, auto partials) {
        std::inclusive_scan(begin(partials), end(partials), begin(partials));
      })
  ;
};

（Page 130）

本页在对 partials 数组执行 inclusive_scan 的 ex::then 块中，添加了 return send_values(input, std::move(partials));。这意味着更新后的 partials 数组（现在包含全局正确的块级前缀和）连同原始 input 一起被发送到下一个阶段。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
    | ex::then([=] (std::random_access_range auto input) {
        std::vector<std::range_value_t<decltype(input)>> partials(tile_count + 1);
        partials[0] = init;
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
        [=] (std::size_t i, auto input, auto partials) {
          auto tile_size = (input.size() + tile_count - 1) / tile_count;
          auto start     = i * tile_size;
          auto end       = std::min(input.size(), (i + 1) * tile_size);
          partials[i + 1] = *--std::inclusive_scan(begin(input) + start, begin(input) + end, begin(input) + start);
        })
    | ex::then([] (auto input, auto partials) {
        std::inclusive_scan(begin(partials), end(partials), begin(partials));
        return send_values(input, std::move(partials));
      })
  ;
};

（Page 131）

本页可视化了对 partials 数组执行 std::inclusive_scan 后的结果。
原 partials = [abc, def, ghi] 变成了 partials = [abc, abcdef, abcdefghi]。
- abc 保持不变（因为它是第一个块的累积和）。
- def 变为 abcdef (即 abc + def)。
- ghi 变为 abcdefghi (即 abcdef + ghi)。
这些是每个块末尾的最终全局累积和。

（Page 132）

本页在处理 partials 数组的 ex::then 块之后，添加了一个新的 ex::bulk 块。这个新阶段将使用更新后的 partials 数组来修正每个数据块中的元素。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
    | ex::then([=] (std::random_access_range auto input) {
        std::vector<std::range_value_t<decltype(input)>> partials(tile_count + 1);
        partials[0] = init;
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
        [=] (std::size_t i, auto input, auto partials) {
          auto tile_size = (input.size() + tile_count - 1) / tile_count;
          auto start     = i * tile_size;
          auto end       = std::min(input.size(), (i + 1) * tile_size);
          partials[i + 1] = *--std::inclusive_scan(begin(input) + start, begin(input) + end, begin(input) + start);
        })
    | ex::then([] (auto input, auto partials) {
        std::inclusive_scan(begin(partials), end(partials), begin(partials));
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
        [=] (std::size_t i, auto input, auto partials) {
        })
  ;
};

（Page 133）

本页定义了新 ex::bulk 块的结构，它再次计算 tile_size、start 和 end，为最终的元素修正做准备。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
    | ex::then([=] (std::random_access_range auto input) {
        std::vector<std::range_value_t<decltype(input)>> partials(tile_count + 1);
        partials[0] = init;
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
        [=] (std::size_t i, auto input, auto partials) {
          auto tile_size = (input.size() + tile_count - 1) / tile_count;
          auto start     = i * tile_size;
          auto end       = std::min(input.size(), (i + 1) * tile_size);
          partials[i + 1] = *--std::inclusive_scan(begin(input) + start, begin(input) + end, begin(input) + start);
        })
    | ex::then([] (auto input, auto partials) {
        std::inclusive_scan(begin(partials), end(partials), begin(partials));
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
        [=] (std::size_t i, auto input, auto partials) {
          auto tile_size = (input.size() + tile_count - 1) / tile_count;
          auto start     = i * tile_size;
          auto end       = std::min(input.size(), (i + 1) * tile_size);
        })
  ;
};

（Page 134）

本页完善了最终的 ex::bulk 块，它对每个块内的元素进行迭代，并加上前一个块的全局累积和。std::for_each 循环遍历当前块的元素，并将 partials[i]（即前一个块的最终累积和）加到每个元素 e 上，从而完成全局 inclusive_scan。

inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
    | ex::then([=] (std::random_access_range auto input) {
        std::vector<std::range_value_t<decltype(input)>> partials(tile_count + 1);
        partials[0] = init;
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
        [=] (std::size_t i, auto input, auto partials) {
          auto tile_size = (input.size() + tile_count - 1) / tile_count;
          auto start     = i * tile_size;
          auto end       = std::min(input.size(), (i + 1) * tile_size);
          partials[i + 1] = *--std::inclusive_scan(begin(input) + start, begin(input) + end, begin(input) + start);
        })
    | ex::then([] (auto input, auto partials) {
        std::inclusive_scan(begin(partials), end(partials), begin(partials));
        return send_values(input, std::move(partials));
      })
    | ex::bulk(tile_count,
        [=] (std::size_t i, auto input, auto partials) {
          auto tile_size = (input.size() + tile_count - 1) / tile_count;
          auto start     = i * tile_size;
          auto end       = std::min(input.size(), (i + 1) * tile_size);
          std::for_each(begin(input) + start, begin(input) + end, [&](auto &e) { e += partials[i]; });
        })
  ;
};

（Page 135）

并行扫描

该幻灯片通过视觉方式展示了并行包含扫描（inclusive scan）的操作过程。它将输入数组（从 'a' 到 'i'）划分为多个块，对每个块独立执行 std::inclusive_scan 操作，产生局部偏导数（partials）。然后，对这些局部偏导数再次执行 std::inclusive_scan 以获得累积偏导数。最后，将累积偏导数添加到后续块的局部扫描结果中，以获得完整的包含扫描结果。

（Page 136）

异步包含扫描适配器代码

这两页幻灯片展示了使用 C++ senders/receivers 适配器实现的 inclusive_scan_async 内联 constexpr 代码。该代码定义了一个异步的包含扫描操作，它：
- 首先使用 ex::then 创建一个 std::vector 来存储偏导数（partials），并用 init 初始化 partials[0]。
- 接着使用 ex::bulk 操作将输入数据划分为多个分块（tile），对每个分块执行 std::inclusive_scan，并将每个分块的最后一个元素的扫描结果存储到 partials 数组中。
- 然后再次使用 ex::then 对 partials 数组进行 std::inclusive_scan，以计算累积的偏导数。
- 随后再次使用 ex::bulk 操作，根据计算出的累积偏导数来调整每个分块的扫描结果。具体来说，对于每个分块，它会将之前计算的累积偏导数 partials[i] 加到其所有元素上。
- 最后，使用 ex::then 返回处理后的输入数据。

这两页的内容完全相同。

inline constexpr scan_sender adaptor auto
inclusive_scan_async = [](ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
  return last
  | ex::then([=](std::random_access_range auto input) {
      std::vector<std::range_value_t<decltype(input)>> partials(tile_count + 1);
      partials[0] = init;
      return send_values(input, std::move(partials));
    })
  | ex::bulk(tile_count,
      [=](std::size_t i, auto input, auto partials) {
        auto tile_size = (input.size() + tile_count - 1) / tile_count;
        auto start = i * tile_size;
        auto end = std::min(input.size(), (i + 1) * tile_size);
        partials[i + 1] = *--std::inclusive_scan(
          begin(input) + start,
          begin(input) + end,
          begin(input) + start
        );
      })
  | ex::then([](auto input, auto partials) {
      std::inclusive_scan(partials.begin(), partials.end(), partials.begin());
      return send_values(input, std::move(partials));
    })
  | ex::bulk(tile_count,
      [=](std::size_t i, auto input, auto partials) {
        auto tile_size = (input.size() + tile_count - 1) / tile_count;
        auto start = i * tile_size;
        auto end = std::min(input.size(), (i + 1) * tile_size);
        std::for_each(begin(input) + start,
                      begin(input) + end,
                      [&](auto& e) { e = partials[i] + e; });
      })
  | ex::then([](auto input, auto partials) { return input; });
};

（Page 137, 138）

麦克斯韦方程组 (Maxwell's Equations) 求解

该幻灯片展示了使用 C++ Senders & Receivers 模式来调度和计算麦克斯韦方程组的模拟。代码结构清晰地描述了迭代计算过程：
- maxwell_eqs 函数接收调度器（scheduler）和计算资源（compute, grid_accessor A）作为参数。
- repeat_n 用于执行外层和内层迭代。
- schedule(compute) 表示在计算设备上调度任务。
- bulk(G.cells, update_h(G, hx, hy)) 和 bulk(G.cells, update_e(time, dt, G)) 分别用于批量更新磁场 h 和电场 e。
- halo_exchange(G, hx, hy) 处理边界数据的交换。
- transfer(cpu_serial_scheduler) 可能用于将结果传输回 CPU。
- then(output_results) 在计算完成后输出结果。
右侧图像显示了麦克斯韦方程组模拟的波场传播的可视化结果。

（Page 139）

麦克斯韦方程组性能加速

这一系列幻灯片展示了麦克斯韦方程组模拟从单 CPU 线程扩展到多 CPU 线程、单个 GPU、多个 GPU 乃至 GPU 集群的性能加速情况。所有测试均基于 maxwell_eqs 函数，仅通过改变调度器（scheduler）参数实现不同硬件的调度。

• 单 CPU 串行 (Page 140):

  • 基础性能：Speedup 为 1。
  • 硬件：2x 64c EPYC 7742。
  • 代码：sync_wait(maxwell_eqs(cpu_serial_scheduler), ...);
    

• 多 CPU 并行 (Page 141):

  • 性能：Speedup 达到 92。
  • 硬件：2x 64c EPYC 7742。
  • 代码：sync_wait(maxwell_eqs(cpu_parallel_scheduler), ...);
    

• 单个 GPU (Page 142):

  • 性能：Speedup 达到 468。
  • 硬件：单个 NVIDIA A100 GPU。
  • 代码：sync_wait(maxwell_eqs(single_gpu_scheduler), ...);
    

• 多个 GPU (Page 143):

  • 性能：Speedup 达到 931。
  • 硬件：2x NVIDIA A100 GPU。
  • 代码：sync_wait(maxwell_eqs(multi_gpu_scheduler), ...);
    

• GPU 集群 (Page 144):

  • 性能：随着 A100 GPU 数量从 32 增加到 1120，Speedup 从 3.95 线性扩展到 27.14。
  • 代码：sync_wait(maxwell_eqs(multi_node_gpu_scheduler), ...);
    
    （Page 140-144）

Palabos 碳捕集模拟

该幻灯片展示了 Palabos 框架在模拟砂岩中碳捕集应用上的性能加速。

• 应用描述:

  • Palabos 是一个用于使用 Lattice-Boltzmann 方法进行并行计算流体动力学模拟的框架。
  • 用于通过多孔介质的多组分流动的代码已移植到 C++ Senders and Receivers。
  • 具体应用：模拟砂岩中的碳捕集。

• 性能加速:

  • 随着 A100 GPU 数量的增加，性能呈现线性加速。
  • 使用 32 个 A100 GPU 时，Speedup 约为 2.1。
  • 使用 128 个 A100 GPU 时，Speedup 约为 3.9。
  • 使用 512 个 A100 GPU 时，Speedup 达到 14.6。
    左侧图像展示了多孔介质的结构和流体粒子（红色球体）的模拟。

（Page 145）

C++ 标准算法的演进

该幻灯片对比了 C++ 标准算法在串行、并行和异步执行模型下的演进：

• 串行 (C++98):

  • 传统的 std::vector<T> x{...}; 初始化。
  • 多个 std::for_each(begin(x), end(x), f); 调用是顺序执行的。

• 并行 (C++17):

  • 引入了执行策略（execution policies）。
  • 多个 std::for_each(ex::par_unseq, begin(x), end(x), f); 调用可以使用 ex::par_unseq 策略并行执行，但它们之间仍然是阻塞的，即在一个 for_each 完成前，下一个不会开始。

• 异步 (C++ Senders & Receivers):

  • 引入了 ex::sender 概念，允许构建非阻塞的任务管道。
  • ex::sender s = ex::transfer_just(sch, x) | for_each_async(f) | for_each_async(g) | for_each_async(h); 这行代码创建了一个任务流，其中 f, g, h 可以异步且可能并行地执行。
  • this_thread::sync_wait(s); 用于等待整个异步任务流完成。

（Page 147）

多维数据抽象与 std::mdspan

多维数据抽象的缺失与解决方案

• 现状 (Page 148):

  • 如今，C++ 对多维数据没有合理的抽象。
  • 这意味着处理多维数组（如图像、矩阵、张量）时通常需要手动管理索引和内存布局，增加了复杂性和出错率。

• C++23 的解决方案 (Page 149):

  • 解决方案将在 C++23 中到来：std::mdspan。
  • std::mdspan 是 C++ 标准库中的一个新组件，旨在提供对多维数据的统一视图和操作。

• std::mdspan 特性 (Page 150):

  • std::mdspan 是 非拥有的（Non-owning） 数据结构。
  • 它由一个 指针（pointer） 和 元数据（metadata） 组成。这意味着 std::mdspan 仅仅提供了对现有内存中多维数据的一个“视图”，它不拥有数据本身，而是通过指针指向数据，并通过元数据（如维度、步长）来描述数据的多维布局。
    （Page 148-150）

std::mdspan 介绍

std::mdspan 的核心特性：

• 非拥有：它只包含一个指针和元数据，不拥有其指向的数据。
• 元数据：元数据可以是动态的或静态的。
• 参数化布局和访问：它参数化了数据的布局和访问方式。
  （Page 153）

std::extents 概述

std::extents 是一个模板类，用于描述多维数组的维度信息：

template <std::size_t... Extents>
class std::extents;

（Page 154）

std::extents 的使用示例：

// 模板 <std::size_t... Extents>
// 类 std::extents;

// 示例 1: 动态维度
std::extents e0{16, 32};
// 等价于:
std::extents<std::dynamic_extent, std::dynamic_extent> e1{16, 32};
// 或者使用 std::dextents 别名:
std::dextents<2> e2{16, 32};

// 属性访问:
e0.rank() == 2 // 维度数量
e0.extent(0) == 16 // 第0维的范围
e0.extent(1) == 32 // 第1维的范围

// 示例 2: 静态维度
std::extents<16, 32> e3;

// 示例 3: 混合维度 (静态和动态)
std::extents<16, std::dynamic_extent> e4{32};

// 示例 4: 更多动态维度
std::extents e5{16, 32, 48, 4};

（Page 159）

std::mdspan 模板定义及示例

std::mdspan 的完整模板定义：

template <class I,
          class Extents,
          class LayoutPolicy = std::layout_right,
          class AccessorPolicy = std::default_accessor<T>>
class std::mdspan;

其中：
* I：数据类型（例如 double）。
* Extents：描述多维数组维度信息的 std::extents 对象。
* LayoutPolicy：布局策略，默认为 std::layout_right。
* AccessorPolicy：访问器策略，默认为 std::default_accessor<T>。
（Page 164）

std::mdspan 的使用示例：

std::mdspan m0{data, 16, 32};
// 等价于:
std::mdspan<double, std::dextents<2>> m1{data, 16, 32};

（Page 165）

多维数组布局与 std::mdspan

std::mdspan 基础用法

std::mdspan 是一个模板类，用于提供多维数组视图，支持灵活的布局策略和访问器策略。

template <class I,
          class Extents,
          class LayoutPolicy = std::layout_right,
          class AccessorPolicy = std::default_accessor<T>>
class std::mdspan;

基本的 std::mdspan 构造示例：

std::mdspan m0{data, 16, 32};
// 等效于:
std::mdspan<double, std::dextents<2>> m1{data, 16, 32};

// 对于行主序布局，元素访问方式为:
m0[i, j] == data[i * M + j]

std::mdspan 可以通过 std::extents 或 std::dynamic_extent 来指定维度和大小：

std::mdspan m2{data, std::extents<16, 32>{}};
// 等效于:
std::mdspan<double, std::extents<16, 32>> m3{data};

std::mdspan m4{data, std::extents<16, std::dynamic_extent>{}};

（Page 166, 167）

行主序 (Row-Major AKA Right)

• 在 C++ 和 NumPy 中默认为行主序。
• 最右侧的维度是连续存储的。
• 访问公式和步长 (stride) 示例：

mdspan A{data, N, M};
mdspan A{data, layout_right::mapping(N, M)};

A[i, j] == data[i * M + j]
A.stride(0) == M
A.stride(1) == 1

（Page 168, 169）

行主序与列主序对比

• 行主序 (Row-Major AKA Right)

  • 常见于 C++、NumPy (默认)。
  • 最右侧的维度是连续存储的。
  • 访问公式：A[i, j] == data[i * M + j]
  • 步长：A.stride(0) == M, A.stride(1) == 1

• 列主序 (Column-Major AKA Left)

  • 常见于 Fortran、MATLAB。
  • 最左侧的维度是连续存储的。
  • 访问公式：B[i, j] == data[i + j * N]
  • 步长：B.stride(0) == 1, B.stride(1) == N

（Page 170, 171）

用户自定义步长 (User-Defined Strides)

std::mdspan 允许用户通过 layout_stride::mapping 定义自定义步长：

mdspan C{data, layout_stride::mapping(extents(N, M), {X, Y})};

// 元素访问公式:
A[i, j] == data[i * X + j * Y]
// 对应的步长:
A.stride(0) == X
A.stride(1) == Y

（Page 172）

布局 (Layouts) 的能力

布局 (Layouts) 将多维索引 (i,j,k, ...) 映射到数据存储位置。
任何人都可以定义一个布局。

布局可能：

• 是非连续的 (Be non-contiguous)。
• 将多个索引映射到同一个位置 (Map multiple indices to the same location)。
• 执行复杂的计算 (Perform complicated computations)。
• 拥有或引用状态 (Have or refer to state)。

参数化布局 (Parametric layout) 能够实现通用的多维算法。
（Page 173-179）

std::mdspan 与现有矩阵库的兼容性

以下是一个使用 Eigen 库的多维矩阵作为函数参数的示例：

void your_function(Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& m);

（Page 180）

幻灯片展示了 your_function 函数接受 Eigen::Matrix 类型参数的定义和调用示例。

该函数 void your_function(Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& m); 接受一个动态大小的 double 类型 Eigen 矩阵的引用。

为了展示 std::mdspan 的通用性，幻灯片进一步列出了 your_function 如何处理来自不同库的矩阵类型，包括 Eigen、Boost.Ublas、PETSc、Blaze 和 CUTLASS。这表明尽管函数签名特定于 Eigen，但通过适当的包装或转换，多种矩阵类型可以与类似函数进行交互。

• your_function(Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>{...});
• your_function(boost::numeric::ublas::matrix<double>{...});
• your_function(Mat{...}); // PETSc
• your_function(blaze::DynamicMatrix<double, blaze::rowMajor>{...});
• your_function(cutlass::HostTensor<float, cutlass::layout::ColumnMajor>{...});

当 your_function 的参数类型被修改为 std::mdspan 时，它能够以统一的方式处理上述所有不同库的矩阵类型，凸显了 std::mdspan 作为通用多维视图的强大能力。

修改后的函数签名为 void your_function(std::mdspan<T, Extents, Layout, Accessor> m);，这意味着 std::mdspan 可以封装不同的底层数据存储和布局，使其能够兼容各种矩阵库。
（Page 181-183）

自定义矩阵类与 std::mdspan 的集成

幻灯片定义了一个简单的 my_matrix 结构体，它包含行数、列数和使用 std::vector<double> 存储的扁平化数据。该结构体提供了构造函数、 operator() 用于元素访问以及获取行数和列数的方法。

struct my_matrix {
public:
  my_matrix(std::size_t N, std::size_t M)
  : num_rows_(N), num_cols_(M), storage_(num_rows_ * num_cols_) {}

  double& operator()(size_t i, size_t j)
  { return storage_[i * num_cols_ + j]; }
  const double& operator()(size_t i, size_t j) const
  { return storage_[i * num_cols_ + j]; }

  std::size_t num_rows() const { return num_rows_; }
  std::size_t num_cols() const { return num_cols_; }

private:
  std::size_t num_rows_, num_cols_;
  std::vector<double> storage_;
};

为了使 my_matrix 能够与 std::mdspan 无缝集成，幻灯片展示了如何为 my_matrix 添加一个类型转换运算符，使其能够隐式转换为 std::mdspan<double, std::dextents<2>>。这允许用户像使用 std::mdspan 一样使用 my_matrix 实例。

  // ... (previous code) ...
  operator std::mdspan<double, std::dextents<2>> const
  { return {storage_, num_rows_, num_cols_}; }

private:
  // ... (previous code) ...
};

（Page 184, 185）

使用 std::mdspan 进行多维操作和并行计算

幻灯片展示了使用 std::mdspan 进行三维7点模板操作的示例。它定义了输入 A 和输出 B 为三维 mdspan，并利用 stdv::cartesian_product 和 stdv::iota 生成迭代索引。std::for_each(ex::par_unseq, ...) 表示这是一个并行、无序的循环，用于计算 B 中的每个元素，基于其邻近元素 A。

std::mdspan A{input, N, M, 0};
std::mdspan B{output, N, M, 0};

auto v = stdv::cartesian_product(
  stdv::iota(1, A.extent(0) - 1),
  stdv::iota(1, A.extent(1) - 1),
  stdv::iota(1, A.extent(2) - 1));

std::for_each(ex::par_unseq,
  begin(v), end(v),
  [=] (auto idx) {
    auto [i, j, k] = idx;
    B[i, j, k] = ( A[i-1, j, k] +
                   A[i+1, j, k] +
                   A[i, j-1, k] +
                   A[i, j+1, k] +
                   A[i, j, k-1] +
                   A[i, j, k+1] +
                   A[i, j, k] ) / 7.0;
  });

为了明确指定内存布局，幻灯片对 std::mdspan A 和 B 的定义中增加了 std::layout_left::mapping(N, M, 0)。这允许开发者精确控制数据在内存中的存储方式，而模板操作的逻辑保持不变。

std::mdspan A{input, std::layout_left::mapping(N, M, 0)};
std::mdspan B{output, std::layout_left::mapping(N, M, 0)};

// ... (rest of the code is the same as Page 186) ...

幻灯片还展示了使用 std::span 实现矩阵转置的示例。由于 std::span 是一维的，所以索引计算需要手动进行，将二维逻辑映射到一维存储。

std::span A{input, N * M};
std::span B{output, M * N};

auto v = stdv::cartesian_product(
  stdv::iota(0, N),
  stdv::iota(0, M));

std::for_each(ex::par_unseq,
  begin(v), end(v),
  [=] (auto idx) {
    auto [i, j] = idx;
    B[i + j * N] = A[i * M + j];
  });

相比之下，使用 std::mdspan 实现矩阵转置则更为简洁直观，因为它直接支持多维索引。通过 B[j, i] = A[i, j] 即可完成转置操作。

std::mdspan A{input, N, M};
std::mdspan B{output, M, N};

auto v = stdv::cartesian_product(
  stdv::iota(0, A.extent(0)),
  stdv::iota(0, A.extent(1)));

std::for_each(ex::par_unseq,
  begin(v), end(v),
  [=] (auto idx) {
    auto [i, j] = idx;
    B[j, i] = A[i, j];
  });

幻灯片进一步简化了 std::mdspan 矩阵转置的迭代方式，通过使用 A.indices() 结合 stdr::for_each，避免了手动生成 cartesian_product 和 iota。

std::mdspan A{input, N, M};
std::mdspan B{output, M, N};

stdr::for_each(
  ex::par_unseq,
  A.indices(),
  [=] (auto [i, j]) {
    B[j, i] = A[i, j];
  });

最后，幻灯片展示了将 std::mdspan 操作与 C++ Send/Receive 机制结合，实现异步并行转置的示例。

std::mdspan A{input, N, M};
std::mdspan B{output, M, N};

ex::sender auto s =
  ex::transfer_just(sch, A.indices())
  | for_each_async(
      [=] (auto [i, j]) {
        B[j, i] = A[i, j];
      });

（Page 186-191）

submdspan 功能

submdspan 函数提供了一种从现有 mdspan 中提取子视图的功能。

其签名为 submdspan(mdspan<...> m, SliceSpecifiers... ss) -> mdspan<...>，它接受一个 mdspan 实例 m 和一系列切片描述符 ss，并返回一个新的 mdspan。
（Page 192）

幻灯片详细介绍了 submdspan 的切片描述符类型：

• 单索引 (Single Index)：参数为 Integral 类型，会将 mdspan 的维度降低一级。

Slice Specifier	Argument	Reduces Rank?
Single Index	Integral	✅

• 索引范围 (Range of Indices)：参数可以是 std::pair<Integral, Integral> 或 std::tuple<Integral, Integral>，它不会降低 mdspan 的维度。

Slice Specifier	Argument	Reduces Rank?
Single Index	Integral	✅
Range of Indices	std::pair<Integral, Integral> std::tuple<Integral, Integral>	❌

• 所有索引 (All Indices)：参数为 std::full_extent，表示选择该维度的所有元素，同样不会降低 mdspan 的维度。

Slice Specifier	Argument	Reduces Rank?
Single Index	Integral	✅
Range of Indices	std::pair<Integral, Integral> std::tuple<Integral, Integral>	❌
All Indices	std::full_extent	❌

（Page 193-195）

std::mdspan 和 std::submdspan 示例

std::submdspan 基础用法

• 从 std::mdspan 创建子视图
  原始的 std::mdspan m0 被定义为 3 维，其范围分别为 64, 128, 32。

std::mdspan m0{64, 128, 32};
auto m1 = std::submdspan(m0, std::tuple{16, 24},
                             std::tuple{32, 40},
                             std::tuple{ 8, 16});

（Page 196）
`m1` 是 `m0` 的一个子视图，通过为每个维度指定 `std::tuple{start, end}` 范围来创建。

• 子视图的秩 (Rank)
  创建的子视图 m1 保持了与原始 m0 相同的秩。

m1.rank() == 3

（Page 197）

• 子视图的扩展 (Extents)
  m1 的每个维度上的扩展是根据 std::tuple 中定义的范围计算得出的。例如，对于第一个维度，范围是 [16, 24)，因此扩展为 24 - 16 = 8。

m1.extent(0) == 8
m1.extent(1) == 8
m1.extent(2) == 8

（Page 198）

• 子视图的索引映射
  m1 中的索引 [i, j, k] 对应于 m0 中的 [i + offset0, j + offset1, k + offset2]，其中偏移量是创建 std::submdspan 时每个维度范围的起始值。

m1[i, j, k] == m0[i + 16, j + 32, k + 8]

（Page 199）

std::submdspan 结合 std::full_extent 的用法

• 使用 std::full_extent 和单一索引创建子视图
  除了指定范围外，std::submdspan 还可以接受单一索引和 std::full_extent 来创建子视图。

auto m2 = std::submdspan(m0, 16,
                             std::full_extent,
                             32);

（Page 200） 这个示例中，\1 是 \1 在第一个维度固定为索引 \1，第三个维度固定为索引 \1，而第二个维度保持完整范围的子视图。

• 子视图的秩变化
  当维度被一个单一索引固定时，该维度会从子视图的秩中移除。因此，m2 的秩从 m0 的 3 降为 1。

m2.rank() == 1

（Page 201）

• 子视图的扩展
  m2 只有一个维度，其扩展对应于 m0 第二个维度（std::full_extent 对应的维度）的扩展 128。

m2.extent(0) == 128

（Page 202）

• 子视图的索引映射
  m2 中的单一索引 [j] 对应于 m0 中的 [16, j, 32]。

m2[j] == m0[16, j, 32]

（Page 203）

应用示例：分块矩阵转置

此示例展示了如何使用 std::mdspan 和 std::submdspan 实现一个分块（tiled）的矩阵转置操作。

• 初始化输入和输出 std::mdspan
  定义输入 A 和输出 B 矩阵，以及分块大小 T。
  A 是 N x M 矩阵，B 是 M x N 矩阵。

std::mdspan A{input, N, M};
std::mdspan B{output, M, N};
std::size_t T = ...;

（Page 204）

• 定义外部循环：遍历瓦片 (Tiles)
  使用 stdv::cartesian_product 和 stdv::iota 创建一个外层循环，用于遍历 N x M 矩阵的瓦片。 (X + T - 1) / T 计算了在维度 X 上所需瓦片的数量。

auto outer = stdv::cartesian_product(stdv::iota(0, (N + T - 1) / T),
                                     stdv::iota(0, (M + T - 1) / T));

（Page 205）

• 并行处理瓦片
  std::for_each 结合执行策略 ex::par_unseq 用于并行或非序列地处理每个瓦片。

std::for_each(ex::par_unseq, begin(outer), end(outer),
              [&](auto tile) {
                  auto [x, y] = tile;
                  // ...
              });

（Page 206）

• 计算当前瓦片的实际范围
  在每个瓦片内部，计算当前瓦片在原始矩阵 A 中所对应的实际起始和结束索引，使用 std::min 确保不超过矩阵边界。

std::tuple selectN{T * x, std::min(T * (x + 1), N)};
std::tuple selectM{T * y, std::min(T * (y + 1), M)};

（Page 207）

• 创建瓦片的子视图
  使用 std::submdspan 为当前瓦片创建 A 的子视图 TA 和 B 的子视图 TB。注意 B 的子视图范围是转置的 (selectM, selectN)。

auto TA = std::submdspan(A, selectN, selectM);
auto TB = std::submdspan(B, selectM, selectN);

（Page 208）

• 定义内部循环：遍历瓦片内的元素
  使用 stdv::cartesian_product 和 stdv::iota 为内部循环创建迭代器，遍历当前瓦片 TA 的所有元素。

auto inner = stdv::cartesian_product(stdv::iota(0, TA.extent(0)),
                                     stdv::iota(0, TA.extent(1)));

（Page 209）

• 执行瓦片内元素的转置
  通过内层循环遍历瓦片内的每个元素 [i, j]，并将 TA[i, j] 赋值给 TB[j, i]，完成局部转置。

for (auto [i, j] : inner)
    TB[j, i] = TA[i, j];

（Page 210）

（Page 196-210）

线性代数并行算法示例

幻灯片展示了使用C++标准并行性进行线性代数操作的示例：

• 矩阵-向量乘法示例：
  • 使用std::mdspan定义多维视图，例如 std::mdspan A{..., N, M};。
  • 演示了计算 y = 3.0 A x + 2.0 y 的并行实现。
  • 核心代码为 std::matrix_vector_product(ex::par_unseq, std::scaled(3.0, A), x, std::scaled(2.0, y), y);，其中 ex::par_unseq 指示并行非顺序执行。

（Page 212-214）

• 三角矩阵向量求解示例：

  • 目标是求解 A x = b，其中 A = U^T U。
  • 定义了 std::mdspan A, x, b。

  • 第一步：求解 U^T c = b，并将结果存储在 x 中作为 c。

    • std::triangular_matrix_vector_solve(ex::par_unseq, std::transposed(A), std::upper_triangle, std::explicit_diagonal, b, x);

  • 第二步：求解 U x = c，并用结果覆盖 x。

    • std::triangular_matrix_vector_solve(ex::par_unseq, A, std::upper_triangle, std::explicit_diagonal, x, x);

  • 这两个步骤都利用了 ex::par_unseq 进行并行执行，并使用了 std::transposed、std::upper_triangle、std::explicit_diagonal 等辅助函数来指定矩阵属性。

（Page 215）