STANDARD PARALLELISM

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

目录

约定 (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++的可移植性契约已扩展,包含以下核心要素:

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

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

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

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

C++标准与平台
(Page 13)

实现自由 (Implementation Freedom)

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

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

C++标准与实现自由
实现自由
(Page 15, 17)

C++ 标准行为与并行性

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

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

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++ 标准算法

串行算法 (C++98)

C++98 提供了多种串行算法,例如:

更多串行算法示例:

C++ 标准算法概览(部分列表):
C++标准算法列表

并行算法 (C++17)

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

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

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);
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);
串行与并行For_each对比

更多并行算法示例:

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代表其他)。

Result of transform_reduce - binary representation
(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) { ... });

Input and operations for std::vector with std::for_each
(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) { ... });

Input and operations for std::iota with std::for_each
(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 操作。这暗示了两个独立的并行执行阶段。
并行for_each操作

此页与Page 47的代码相同,但进一步强调了两次 std::for_each 调用之间的阻塞性质。
- 代码功能: 同Page 47。
- 操作图示: 除了Page 47的并行操作图示外,右侧添加了“Blocked”标签和表示阻塞的波浪线,明确指出第一次 std::for_each 完成后,第二次 std::for_each 才能开始,尽管每个 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 链可以并行执行。
变换与并行for_each操作
(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++。

STLBM

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

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倍的加速。
M-AIA加速比和可视化
(Page 53)

异步模型: Senders & Receivers

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

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

当前 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 发送器完成,并获取其结果。

在这些幻灯片中,展示了一个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。

核心概念定义

执行上下文与调度器

调度器是执行上下文的句柄。
(Page 69)

执行上下文示例:
* CPU线程池(CPU Thread Pool):
CPU线程池执行上下文
(Page 70)
* GPU流(GPU Stream):
GPU流执行上下文
(Page 71)
* 当前线程(Current Thread):
当前线程执行上下文
(Page 72)

调度器与执行上下文的关联:
不同的调度器可以对应不同的执行上下文,例如CPU线程池、GPU流和当前线程。
调度器与执行上下文的关联
(Page 73)

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

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

发送器函数属性

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

发送器链式操作示例

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

示例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 消费者

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

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

before | then(f) | after;
(Page 109)
sender auto before_snd = ...;
sender auto then_f_snd = then_sender(before_snd, f);
sender auto after_snd = after_sender(then_f_snd);
Sender的声明与构成 Page 111 上图直观地展示了 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)

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 auto
inclusive_scan_async = [] (...) -> ex::sender auto {
}
(Page 115)
inline constexpr sender adaptor auto
inclusive_scan_async = [] (ex::sender auto last, auto init, std::size_t tile_count) -> ex::sender auto {
}
(Page 116)
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)
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)
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]

并行 inclusive_scan 的分块处理
(Page 122)

本页明确指出每个块内的部分 inclusive_scan 是通过 std::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]。

Partial sums 数组的形成
(Page 128)

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

对 partials 数组执行 std::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)。
这些是每个块末尾的最终全局累积和。

扫描 partials 数组后的结果
(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)参数实现不同硬件的调度。

Palabos 碳捕集模拟

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

Palabos碳捕集模拟及其性能
(Page 145)

C++ 标准算法的演进

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

C++标准算法演进
(Page 147)

多维数据抽象与 std::mdspan

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

std::mdspan 介绍

std::mdspan 的核心特性:

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)

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)

行主序与列主序对比

行主序与列主序对比及内存布局
(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, ...) 映射到数据存储位置。
任何人都可以定义一个布局。

布局可能:

参数化布局 (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 类型参数的定义和调用示例。

your_function 接受 Eigen::Matrix 参数
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 调用不同矩阵类型
your_function 调用不同矩阵类型

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

your_function 接受 std::mdspan 参数
your_function 接受 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() 用于元素访问以及获取行数和列数的方法。

自定义 my_matrix 结构体
自定义 my_matrix 结构体
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 实例。

自定义 my_matrix 集成 mdspan
自定义 my_matrix 集成 mdspan
  // ... (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。

三维7点模板操作
三维7点模板操作
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)。这允许开发者精确控制数据在内存中的存储方式,而模板操作的逻辑保持不变。

三维7点模板操作与布局
三维7点模板操作与布局
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 进行矩阵转置
使用 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 进行矩阵转置
使用 std::mdspan 进行矩阵转置
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。

使用 stdr::for_each 进行矩阵转置
使用 stdr::for_each 进行矩阵转置
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 机制结合,实现异步并行转置的示例。

使用 Send/Receive 进行异步矩阵转置
使用 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 函数签名
submdspan 函数签名

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

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

submdspan 切片描述符 (单索引)
submdspan 切片描述符 (单索引)
Slice Specifier Argument Reduces Rank?
Single Index Integral ✅
submdspan 切片描述符 (索引范围)
submdspan 切片描述符 (索引范围)
Slice Specifier Argument Reduces Rank?
Single Index Integral ✅
Range of Indices std::pair<Integral, Integral> std::tuple<Integral, Integral> ❌
submdspan 切片描述符 (所有索引)
submdspan 切片描述符 (所有索引)
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 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}` 范围来创建。
m1.rank() == 3
(Page 197)
m1.extent(0) == 8
m1.extent(1) == 8
m1.extent(2) == 8
(Page 198)
m1[i, j, k] == m0[i + 16, j + 32, k + 8]
(Page 199)

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

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

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

m2.rank() == 1
(Page 201)
m2.extent(0) == 128
(Page 202)
m2[j] == m0[16, j, 32]
(Page 203)

应用示例:分块矩阵转置

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

std::mdspan A{input, N, M};
std::mdspan B{output, M, N};
std::size_t T = ...;
(Page 204)
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, begin(outer), end(outer),
              [&](auto tile) {
                  auto [x, y] = tile;
                  // ...
              });
(Page 206)
std::tuple selectN{T * x, std::min(T * (x + 1), N)};
std::tuple selectM{T * y, std::min(T * (y + 1), M)};
(Page 207)
auto TA = std::submdspan(A, selectN, selectM);
auto TB = std::submdspan(B, selectM, selectN);
(Page 208)
auto inner = stdv::cartesian_product(stdv::iota(0, TA.extent(0)),
                                     stdv::iota(0, TA.extent(1)));
(Page 209)
for (auto [i, j] : inner)
    TB[j, i] = TA[i, j];
(Page 210)

(Page 196-210)

线性代数并行算法示例

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

Matrix-Vector Product Example
(Page 212-214)

Triangular Matrix Vector Solve Example
(Page 215)