让 C++ 标准的并行机制走向多维¶
标题:Making C++ Standard Parallelism Multidimensional
日期:2026/07/12
作者:Mark Hoemmen
链接:https://www.youtube.com/watch?v=VAwW_s1uEHY
注意:此为 AI 翻译生成 的中文转录稿,详细说明请参阅仓库中的 README 文件。
引言:并行与并发的区别¶
大家好,早上好。早上好。今天我们要探讨的话题是:让 C++ 标准的并行机制走向多维。请大家不要担心,我今天带来的并不是什么庞然大物般的巨型提案。不用太紧张。确切地说,这是一个关于“仍在厨房里酝酿的想法”的演讲,我现在正在培养酸面团的酵头,正在揉面团。也就是说,我还没有把这些想法真正烤熟。今天大家并不是来面包店品尝成品的。
在演讲中,我会先谈谈目前我们在 C++ 标准中对并行的支持,不仅是多维的,也包括一般的并行支持。我们现在有并行算法,有数据并行类型,有并行的权限。然后,我会提供一些理由,说明为什么我们应该在 C++ 中支持并行的多维算法。我们为什么要走向多维?假设我们确实需要它,接下来我就会讨论将多维特性与并行特性结合起来会是什么样子。
在开始之前,我需要明确我所说的“并行(parallelism)”究竟指什么。我要把这种对并行和并发(concurrency)之间的对比与区分归功于 Daisy Holman。我之所以选择“并行”,是因为我想要放宽执行的顺序,以换取性能。而对我来说,“并发”意味着恰恰相反的东西:它意味着我想要约束执行的顺序,我这样做是为了保证正确性,因为我的算法需要它。
当我说放宽顺序时,这可能意味着如果我有两件事情,它们彼此之间可以以任何顺序发生,那就是并行。或者它们可以相对于调用者以任何顺序发生,那是异步(asynchrony)。但在我脑海中,这两者都属于并行。这本质上关乎“权限”。而对于并发,有时我确实需要线程。我需要一个大于一的线程数,因为我在对“向前推进(forward progress)”做出某种承诺。这更像是对事物执行方式的顺序约束。相比之下,并行更多的是关于执行的权限,而不一定必须实际拥有多个线程。因此,在讨论并发时,人们会谈论“先行发生(happens before)”。这是一种顺序约束。而并行则完全相反。
标准 C++ 对并行的现有支持¶
目前的标准 C++ 以四种不同的方式支持并行:异步执行(asynchronous execution)、并行算法(parallel algorithms)、数据并行类型和操作(data parallel types and operations),以及数组(arrays)。是的,数组。
对于异步执行,这既可以意味着显式创建在后台运行的线程(我们可以敲入 std::thread、jthread 或 std::async),也可以意味着创建可能在调度器上异步执行的具有依赖关系的异步任务(这涉及 std::execution 或协程)。
其次是并行算法,我们现在有各种各样的并行算法。在 C++17 中,我们引入了基于迭代器的并行算法。感谢 Ruslan,我们在 C++26 中拥有了并行的 ranges 算法(parallel ranges)。C++26 结合 std::execution 还引入了异步并行算法:bulk、非分块的 bulk 和分块的 bulk。它还引入了并行的线性代数算法。此外,目前管道中还有一些针对多维算法的提案。例如用于 mdspan 的 copy 和 fill,这些是在 LEWG 完成设计审查的并行算法,正在等待 LWG 的措辞审查。还有“多维 for_each”提案,它正在等待设计审查,但已经在四月份的邮件列表中提交了。
我在屏幕上展示了一些例子。for_each_index 是拟议中的多维 for_each。矩阵向量乘积(Matrix vector product)是 C++26 线性代数库的一部分。而 ranges::transform 是我们在 C++26 中拥有的并行 ranges 的一部分。
第三点是数据并行类型和操作。也许你并不一定会将其视为 SIMD(单指令多数据流),但我们一直都拥有它们,因为你可以对整数进行按位操作(bitwise operations)。这些就是并行算法,是字(words)上的按位并行算法。它们从一开始就在那里。我们在 C++20 中通过 <bit> 头文件获得了一些新的操作。
然后是真正的 SIMD。它不是字上的按位操作,而是对一组字(字的区块)进行的按元素操作(elementwise operations)。右边展示了一个 SIMD 并行算法的代码示例。你不需要读完所有的代码,关键点在于这些机制自带了类型和操作。操作不仅包括对 SIMD 类型的算术、数学计算,还包括将它们从内存中取出并放回去的方法:加载(loads)、存储(stores)、收集(gathers)、分散(scatters)。
并行算法和 SIMD 结合在一起,开始构建一个并行处理数据的系统。全局算法是全局的。并行算法是全局的,它们处理你的整个数据结构,在最外层的循环上工作。而数据并行类型是局部的。它们作用于你数据的一小块。我们在稍后的演讲中会更多地讨论全局与局部之间的这种交互。但记住这一点很重要:并行算法处理的是整个东西,外层循环,整个数据结构。而 SIMD 类型处理的是我现在正手头处理的当前那一小块。
最后是数组。我们有各种不同的数组。有原生数组(raw arrays),有内存分配,有标准中表达内存连续性的各种方法,有很多这样的方法。有像 vector 和 array 这样的容器,还有像 span、string_view 和 mdspan 这样的视图。有连续迭代器(contiguous iterator),连续范围(contiguous range)。甚至还有 valarray。人们喜欢说 valarray 的坏话。但是确实有人在使用 valarray。所以每当我听到有人说它的坏话时我总觉得心里过意不去,因为这会让其他使用它的人感到难受。所以请大家对它温柔一点。好吧,虽然我确实不喜欢它的设计,但它依然是人们在使用的一个东西。mdspan 有访问器(accessors),这些也是表达连续性的方式。此外还有表达超对齐(over-alignment)的方式,比如 assume_aligned 或判断是否充分对齐。
我刚才列举了这么多关于数组的东西,听起来可能有点奇怪,因为你通常不会把数组和“并行”联系起来。但事实证明,它们与并行有很大的关系。它们是实现并行的赋能技术(enabling technology)。我将在这个演讲中对此进行详细阐述。
数组及其支持并行的三大特性¶
数组实际上揭示了三个支持并行的属性:不相交性(Disjointness)、局部性(Locality) 和 耦合(Coupling)。我会在接下来的演讲中逐一讨论。
当我们创建和使用数组时,我们会做出许多不同的假设。当你创建一个数组时,你假设如果我有不同的索引(indices),那么对应的元素也是不同的。这是一个合理的假设。对于单个数组的情况,我称之为不相交性:即不同的索引对应不同的元素。对于两个数组,标准向我们承诺,如果我通过默认方式创建两个对象,或者我创建两个内存分配,它们在内存中是不会重叠的。这里我创建了数组 x 和数组 y,我知道这两个数组在内存中不重叠。
当创建一个数组时,我也知道它具有某种形式的局部性。例如,我知道如果我有一个元素,我获取它的下一个索引(给索引加 1),那么下一个元素在内存中也是紧挨着的。这被称为空间局部性(spatial locality),或者说它启用了空间局部性。我还知道,如果我在一个数组中获取一个索引,我知道内存中的那个地址通常不会改变(当然,vector 可能会重新分配内存,前提是我们谈论的是同一次分配)。这个不变量启用了时间局部性(temporal locality)或重用。
然后是我称之为 耦合(coupling) 的东西。耦合主要关乎诸如对齐(alignment)以及内存与局部数据结构之间的关系。因此,知道数组在分配时是对齐的,并且能够分配超对齐的数据,这就是我所说的耦合。我将在演讲的后面详细说明。
如果你仔细观察代码,你可能会注意到这一点。我不仅仅是使用 < 和 <= 来比较我的指针,因为如果我使用标准的、内置的指针算术运算符,那只定义了一个偏序(partial order)。如果分配是不相关的,那么比较的值是未指定(unspecified)的。但事实证明,如果你使用标准的函数对象,比如 std::less 或 std::less_equal,你就可以比较指针,并且它们之上存在一个全序(total order)。这就是你如何决定指针是否重叠的方法。
不相交性与内存重叠¶
不相交性赋予了重新排序写入操作的权限。在并行方面,不相交性实际上关乎的是向数组写入数据。因此,对于单个数组,不相交性意味着每个元素都有一个唯一的索引,所以不同的索引指向不同的元素。这听起来微不足道,通常当你分配一个数组时,不同的索引自然代表不同的元素,对吧?
然而,描述数组有不同的方式。例如,在基本线性代数子程序(BLAS,C 和 FORTRAN 中用于线性代数操作的标准)中,有一种叫作……哦,那不是 Axby。那是另一个例程。那是,让我看看能不能想起来。GenV,抱歉,我拼错了。是 GEMV,这是 BLAS 中的基本矩阵-向量乘法。谢谢。它计算 y = beta * y + alpha * A * x,其中 y 是向量,beta 是标量,alpha 是标量,A 是矩阵。
BLAS 为向量接受整数类型的步幅(strides)。即内存中连续元素之间的距离(以字为单位)。你不能在这里填入 0。但是步幅为 0 在数学上是有意义的。它的作用是,如果我以只读方式访问一个数组,并且步幅为 0,那么数组中的每个索引都指向内存中的同一个地址,指向同一个元素。这是表示一个所有元素都相同的向量的绝佳方式,在数学上完全合理。例如,如果我想计算矩阵 A 各行的和,一个好方法是进行矩阵-向量乘法,其中向量全是 1。因此,我不需要在内存中存储 1, 1, 1, 1,我只需要在那里存一个 1,并使用一个步幅为 0 的数组。
但是 BLAS 不允许我这么做。但我们可以想象这样做。使用类似 mdspan 这样的工具,你就可以表达这一点。而且,对于用 Python 编程、在 Python 中做深度学习或其他数值计算的人来说,流行的格式(如 DLPack)通常支持数组的零步幅。如果我可以有零步幅,我就失去了不相交性。就像我之前说的,不相交性关乎重新排序写入操作的权限。对于两个数组,它关乎的是询问数组是否重叠。是否存在这样的索引:来自这个数组的一个索引和来自不同数组的另一个索引指向内存中的同一个元素?
如果我知道不相交性(两种不相交性都知道),那我就获得了对数组上的写入进行并行的完全权限。我可以向量化(vectorize),这意味着重新排序、为数据块(比如针对 SIMD 指令)进行优化;我也可以使用并行算法进行并行化。我在这里展示了 memcpy 和 memmove。它们之间的区别在于 memmove 允许数组重叠,而 memcpy 不允许。这意味着 memcpy 是一种并行算法,因为用户通过它赋予了算法以它喜欢的任何方式重新排序写入操作的权限。所以 memcpy 对于不相交的数组是没有问题的。你看,不相交性是非常有用的,即使你从未创建过线程或关心过多个线程。只要你调用 memcpy,你实际上就在关心并行。
我提到了关于不相交性的这两个问题。C++ 已经有一些方法来表示这些。关于我在分配数组时的语言规则:当我有不同的分配或创建不同的对象时,它们默认在内存中是不重叠的。理论上,你可以通过一种方式在 C++ 中声明你假设两个数组是不相交的。你现在可以使用 [[assume]] 属性声明来实现。但在实践中,我们看到的编译器并不能真正利用这一点。问题在于,是的,我做出了这样的代数陈述:我有这个指针和另一个指针以及一个长度,我承诺这个小于那个,并且这个小于那个。这太好了,这是代数。但编译器怎么知道我这是什么意思呢?它不知道。我只是说了一些数学表达式,编译器必须以某种方式对它进行模式匹配,弄清楚我想表达的是关于内存重叠的事情。这挺难的。
如果我能直接告诉编译器就好了。你可能知道 restrict 关键字。它不是 C++ 标准的特性。但是很多 C++ 程序员都在使用它,因为它是从 C 标准借用来的。而且它是 C++ 编译器中一个流行的语言扩展。所以编译器允许你使用它,而且声明它也有其特定的含义。通常当用户使用 restrict 时,他们希望它意味着类似:“我承诺这个指向一个数组,那个指向一个数组,并且这两个数组不重叠。”在编译器层面,它代表了关于指针表达式复杂得多的含义。C++ 标准一直抗拒采纳它,这有一段历史。Hal Finkel 等人在 2014 年提交的“别名集(alias sets)”提案谈到了这段历史。
但如果你能使用 restrict,编译器就能进行相当程度的优化。例如,在左边,我展示了一个反向复制(copy backwards)算法,我有一个输入数组 Q,一个输出数组 P 和一个长度。我要把数据从 Q 复制到 P,但由于某种原因我是反向复制的。如果在(且仅在)你放了 restrict 的情况下,至少我测试过的少数编译器,会从字面上把它变成 memcpy。你查看生成的代码,里面会有一个 memcpy 调用。
还有其他的例子,比如右边的代码,我有两个相互交错(interleaved)的数组 T 和 S。我从 S 复制每隔一个的元素到 T。我的 S 数组每隔一个元素跳过一次,而且 S 数组的起点在 T 数组起点之后的一个位置。所以 T 数组的元素是蓝色的,S 数组的元素是橙色的,它们互不重叠。元素之间不重叠,但它们在内存中共同占据了连续的块。如果我没有一种讨论重叠的方式,我就无法真正声明这两个东西不重叠。
实际上我想做的是划分(partition)。我想声明式地说:“我有一个数组,我把它分成不相交的块,而且我知道如果我写入这个块,我就不会写入那个块,反之亦然。”但通常当我拥有数组时,我能对它们做的操作是切片(slice),而不是划分(partition)。两者的区别在于,当我获取数组的切片时,我是在查看元素的子集,而我的子集视图与原始视图同时存在。现在我就有了两个视图,它们的元素重叠了。我制造了别名(aliasing)。但我并不想这样做。
我在左边展示了用 span 表达的快速排序(quicksort)的例子。我获取我的输入 span,对其进行分区,在中间选取一个枢轴(pivot)p,然后我说,好吧,左边的一切我取一个视图,右边的一切我取一个视图。现在我可以调用快速排序,我可以并行调用它,因为我知道这些东西是不相交的。然后我返回,我可以在此基础上递归。所以即使我知道 x_left 和 x_right 这两个分区不重叠,它们却与父视图重叠。所有这些生命周期都存在,它们同时存在。
我真的希望能进行分区,能有一种方式说:“好吧,这些原始父视图的子视图是不相交的,并且我承诺在使用子视图时,我永远不会去触碰父视图。”在 C++ 中,你可以通过创建鼓励这种程序员纪律的库来勉强表达这一点。比如,我展示了一个 apply_left_right 函数,它的作用是获取一个枢轴和某个视图,在中间将其分割,然后先执行左侧,再执行右侧。你不能“看到”我创建了子视图的这一事实。然后我可以用这种方式表达快速排序。在这个快速排序版本中,当我的父视图处于活动状态时,它不让我看到子视图的存在。但它们仍然存在。所以这是一种库层面的纪律。
在 Rust 中,你可以写这样的代码,借用检查器(borrow checker)会保护你。当我调用 Rust 中的 split_at_mut 进行分区时,我创建了这个分区 x_left 和 x_right,当 x_left 和 x_right 存在时,我不允许写入 x。语言本身阻止了我这样做。然后我就可以并行运行这个。如果我也能在 C++ 中这样表达就好了。
(此时有观众提问:Mark?你能回到上一张幻灯片吗?apply_left_right 里面的 template for 是在做什么?)
哦,是的。我之前没有解释 split_left_right。这只是某种假设的代码。但你可以想象 split_left_right 创建了一个视图 x_left 和一个视图 x_right,并返回包含这两个视图的元组(tuple)。因为它是元组,所以它具有编译时已知的大小,因此我可以遍历元组的元素。我在这里所做的就是为该分区的每个元素调用输入函数。但我是按顺序调用的。这也可以不用 template for,而用某种允许你这样做的并行构造。是的,谢谢你的问题。
总之,如果我能在 C++ 中表达分区就好了。这有点难。表达不相交性的另一种方法是,例如 Hal Finkel 领导的 2014 年的“别名集”提案。如果我们能教它认识标准库里的东西就好了,这样它就能知道 span 代表内存中的一块连续区域。但如果我有了分区特性,我其实就不需要这个提案了。如果语言免费给了我这种保证,我也不想去用 [[assume]] 之类的东西。但如果必须使用假设类的手段,我也会去用。关于我们到底有多需要这个,目前还是个问题。我们稍后会重新探讨这个话题。
空间局部性与时间局部性¶
现在,我将从不相交性转向空间局部性。空间局部性来自于连续性(contiguity)。连续性意味着,如果我在数组中获取下一个索引,它就是内存中的下一个元素。这实际上指的是虚拟内存,而不是物理内存。我未必能控制物理内存。但虚拟内存是物理内存的一种提示(hint)。你自然希望虚拟内存中的下一个地址最好也是物理内存中的下一个地址,从而获得相应的性能优势。
这对所有语言的程序员来说都非常重要。事实上,对于 C++ 用户来说,它重要到人们对最初的 C++98 提交的最早的一批缺陷报告(DR)之一就是:“我希望 std::vector 能够保证内存的连续性。”它在 C++98 中没有做到这一点。所以它作为针对 C++03 的 DR 被修复了。
计算机以各种方式针对连续访问进行优化。我在右边展示了 IBM 7030(即 Stretch)的例子,它是 1961 年在洛斯阿拉莫斯(Los Alamos)首次安装的。那是交叉内存(interleaved memory)的早期采用者,这意味着内存有不同的库(banks)。当你连续加载地址时,你会从一个库转到下一个库,再到下一个库,然后再循环回来。所以这是一种四个库一组的内存分块。你可以有多个库以获得更大的容量,但在某个限制内,它是一个四库交错系统。这种交错的目的是隐藏每个存储库的内存延迟。你可以开始第一次访问,在等待数据的时候,你可以开始下一次访问并继续下去。因此,交叉内存是隐藏延迟的一种方法,从而针对连续访问进行了优化。
你可能熟悉缓存(caches)的工作原理,知道有缓存行(cache lines)。缓存行的意义在于摊销(amortize)连续访问的延迟。如果你知道自己正在进行连续访问,你只需要为整个缓存行支付一次代价,而不是为每个字支付一次代价。我们的 GPU 有一种叫做合并访问(coalesced access)的技术,它针对不同线程共享连续块的情况进行了优化。它有点像缓存行。还有预取(prefetching)。我们的 GPU 支持异步 memcpy,一些其他系统也支持。内存页面(memory pages)以及转换后备缓冲器(TLB)的工作方式,即使说不上“优化”,至少也是偏好连续访问的。如果你真的在大范围的内存中跨步访问(stride),你会为此付出代价的。你很难找到故意不针对连续访问进行优化的计算机。我能想到的可能只有 Tera 以及后来的 Cray MTA 或 XMT,它们会随机散列你的内存访问,以避免过度频繁地访问同一个内存模块。那些架构是为了指针追逐(pointer chasing)、树和链表等结构而优化的。这些情况有点不寻常,它们确实存在,但大多数架构都是针对连续访问优化的。Stretch 是一个非常迷人的历史见证。它非常超前,非常有野心。尽管它在实现如“设定未来计算机的模样”这类直接目标上并没有成功,但它仍然是一项惊人的成就。
我已经谈论了空间局部性,针对连续访问的优化。时间局部性是针对重用(reuse)的优化,即时间上的局部性而不是空间上的局部性。我非常喜欢一篇关于“高等研究院(IAS)机器”(后来被称为 MANIAC)的论文中的这句话,该机器就像是战后计算机架构设计的模板。那篇论文甚至谈到了内存层次结构(memory hierarchies),因为你可能拥有容量巨大但速度很慢的内存。你的大内存很慢,你的快内存很小,你不可能鱼与熊掌兼得——既快又大。这是不存在的。即使在那个年代,你也不可能得到既快又大的内存。
我记得我导师的导师曾谈论过在磁芯内存(core memory)存在之前的核外算法(out-of-core algorithms),所以称其为“out-of-core”有点时间错乱。他有一个转鼓(rotating drum)作为慢速内存,而快速内存简直就是——我记得是选择管(selectron),看起来有点像阴极射线管的东西,它有一个监视器,你可以看到亮点,那些亮点就是你内存里的比特位。那相当疯狂。但是转鼓的问题在于,哦,谢谢,它会不停地转。所以如果你错过了你的数据位,你错过了你的位置,你就必须等待这个机械装置转回来。所以你需要承受这种延迟。你当然不希望每次访问都要承受这种延迟。
我展示了一些参与了 IAS 机器的人员。我觉得这在普林斯顿 IAS 机器前拍的合影非常经典。Alan(可能指台下的某位观众),你可能认出右数第二个人,他最近还在一部电影中被描绘过。但我喜欢这张照片是因为约翰·冯·诺依曼(John von Neumann)在笑。人们通常喜欢展示他的洛斯阿拉莫斯证件照,他在那里看起来很暴躁。但我也会这样,因为当你刚开始在实验室工作时,你必须去办一堆文书工作、到处跑,我能想象他肯定很累。而在这里他处于他最如鱼得水的环境中,所以我很高兴。
内存重用关乎针对内存层次结构进行优化。内存层次结构在各个层级的并行处理中都很重要,并且它倾向于奖励一种块级重用模式(chunk-wise reuse pattern)。我想充分利用我那虽小但速度快的内存,因此我将我的局部数据块调整到适合放入我那小而快的内存中的大小。我把它加载进来,尽可能多地用它进行计算,当不能用它做更多事情时,我就把它放回去,为下一个块腾出空间,从而充分榨取小而快内存的价值。
人们提到向量计算机时,往往不怎么联想到它,他们往往将其与基于缓存的计算机进行对比。人们会说,这是向量机的做法,那是缓存机的做法。然而,即使是 Cray-1 也有内存层次结构,并且它同样针对时间局部性以及空间局部性进行了优化。例如,Cray-1 有一个巨大的用于向量数据的寄存器文件,如果你足够聪明,如果指令的延迟计算得当,你可以进行一种叫做“链接(chaining)”的操作。它的意思是,你从内存中加载向量数据,将其放入向量寄存器中,然后你可以跨多个算术操作或多个向量指令重用它。这就是(Cray-1)配备向量寄存器文件(而不是仅仅从内存中加载东西,进行算术运算,然后再放回内存)的理由之一,链接就是为了这个。
所以 Cray-1 实际上有一个内存层次结构。这是从早期的内存-内存向量机器(如 CDC Star 100)中吸取的教训。Cray-1 也有所谓的即时寄存器文件(immediate register files),比如 B 和 T 寄存器文件,用于标量和地址。将数据从它们传输到主寄存器的延迟非常低。它们还支持与内存之间异步的块传输。看到这些想法在历史中一次又一次地出现,真是有点疯狂。但这恰恰说明了,时间局部性就是一切。
我一直在谈论局部性和重用。这其实是个伏笔——我为什么要讲这些?我完全还没提到“多维”这个词。事实证明,如果我执行一个一维算法,我只是在数据上遍历一次,我能获得的重用量受限于我对数组或范围内的每个元素执行了多少工作。我在这里展示了一个 ranges::transform 的例子。如果我只是做一个加法和一个乘法,那工作量并不大。这意味着我算法的成本受到带宽的限制,也就是我的内存有多快。而这远比我的计算机算术运算速度要慢得多。
相比之下,有很大一类多维算法中存在着更多的潜在重用和可实现重用,即使我每个元素的计算量并不大。我在右侧展示了一个稠密矩阵的矩阵乘法,但还有其他一些算法具有类似的重用属性。
耦合与 SIMD¶
接下来我要谈谈耦合(coupling)。我谈过不相交性,也谈过空间和时间局部性。这些主要都是抽象的概念。你可以编写算法来利用它们。你不一定非得使用(尽管有时会有所帮助)非标准、非常特定于系统的库或指令。你可以编写尊重不相交性、空间和时间局部性的可移植的 C++ 代码。这是 C++ 的惊人之处:你真的可以在其中表达所有这些东西。而且你甚至不需要在代码里写死依赖于硬件的常量。
但并行化的一些特性是与内存耦合在一起的。当我说耦合时,我的意思是什么?耦合的一些表现是对访问的限制。比如我有一个东西,它感觉就像应该是一个数组,但它实际上并不是,因为我无法获得对单个元素的引用(reference)。超对齐(Overalignment)是耦合的另一个症状。当访问内存时,我只能以某种特定的超对齐方式进行,或者这是被强烈推荐的方式,如果不这样做就会很慢。此外,内存访问的大小与并行操作的规模之间也存在耦合。还有就是局部内存(local memories)。有时我能看到它们并对它们编程,有时它们在编译器的控制之下,我看不到。
耦合的影响体现在取决于硬件的选择上。具体来说,我希望这些依赖硬件的选择只是像整数常量(大小、维度)之类的东西。但有时候,很不幸,我不得不采用算法的不同变体,也就是完全不同的算法。
我之前谈到的耦合的一个例子是对字(words)进行按位操作。那些是并行算法,是你字数据上的数据并行算法。谈到访问限制,对吧?我没法单独访问一个位。我没法直接更改一个位。我可以更改整个字,我可以通过掩码把它屏蔽掉,让它看起来像是我更改了一个位。我可以读取一个位,我可以写入整个字,以这种方式实现对某个位的写入,但我真正在访问的是整个字。我无法直接写入那单独的位。我也无法观察位的存储顺序。我知道有位,我可以假设存在一个概念上或算术上的存储顺序,但我并不知道物理上的顺序究竟是什么。而且当我对比特位进行内存访问时,它们是按字符(characters)对齐的,是按 8 位组对齐的。实际上它们的对齐要求比这个更高,但抽象上是字符对齐、8 位对齐。这又与字上可用的并行性相耦合。所以我对字进行位操作,而不是对位。字与内存是相关的。我将字从内存中加载和存储。我在字上工作并放回去。我不打算……无法访问字的部分原因在于那里有一种局部内存,也就是我电脑里的寄存器。我把字加载到寄存器中,在寄存器上操作,然后再把它们放回去。寄存器就是一种局部内存。就像我们刚才谈到的 Cray-1 拥有一个寄存器文件一样。
但我通常不需要对着寄存器编程,除非我在可移植的 C++ 之外编写汇编语言。编译器负责处理寄存器,我处理的是值(values),我对值编程,编译器会尽力将这些值放入寄存器中。
其中也涉及到一些常量。也许在你的电脑上,32 位字的操作更快。所以在进行并行的按位操作时,我应该使用 32 位的字。也许在另一台电脑上是 64 位的字更快。这些就是常量,我必须了解硬件。如果我在编写完全可移植的代码,我可能会将其作为模板参数或函数的编译时参数。我也许想要自动调优(auto-tune)或测试它,或者准备好不同的变体。在某些情况下,如果你真的是……我不是那种超级痴迷于位级并行的书呆子。我平常不干这个。但有些人真的很喜欢使用位级并行,他们对“位”了如指掌,他们能顿悟位的奥秘。他们知道在某些硬件上你应该这么做,在另一些硬件上你应该那么做。我所说的这么做、那么做,指的是不同的算法,是为实现相同目标而设计的不同的并行算法。尽管它们看起来都是加号、乘号、与(and)和或(or)。它们依然是不同的算法。因此,它们是算法的变体,我可能需要在不同的硬件上以不同的方式使用它们以获得更好的性能。
我又有一个来自 Cray-1 的耦合示例。我刚才谈到了链接(chaining)。好。我谈到了 IBM Stretch 的交叉内存,Cray 也有交叉内存,一般在正常安装下,你会有以 16 为一组的存储库(banks)。当你连续加载时,每个存储库提供一个字。并且你可以每 4 个时钟周期访问同一个存储库一次。因此,如果你以 8 的倍数为步幅进行加载或存储,你实际上会遭受带宽惩罚(bandwidth penalty)。换句话说,每个字都将招致延迟惩罚,但对于你整个数组或向量来说,你会遭受带宽惩罚。如果你的步幅是 8,你的带宽惩罚是两倍,即 50% 效率。因为你不再能够完全重叠隐藏存储库的延迟,而是要遭受 2 个周期的延迟。
我谈到过的链接机制所在的局部内存,通过它你可以……你实际上可以在一个时钟周期内产生多条指令的结果,方法是在前一条指令还在执行中间时就开始启动下一条指令。因此,通过使用这种局部内存,即用于向量的寄存器文件,你可以每个时钟周期产生超过一个结果。但是,内存和链接之间存在耦合。因为如果你有一个步幅是 8 的倍数,你就无法进行链接。在 Cray-1 上你会错过链接槽(chain slot),因此你将无法实现这种指令的重叠。当时的技巧是填充(pad)你的二维数组。你会使得步幅与存储库的数量互质(relatively prime)。这也是一种整数常量,对吧?你必须了解该硬件的一个填充因子(padding factor)。
有趣的是,这实际上对我们现在的 GPU 也很有用,因为 GPU 上的小型共享内存单元(shared memory units)也有存储库的分块(banks)。人们通常使用的技巧之一就是对数组进行填充,以帮助避免这些 bank 冲突(bank conflicts)。我也曾在像 Sun Niagara 这样缓存关联度非常低的架构上看到过类似的情况。有时你会遇到类似 bank 冲突的效应,导致大量的缓存未命中(cache misses)。所以填充数组也能帮助解决那个问题。这也是一种耦合。
现在,回到看起来像按元素(element-wise)的数据并行操作,我会谈谈 SIMD。在使用 SIMD 时,你确实会遇到访问限制。你的 simd_vec,你的 simd_mask,它不是一个数组。你无法获取对某个元素的引用(reference)。这不是一个合法操作。你无法获取指向某个元素的指针。虽然有一个方括号操作符 [],看起来像是在访问数组,但那只是一层皮。你得到的是一个值,而不是引用。没有 .data() 成员。有一个迭代器,但它不是连续迭代器。
之所以这样做是有优化层面的原因的。因为一旦我获取了一个指向这个 SIMD 结构中间某个东西的指针,那个 SIMD 结构就必须存在于内存中。因为我需要有一个可以引用的地址,我需要能够写入它,读取它。于是它在内存中就变成了一个数组。我们称这种现象为物化(materialization)。物化意味着编译器再也不能把这个东西当作一个纯粹的值来推理了。它被强制要求在内存中腾出空间并将其放在那里。而物化意味着你不能直接把它留在 SIMD 寄存器中。现在它占据了内存,占用了你缓存中的空间。现在你必须从内存中将它加载到向量寄存器中来操作它。因为很多时候,SIMD 指令只对生存在 SIMD 寄存器中的东西起作用。这里 SIMD 的局部内存就是寄存器。
这其中也有大小限制。标准对 SIMD 有一个最大尺寸限制。这是一个由实现定义的数字,比如 64 是它的一种最小保证。所以我不可能拥有大小为 1000 的 SIMD 类型。大概率不可能。而且我可能甚至也不想要大小为 1000 的 SIMD,因为那不是 2 的幂。很多 SIMD 指令喜欢 2 的幂,并会针对这些数字进行优化。
我也提到了局部内存作为耦合的另一个症状。用户声明 SIMD 变量,这就是你使用 std::simd 库的方式,编译器会尽力将它们放入 SIMD 寄存器中。如果我强制物化,那么它就做不到。但假设我不这样做。编译器就会将那些 SIMD 寄存器作为输入和输出用于 SIMD 指令。所以它们就像小型的局部内存。就像 Cray-1 拥有一个寄存器文件一样,你的拥有 SIMD 指令的计算机也有一个用于 SIMD 值的寄存器文件。指令必须使用这个局部内存。
SIMD 操作还偏好或要求以超对齐(over-aligned)的方式访问内存。通常需要对齐到 SIMD 的宽度。虽然并非在所有架构上都必须如此,但如果它是超对齐的,至少会运行得更快。那个 SIMD 宽度就是你需要对其编程的整数常量。我在右侧展示了一个例子。有一个 simd_vec 的 float,我有一组 SIMD 宽度的 float 数据,我声明了超对齐。我告诉系统,“嘿,SIMD,这个东西是以 SIMD 宽度 × alignof(float) 来超对齐的。”然后当我进行部分加载(partial load)时,意思是“从这个地址加载”,并承诺该地址是超对齐的。那里有一个小惯用法,涉及一个叫做 iota 的东西(那个希腊字母)。iota 是 0, 1, 2,以此类推。对于 SIMD iota,这是直到向量大小。所以它是 0, 1, 2, 3,等等,直到 vector_size - 1。这是用于掩蔽(masking)越界元素的一种常见 SIMD 模式。我们在 C++26 中实际上并没有直接提供这个功能,但你可以使用生成器构造函数(generator constructor)自己实现它。
SIMD 本质上是一个用于创造局部性机会的系统。当你手里有一些分散的杂乱数据(hogs,双关语)时,你可以通过把它们变成 span(同音 spam,午餐肉)来增加它们的空间局部性。对不起,我觉得这个笑话很有意思。所以我谈到了这种模式:我遍历我整个全局数据结构,我加载进一小块,在它们上面进行操作,榨取出局部性的好处,然后再把它们放回去。SIMD 自带了这类功能,例如,SIMD 自带了收集(gathers)和分散(scatters)。收集操作将我散落在内存中、不连续的数据拉取进来,把它们变成一个连续的块。这是一个非常重要的操作,以至于有些计算机为它配备了硬件加速。我还谈到了带掩码的加载和存储(masked loads and stores)。它们的作用是,我只加载某些掩码为 true 的元素,在掩码为 false 的地方,我只是用默认构造的东西(在这个例子中是零)把它们填满。收集和分散是处理内存中非连续的数据。掩码则是处理非连续的控制流。比如,使数据访问变得非连续的 if 语句。这有点像 copy_if 算法,由于条件语句的存在,输入可能是不连续的,但输出会被打包成连续的形式。这实际上引出了一点,即我在所有这些小的 SIMD 块上全局地并行循环。然后我在这些块上以局部的数据并行方式工作。
如果我现在想这样做,我有两个选择:我可以自己手写(roll my own),或者我可以使用 SIMD 算法提案,我稍后会谈到它。但如果我必须自己编写全局算法,首先我必须断言超对齐。所以我可以使用断言超对齐的各种方法(assume_aligned)以及用于超对齐的 SIMD 标志。然后我必须在一个索引范围上执行并行算法。在这种情况下,我以 SIMD 宽度的索引组跨步前进(striding)。而且我的范围是一个索引范围,所以我必须做索引算术。我必须创建一个掩码以确保我不越界。我还必须做加载、存储和向量相关的事情。
Matthias Kretz 已经提出了将 SIMD 与并行算法集成的提案。它的工作方式是通过基于执行策略(execution policy)的重载。将执行策略更改为 SIMD 执行策略(SIMD execution policy),会改变算法的行为。它改变了算法处理的元素类型,使得它看起来像一个 std::simd_vec。这非常棒,因为它处理了我刚才在示例中展示的你必须自己做的所有事情,对吧?它处理了超对齐,确保了访问的不相交性,并处理了所有的循环索引。但如果你想使用这个 SIMD 算法,你必须找到一个能适应这些算法模式之一的算法。而且,标准中的许多 range 适配器(adapters)并不满足连续范围(contiguous range)的要求。所以你可能从 vector 或 array 等在内存中连续的东西开始。但一旦我开始对它应用 range 适配器,我就失去了连续范围特性。然后我就不能对那个东西进行 SIMD 的加载和存储了。
所以关于这个提案还有一些悬而未决的问题。它需要针对 parallel ranges 算法进行更新。而且它只讨论了顺序 SIMD 迭代。所以还需要为如何集成进行设计……
(Ruslan 举手提问)
Ruslan: 如果我可以在这里快速评论一下,我脑子里正在思考,把并行算法(包括 parallel range algorithms)耦合起来,并将 SIMD 作为策略(policy),这是否是个好主意。也许它应该是一组独立的算法。但这已经是一个很老的论文了,它肯定需要被修改。
Mark: 是的,我不确定人们是否需要这个……他们肯定想要某种子集。
Ruslan: 这几乎就是我们在开始着手之前需要考虑的所有 80 多个问题。
Mark: 确实。是的。
是的,我会……我其实正要谈到这一点。存在各种表达代码中并行性的方式,有一个连续的图谱。正如 Ruslan 提出的,SIMD 算法在这个图谱中处于什么位置是个问题。这个连续体是关于我所说的“隐式(implicit)”与“显式(explicit)”的对比。
如果我在算法中是完全隐式的,我只说,我有一个数组,我把它加到另一个数组上。这天然就是并行的,这就是并行的。在完全隐式的模式中,系统会全面解决到底由谁将工作项(work items)映射到执行代理(execution agents)的所有细节。而且系统会说,我只声明工作项之间在局部上的交互。也就是说,通过将两个数组相加并声明那是并行的,我只是在说,逐元素的对应项对齐了。我不需要谈论事物之间的全局交互。使用隐式并行性,是由系统来确保局部的不相交性,即确保我不会以一种糟糕的并发方式去写入事物。但尽管如此,总是需要由用户来承诺我的数组是不相交的。所以在基于数组的编程模型之间实际上有一个连续统一体。SIMD 就是某种隐式并行机制(在基于数组意义上)的一个例子。
你可以想象 valarray。valarray 曾经有过并行的实现。人们谈论过……我认识 Alfredo,他可能在观众席中……哦,你好。你谈到过……你写过一个并行数组库,对吧?是的。那就是一个例子。或者在其他编程语言中,那是表达并行的常规方式。而处于连续体最远端的“显式”并行版本,就是我仅仅去手动创建线程,由我自己处理数据到线程的映射。但在这两者之间有一些很有趣的东西。
标准并行算法具有“按元素并行(element-wise parallelism)”的特性。我不对整个数据结构全局地工作,而是逐个元素地工作,除了当前元素,我不必考虑其他元素,而且我也不被允许去访问它们。Matthias 的 SIMD 算法提案就是这种例子。
基于分区(Partition-based)的并行,可能不是我们在标准库里所熟知的东西。我之前展示过 Rust 中 split_at_mut 和 join 的例子,其中由我自己来创建分区。我声明,“嗯,我想要那边的那一块和这边的这一块”,然后由我决定如何进行分区。所以我是在自己处理工作项,或者自己处理负载均衡。做分布式、大规模并行计算的人通常用分区的思维思考,这对那些人来说更自然。
std::execution 中的块组块(bulk chunked)确实创建了小的分区。确切地说,是系统创建了分区然后交给你。当我提到 for_each_index 时,它引入了一种间接层。你不必在工作项(work items,如同标准算法那样)上进行并行,而是在进入那些工作项的 索引(indices) 上进行并行。然后你可以决定如何将这些索引映射到工作项。而且你可以访问当前元素之外的东西。例如,如果你想进行最近邻计算,或者你有时候由于某种原因想要读取那边的数据。所以你要靠自己来确保那样做是安全的。
正如我所说,Matias 提议的 SIMD 算法是“处理每个元素(for each element)”风格的算法。但对于那些以编写 SIMD 代码为生的人来说,这其实并不太自然。以编写 SIMD 为生的人(尽管这有点过度概括的风险),他们其实不相信库。他们几乎完全不信任库。他们想自己触摸那些比特位。他们想自己来做。他们想自己写索引。而你现在的提案阻止了他们这样做。他们不确定能否获得他们想要的所有性能。
但是,我们在过去两个 C++ 标准中取得的许多进展,让编写“对每个元素(for each element)”类型的并行算法变得更加容易。因此,我认为我们不应该害怕去尝试它们。即使我们最终选择了 for_each_index 的设计,我们也不应该害怕去探索划分数据和表达并行性的不同方式。
为什么要走向多维?¶
现在我将插入一段关于多维性的讨论。支持多维算法和多维性表达的一大动力是:一些多维算法,很大一类多维算法,具有更多潜在的和实际可实现的重用,即时间上的重用,时间局部性。如果你只从一维算法开始,是很难看出这一点的。很多一维(扁平范围的)算法,尤其是那些在数据上只扫一次的算法,它们访问 N 个元素(读、写),并且对每个元素做 1 个单位的工作。所以总工作量也是 O(N) 的。那么潜在的重用率是多少?那就是,我有多少工作要做,我必须读写多少数据?因为最好的情况是我获得了完美的重用、完美的时间局部性,那么我可能获得的最大重用量上限,就是工作量除以数据量(work / data)。对于许多标准算法来说,这个比值只是 1。做了 N 个工作,读了 N 个数据,重用率为 1。
这种算法的一个例子是向量的点积(dot product),对于标准算法来说,那对应的是二元 transform_reduce。我遍历两个范围,逐个元素地相乘,然后把结果加起来。我们可以将其与增加维度后的情况进行对比。
我在右侧展示了矩阵-矩阵乘法。如果你熟悉求和公式,我是将 A 的 (i, k) 元素与 B 的 (k, j) 元素相乘后求和(注:演讲者口误,讲成了 i,j 和 k,j,但意思传达到了)。A 是二维的,B 是二维的,但这里有三个循环维度。最后的结果也是二维的。因为其中一个维度被加和了,被归约(reduced)了。因为有三个维度,每个维度大小都是 N,我真的有 N^3 的工作要做。当然有其他制定矩阵乘法的方法,它们也具有良好的局部性属性。我就不提 Strassen 或 Winograd 或其他那些算法了,我们只限制在我们探讨的这种 O(N^3) 复杂度的矩阵乘法上。所以这里有 N^3 的工作量,但只有 3 * N^2 的数据量。因此,我潜在的重用率是 N^3 除以 N^2,约等于 N。这比点积的重用率要大得多。而且它随着数据的增加而增长。我的问题规模越大,我能获得的重用就越多。
但并不是所有的多维算法都是这样的。在中间我展示了矩阵-向量乘法。矩阵-向量乘法就像我试图用图表展示的那样,它有点像是矩阵乘法的一种特例,其中一个对象是 1 阶(秩为 1)而不是 2 阶的。就像对于矩阵向量乘法,即使我在输入向量中获得了一些重用——我可以把它想象成 x 与 A 第一行的点积,x 与 A 第二行的点积,所以我在 x 上获得了重用——但它仍然只有 O(1) 的潜在重用。我有 N^2 的数据,做了 N^2 的工作。前面我谈到了 Cray-1 上的链接(chaining)。像矩阵向量乘法这样的操作,它们是系统擅长处理的,因为有时间局部性可以利用,只不过那不是算法层面的时间局部性(algorithmic temporal locality)。
有很多这样的算法。你可能会想,“哦,矩阵,线性代数,那是很特殊的东西。”其实有很多算法都是类似的。事实证明,许多看起来像三个嵌套循环的算法都是这样。即使我有更多的嵌套循环,比如张量收缩(tensor contractions),它们本质上也就是矩阵乘法,或者说我进行了一系列矩阵乘法运算。我可以重新排列索引把它们组合起来,它们本质上依然是矩阵乘法。很多东西都是矩阵乘法。
我之前谈论的是潜在的重用,但实际的重用受限于我局部内存的大小。这就是为什么你要将算法设计为装满你的局部内存,读入一个数据块,尽可能多地用它进行计算,然后再放回去。你受限于你那容量有限的高速内存、局部内存或缓存。我们所谓的 3D 嵌套循环算法很多都表现出这种行为。有一个理论上的局部性下界,而且很多算法实际上能够达到这个下界。实现这一目标的方法就是将你的问题打破成一块一块的“瓦片(tiles)”或分块(chunks)。有一篇文章提供了相关的理论基础和几何直觉,供大家参考。
C++ 中多维特性的现状与未来¶
现在我将谈谈 C++ 对多维性的支持:我们现在有什么,以及我们想要什么。我会从非多维的情况开始逐步构建。
我之前提到过,C++ 目前以四种方式暴露一维并行性。有异步执行,有并行算法,有数据并行类型和操作,还有数组。我从数组开始。在 C++ 标准中我们已经有了很多数组相关的设施。在 C++23 中有了多维数组视图 mdspan。在 C++26 中有了切片(slicing),所以我们可以表达访问 mdspan 的不同方式。但我们还没有多维容器(multidimensional containers)。谢谢提醒。我们还没有容器或者进行多维分配的方法。有一个 mdarray(多维数组)提案,这是一个好的开始,但我认为它还需要一些工作。我待会儿会谈到它。
我把异步执行和并行算法放在一起看。是的,你可以对一个多维数组进行排序(sort),这也可以有确切的意义。但这更像是在一个扁平数据结构上做的事情。所以这并不意味着我们需要跑到标准委员会说:“我们必须拥有多维 sort 算法,老兄,我们必须拥有它!”不,我们不需要那个。所以有很多工作我们并不需要做,这是件好事。
但是,随着多维度的到来,也伴随着新类型的算法。特别是,你有了关于维度(extents)子集的选择。例如,我可能有一个矩阵,我想计算行和或者列和。我不想把所有的元素都加起来。我想从一个秩为 2(rank 2)的东西变成一个秩为 1(rank 1)或秩为 0(rank 0)的东西。或者我从秩为 3 降到秩 2 或秩 1。我将这些称为降秩归约(rank reducing reductions)。还有张量收缩(tensor contractions),它是矩阵-矩阵乘法的推广形式。这就像多维版本的二元 transform_reduce。所以降秩归约就像是一元的 transform_reduce。张量收缩就像二元的 transform_reduce。
还有其他一些更奇怪的算法。有可能在扩展范围的子集上进行扫描(scans)。也许你混合了扫描、归约和 for_each 操作。也许你有非规则的多维数据,比如它不仅仅是不能用 mdspan 表达,而是它是参差不齐的(jagged),或者是分段扫描等类型的操作。这里有太多事情可以做。不要试图一次性解决所有问题。不要提出大包大揽的混合法案式(omnibus)提案。我之前还讨论过常规的 SIMD。那多维 SIMD 是什么?那到底是个啥?
说回 mdspan,mdspan 是 C++23 引入的多维数组视图。它推广了 span 的概念,包含三个部分:
数据句柄(data handle):推广了“指针”的概念。这是一种非唯一拥有的东西,通常是不拥有的,但也可以是共享指针。
布局映射(layout mapping):将你的多维索引映射到一个一维的偏移量。
访问器(accessor):接受数据句柄和来自布局映射的偏移量,然后返回一个对元素的引用。
我在右边黑色方框里展示的是在 mdspan 中数组访问操作符看起来的样子。字面意思就是这样:我获得访问器;我获得数据句柄;我根据我的多维索引评估布局映射,得到一个偏移量,传入访问器的 access 成员函数。出来就是一个引用。
我之前谈论过不相交性和局部性。mdspan 的布局表达了这些属性。它们表达了不相交性和空间局部性。这就是布局映射。你可以向布局映射查询这些属性。布局映射自带了查询接口:is_X 或 is_always_X。is_always_X 是一个编译时的问题,而 is_X 可能是一个运行时的问题。你可以询问一个布局是否是唯一(unique)的。这意味着它是单射的(injective),一对一的。所以每个索引都映射到一个唯一的元素。这立刻就代表了不相交性。
你可以询问布局映射是否是穷尽的(exhaustive)。这意味着它是满射的(surjective),它覆盖了整个一维范围。但它不一定是唯一的。右侧展示的一个穷尽布局的例子是唯一的,但它不必如此,可能会有重复。你还可以询问布局映射是否是跨步的(strided)。跨步意味着它是索引的线性函数。因此,我可以计算索引与步幅(整数)的点积来评估跨步布局。所以这些特性映射到我们之前谈论的数组的并行属性。不相交性和唯一性,意味着同样的东西。连续(Contiguous)意味着:唯一 + 穷尽 + 跨步,这三个加在一起就给了我连续性。
不过这里有个需要注意的限制(caveat),因为在 C++ 标准中,连续性也暗示了存在内存和一个指针。但使用 mdspan 访问器,我并不必须拥有内存或指针。也许你们中有些人看了那场非常精彩的 CppCon 演讲,演讲者使用 mdspan 来表达机器人手臂伸向远方并从小盒子里抓取啤酒瓶的动作。那里根本没有传统意义的内存。我也许可以说那里存在“物理内存”,只不过那里装的是一箱啤酒瓶。所以它不一定非得是内存。连续暗示了内存。所以我在这里是比较宽松地使用这个词的。
有时候我想问,是否部分维度(extents)是连续的。这对于实现(比如)多维数组复制可能会很有用。你可以通过确认特定的布局是否具有该属性来询问。通过知道该布局是跨步的,并且其中一个步幅是 1 来确认。
布局映射具有这些查询接口。你可以问所有这些问题,这非常棒。但是访问器(accessors)却没有这些查询接口。我的意思是,我不仅不知道它是否在内存中,在用访问器时我甚至不知道不同的索引是否指向不同的元素,因为标准里没有任何规定说这必须是真的。虽然那样会很奇怪。我想说,惯用的(idiomatic)mdspan 访问器确实具有唯一性属性。它们具有不相交性。因为你可以直接在布局映射中表达非唯一性。mdspan 的设计非常刻意地将布局映射和访问器解耦了。非常刻意,是基于他们以往设计的经验。所以如果在访问器中表达非唯一性,我觉得会很奇怪,你根本不需要那么做,只会让事情变复杂。但标准并没有强制你编写符合惯用的访问器。所以也许我们应该考虑为访问器添加查询接口。或者我们甚至应该考虑直接规定访问器拥有唯一的元素。我倾向于后者而不是添加查询接口,因为我能拿一个非唯一的访问器干什么呢?那太怪异了。
对于一维的情况,我谈到了切片(slicing)和分区(partitioning)。在一维中,我甚至无法区分这两件事(如果你不涉及生命周期的话)。在没有 restrict 或类似的声明时,编译器很难推断出交错的用户数据是不相交的。我展示了那两个交错数组步进访问的例子。如果没有 restrict 关键字,编译器很难处理这个问题。事实证明,当你编写多维算法时,你更有可能创建这种交错的布局。例如,矩阵分解的一个经典模式是原地(in-place)分解,方法是将其划分为一个 2x2 的块矩阵。得州大学奥斯汀分校的 Flame 项目针对这种算法有一系列非常漂亮的出版物,使用了非常优美的计算机代数系统。通常那是一种表达最优矩阵分解的惯用方式。而且当你这样做时,你选择什么布局并不重要。我在右侧展示了一个列主序(column major)布局(即 layout_left),其中一列中连续的索引在内存中是连续的。但如果是行主序布局你也会遇到同样的问题。你可以编写一个布局,使每个分区在内存中都是连续的。但对用户来说,那真的不自然。所以在这里你遇到交错布局(interleaved layout)的问题了。
我在这里展示了一个带部分主元选择的 LU 分解(LU factorization with partial pivoting)的例子。这只是说明,这种分区是一种表达我们想要表达的那类算法的自然方式。所以在这里我展示了这个分区操作。因为在 C++ 中我没有原生的分区语法,所以我必须在这里使用切片(slices),所以我正在使用 submdspan。但这就是你现在在 C++ 中编写矩阵分解的方式。
而我想要分区功能。我不想用这些与它们的父对象同时存活的切片。我想要某种方式来表达这些部分是不相交的。我还不确定那到底该是什么样子。但我认为我们需要更深入的理解:这真的有多重要?我想向标准添加新东西来表达这些切片的不相交性吗?如果我向标准添加新东西,它需要能理解 mdspan 和 span。我不能只是让它处理裸指针,因为正如右侧示例所示,一个指针和一个大小,并不足以表达我们想要表达的那种不相交性,尤其是对于多维算法而言。因此,任何声明语法都需要理解 mdspan。
我之前谈论过 mdarray 提案。那是个可爱的提案,因为你可以避免“由于我需要声明一个布局映射,以便知道我需要多少存储空间,才能去创建一个 mdspan”这种窘境。目前所有的例子都显示这感觉很蠢,为什么我不能直接拥有一个能做到这一点的容器呢?而这个提案之所以有点停滞,原因之一是,如果我对一个 mdarray 移动(move)操作,会发生什么?
mdarray 并不是像 vector 这样的独立容器。它是一个容器适配器(container adapter)。它包装了你的底层连续容器。而那可能是像 vector 一样动态分配的东西,也可能是 std::array。所以,移动后的行为取决于它表现得像 vector 还是像 array,但更糟糕的是,在标准中,连续容器并没有一个统一指定的“被移动后(move from)”行为。假设我有一个包含一些静态边界和一些动态边界的 mdarray,或者它有纯静态的边界。这就像是最糟糕的情况。然后我对其执行 move 语义。我希望那些维度大小(extents)在 move 之后变为零。但它们不可能变为零,因为它们是编译时的常量!我无法将编译时的数值改为零。因此,如果我有一个纯静态维度的 mdarray,它应该总是表现得像 std::array。std::array 的维度在被移动后是不会改变的。你只是移动了里面的元素。
由于这种对“被移动后”行为缺乏明确规范,使得在泛型代码中推理 mdarray 变得非常困难。将 mdarray 当作某种分配器来使用是一种常见模式。你创建一个 mdarray,然后对其进行切片分区并在其上做事情。然后你对其各个块调用函数,之后你将其移动到流水线的下一个环节。这样你只需做一次分配,然后你可以在不同的方式中重用它。但如果我不知道被移动后的行为是什么,我就写不出这种代码。
mdarray 的另一个问题是它混淆了小规模和大规模操作(small and large)。当我创建一个 std::array 时,那通常是比较小的东西。我无法分配很大的因为那样会耗尽栈空间。我在那里有编译时常量。这是一个局部的、小型的数据结构,作用于局部且小型的物体。而“大规模”更像是 std::vector。我需要一个分配器(allocator)。它是动态分配的。我希望 move 是移动底层存储资源,而不是逐个复制。构造一个大 vector 并将一个 vector 的内容复制到另一个 vector,这是全局的操作,这是算法。这是并行算法。因为我有很多元素,我把它们复制到那边,我都不知道具体有多少。我也许想并行运行那个复制。
目前的 mdarray 设计有一个接收 mdspan 的构造函数。那里面隐藏了一个可能极大的并行算法操作。因为数据可能很多,而且那里面隐藏的并行算法并不关心我的数据最终想要生存在哪里。我的 mdarray 可能会使用一个自定义的分配器。内存可能被分配在 GPU 上。它可能被分配在一个我需要快速将数据复制过去的地方。而这一切都被隐藏在这个构造函数里,我无法控制它。这很难推理,也很难编写这种代码。也许我甚至希望异步地运行它。但是容器并不是异步的。这感觉很奇怪,对吧?容器通常给人一种非常同步的感觉。你创建一个容器,然后你就能立即访问它的元素。到底什么是“异步的容器”?比如在我们的 Rapids RAP 库中,我们部署了一个带有一些异步 CUDA 流的 mdarray 版本。复制是异步的,并且它们沿着 CUDA 流被同步。那确实能工作,但把这种行为硬塞进容器语义里真的非常尴尬。
但有趣的是,用户偏偏想要那个看起来像 std::vector 的东西。他们就是想要它。即使在架构上它是错的,他们依然想要。所以这里存在一种斗争。我们想给用户他们想要的东西,但他们想要的东西是不恰当的。但这毕竟是“顾客永远是对的”的情况,你要怎么办?我不知道。
我个人希望 mdarray 成为一个真正的容器,而不是一个适配器。最初将其设计为容器适配器的部分原因是编写标准措辞图一时方便(expediency of wording)。但结果发现那反而让措辞变得更加困难,这就很搞笑了。我认为它应该是两种不同的东西。应该有一种全静态维度的容器,表现得像 std::array。然后还应该有动态分配的版本,表现得像 vector(也许可调大小,也许不可调大小,我不确定)。而且我永远都不想要一个动态分配的多维数组却非要表现得像一个简单的容器。我永远都不想要那个东西。所以我不会去写那个提案。如果你们中有人想要的话,你们可以去写那个提案。如果有人想要那个,我也不会反对。但我想要的是那种静态的版本。
多维 for_each 提案与开放性问题¶
我谈论了这么多,现在我要谈谈多维 for_each。我们已经有了一个提案。它结合了许多前人的努力。Ruslan 在并行的 ranges 算法上的工作。David 修复了 ranges 的工作,以便它们可以与并行算法一起工作,谢谢。并且你现在可以将它们分解为迭代器(iterators)和哨兵(sentinels)。你看,这里的很多东西直到最近才真正起作用。还有 views::iota 提案,那真是太棒了。它修复了 views::iota,这样如果我把零作为下界,它不会把我的范围变成一个推导的 int。所以那个 views::indices,那太棒了。你必须用那个。
好吧,我在左边展示了多维的 for_each。这是一个三维循环。我完全不想敲出左边那长长的一串代码。我不想。而且这里面牵扯到某位法国哲学家(指笛卡尔 Descartes)。但你连他的名字都拼不对,这是搞什么?这个 cartesian 到底是个什么鬼东西,对吧?而且作为用户,你必须手动去解包(unpack)这个索引元组。而且它只是默认行主序(row major order)的。你可以通过各种手段按摩你的笛卡尔积来获得列主序或其他东西,但那太痛苦了。
还有其他问题。比如如今的 ranges,与 transform_view 以及类似的东西一起使用时,在某些情况下——甚至当 lambda 闭包本身是平凡可复制(trivially copyable)的时候——它们却不是平凡可复制的。比如,我有一个 lambda,我按值捕获一个 int。那是平凡可复制的,只是不能被赋值。但我把它放进 views::transform 里,现在它就不再是平凡可复制的了!那真没劲,因为这样我就不能在加速器(accelerators)上运行它了。我就不能把它 memcpy 复制过去了。所以 Ruslan 正在主导一个提案来修复这个问题。
此外,由于 ranges 库并不是真正为高度可定制而设计的,我无法通过特化代码来进行可移植的优化。如果我手里有一个 cartesian_product_view(笛卡尔积视图),我无法获得它的子视图(sub views)。因此,我无法对笛卡尔积编写我自己的算法定制。我想编写的其中一个定制可能是:我注意到我拥有一个这些索引视图、iota_view 的笛卡尔积。这是一个多维循环。我为什么不能把它拆开,把循环边界提取出来,然后将其重写为一个嵌套的 for 循环,并在上面加上我的 OpenACC 预编译指令(pragmas)或 OpenMP 预编译指令,或者编译器特定的扩展指令呢?因为编译器知道如何很好地优化那种形式的代码。如果我能深入到笛卡尔积的内部,获取所有的范围,我就可以编写做这些事的特化代码。但我无法做到,因为除非我使用非标准的方法,否则我无法进入笛卡尔积内部。
假设我是标准库的实现者。我非常想编写那样的特化代码。但那样的话,我就需要在我的笛卡尔积内部放入一些非可移植的小钩子(hooks)。然后我的标准库就变成了专属于我自己的“特制版”。然后我的编译器就只能和我自己的标准库一起工作。并且我必须重新实现整个 ranges 库才能做到那种特化优化。那毫无乐趣可言,我不应该被迫去这样做。
以上这些问题,就是我们提出多维 for_each_index 提案的动机。我们不用去写那一长串难以优化和难以特化定制的东西,而是写一种很短的表达方式。它接受一个布局映射作为参数,因为那为我提供了一个首选的迭代顺序提示(hint)。
我之前稍微谈到过循环分块(loop tiling)。在基于编译指令和编译器扩展(如 OpenACC 或 OpenMP)的模型中,我可以告诉编译器我想要将循环打碎成一个个小块。那是为时间局部性而重组循环的一种方法。但如果你使用的是一维(1D)的方式去访问多维数据,这样做就不太有用了。我稍后会讲这一点。为了回到我们提到的“隐式与显式并行的层级结构”话题:对于多维 for_each,我可以从基于数组的方法开始。但如果我能得到编程语言级别的支持,基于数组的方法会好得多。
我在左下方紫色的框里展示了 Z-level 编程语言中的一个例子。那是一个早期的基于数组的并行编程模型。在这个模型里,如果你想说:“我想把索引范围移到北边”,你可以直接说 north。太棒了。如果我(用目前的 C++)像右边一样有多维 for_each_index 提案中的索引,我可以直接通过直接访问数组来表达这一点。虽然这可能在库的更高级别上更难优化,但至少我可以说出我想表达的逻辑。
针对多维算法的按元素操作(Element-wise)可能会有点奇怪。就像我试图在橙色框里展示的那样,这种模板模板操作(stencil operation)的 element-wise 版本会变得非常难看。它并不是一种自然契合的表达。所以你要么走基于数组的路线,要么走基于索引(index-based)的路线。而我们在这个提案里选择了基于索引的路线。
我可以谈论一下多维块(multidimensional bulk)算法。如果是作为我们讨论过的同步算法的异步版本,这没问题。但实际上除了多出维度并没有什么新东西。bulk 的最初意图实际上就是处理多维的。所以如果你回头看看早期的提案版本,那里之所以被称为 shape(形状)而不是索引或大小,就是因为它本意指代一个形状。把它带回来可能并不难,但我认为也许我们应该先专注于处理一维异步算法。就像 Bryce 的提案中所讨论的那样。
还有其他的并行多维算法。我们谈过 copy、fill、for_each_index。你还可以想象出一些其他的。当你有了多维数组,你自然会催生出新种类的算法。因为你有了更多的边界范围,以及执行并行的更多方式。那些是基于数组的全局算法。也许你会想使用专门的数组库来处理这些。目前尚不清楚它们最终进入标准时应该呈现什么形式,但业界有很多这样的库可以选择。Connor Horkster 的 parrot 就是个不错的库。Alfredo 的库也是另一个好选择。所以这里对于基于数组的库有很多选择,我不确定什么才真正属于标准之内。
多维 SIMD 与 QTile¶
在结束之前,我想谈谈多维 SIMD。它就像是 SIMD,只是变成了多维的。听起来好像有点搞笑。我所说的多维是什么意思?正如我之前提到的,多个维度会带来新类型的算法。有降秩归约(rank-reducing reductions),有诸如矩阵-矩阵乘法的收缩计算(contractions)。有行和列的求和,这些是在各种子边界范围上的归约操作,它们会产生一个秩(rank)较低的结果。还有可以重塑(reshape)数据形状的方法,因此不是一个 4 x 2 的形状,它看起来像一个 1 x 8 的形状。或者我可以取转置,把它翻转过来。
但因为这是 SIMD,它必须受到某些限制。它不是一个纯粹的数组库,它是一个 SIMD 库。所以也许你想限制其维度只能是某种特定的形状,也许你希望它们必须是 2 的幂。也许你希望加载或存储操作必须是对齐的(aligned)。也许你只想支持特定种类的元素类型,或特定的归约操作。你希望数组的布局是不透明的(opaque)。你不希望用户能够直接获取指向其内部的指针,那是不酷的。
我们现在刚刚发布了一个叫做 QTile 的东西。那是多维 SIMD 的一种体现或版本。QTile 是一个基于数组编程语言的家族。我们在 12 月发布了 Python 版本。之后会有 Julia 版本和 Rust 版本。甚至还有一个 Basic 版本,那是 Bryce 帮忙发布的,我们做那个纯粹是为了好玩。但我想说的是,为它编写不同的语言绑定并不难。而且 C++ 版本很快就会推出。
这些都是基于一种虚拟机(virtual machine)、一种虚拟指令集的。它的设计初衷是可移植性。其理念是,如果你喜欢用 APL 编程并且你想要一个 QTile 版本的 APL,你可以实现一个能将你的 APL 编译为 Tile IR 的编译器。然后你马上就能得到 QTile APL。QTile 的工作方式很像 CUDA。用户用宿主语言的一个子集编写 tile kernel(瓦片内核)。这些 kernel 在 GPU 上运行。然后你就能获得全局算法的并行化。你显式地获得了全局的并行性。它就像多维的 bulk,你在一个代码块构成的网格(grid of blocks)上运行。在每个 block 内部,代码看起来是顺序执行的,但它实际上在使用这些隐式并行的 tiles(瓦片块)。
所以,其中有一些看起来像 mdspan 的东西,它们是全局的输入和输出数组。我们没有直接使用 mdspan,因为在这些东西上有诸多限制。你不能直接访问里面的元素,你必须通过视图(views)去访问这些元素。这里的 Tile 也就是我们所说的多维 simd_vec。你可以创建不同种类的 tiles。你可以推理形状,它们全都是编译时的产物。这里的 Views(视图)用于声明去执行经过对齐的、平铺好的(tiled)加载和存储操作的意图。你不用 views 也可以做加载和存储,但 views 能够帮助声明这种访问意图,因此它允许编译器开启一些优化。也会有其他种类的 views,它们将支持平铺加载和存储之外的其他访问模式(例如针对模板计算的 strided 模式)。
BID 是多维块的索引(block index)。那是某种索引范围构造,它的存在是为了帮助我们优化 for 循环。但这并不是强制需要的。里面还有矩阵乘法和其他类型的计算。
未来将会有很多关于 QTile 的演讲。C++ 版本即将发布。Bryce 和 Steven Jones 他们的演讲展示了 QTile C++ 的用例。
关于命名,我可能想谈谈为什么我们称其为 tile(瓦片)。或者说,你打算把多维的 simd_vec 叫什么?拜托不要叫它 mdvec,千万别。而且它也不是原子的(atomic)。你脑海中想象的那个东西,它并不是一个普通的数组,它只是负责对别的东西进行划分(partition)而已。叫它 cell(细胞/单元格)也许勉强可以,但它又不是电子表格程序,它没有那种模型。tile 这是一个很好的名字,因为它划分并覆盖了某些东西。brick(砖块)也不错,你可以用它来构建结构。所以大家好好想想命名。你可能会想要全局的多维 SIMD 算法,我们可以考虑。
掩蔽操作(Masking)在遍历多维 SIMD 对象时非常重要。因为你总会遇到“尾巴(tail)”问题。使用 1D 的 SIMD 时,你的尾部只是一维的东西,它的大小永远不会超过你的 SIMD 最大宽度。但是,你拥有的维度越多,你拥有的“尾巴”就越多(例如多维边界的零散元素)。这使得掩蔽变得尤为关键,不仅仅是为了编写算法的便利,实际上这也是为了性能。
最后我想谈谈的是,你拥有所有这些“耦合(coupling)”,对吧?SIMD 归根结底全是关于将硬件底层、内存操作与微小的数据并行体耦合在一起。所以里面充满了很多整数常量(integer constants)。也许有各种算法的变体,也许你需要 在运行时(runtime) 做出这些选择。因此,你需要能在运行时进行编译。如果能有即时编译(JIT, Just-In-Time compilation),那就太好了。Hal Finkel 很久以前就曾有过在 C++ 中进行即时编译的提案。我认为,也许现在的**反射(reflection)**特性改变了我们思考该语法的方式。请大家思考一下可以用于这方面的各种不同的表示法。
总结与问答环节¶
总结一下,今天我给大家留下了很多悬而未决的问题(open questions)。我们谈到了 restrict 关键字、不相交性、多维数组。多维 SIMD 究竟需要多少编译器级别的支持?我们的 QTile 本质上就是一个编译器,它严重依赖编译器的支持。它有可能纯粹作为一个库的特性存在吗?但与此同时,由于我个人的背景不是搞编译器的,我以前很害怕写编译器。但也许我现在不需要害怕了。现在写一个编译器,比我当年成长(我把读研称为成长)的时候要容易得多了。所以也许我们不必再对编写编译器感到恐惧。即时编译(JIT)也是一个巨大的开放性问题。
所以总结来说,我展示了当前标准支持一维并行的所有方式。然后很自然地把那些对应到了多维并行的支持上。正是这些方式启用了不同的属性特性。因为并行不仅关乎“并行算法本身”,它更是关乎具备去进行并行的能力。这正是为什么我花了大量篇幅去谈论不相交性、局部性和耦合。
所以,对于想要起草提案或者修改提案的人来说,还有很多事情等着我们去做。还有很多东西等着去探索。但我希望确保大家不要迷失,不要失去对“多维性”这个全局大图景(big picture)的关注。我做这个演讲的部分目的,就是为了传达这个宏观的图景。
我非常喜欢埃米莉·狄金森(Emily Dickinson)这句关于搬家的名言。如果你曾经搬过家,你可能会对她的感受产生共鸣,那种拼命想要弄清楚所有东西都在哪里的感觉。
“I took it a time in the memorandum of my several senses. And also of my hat and my coat and my best shoes. But it was lost in the melee. And I’m out with lanterns looking for myself.”
(我曾在感官的备忘录中仔细清点,还有我的帽子、我的外套和我最好的鞋子。然而,它们都在混乱中遗失了。于是我提着灯笼走出去,去寻找我自己。)
所以,大家请看好你们的“帽子”和“鞋子”(译注:不要在多维并行的复杂技术细节中迷失方向)。谢谢大家!
是的,我不想阻挡大家去喝咖啡的步伐。很抱歉我拖堂了几分钟,但如果大家有任何问题,我很乐意回答。
Ruslan 是第一个提问的。
(问答环节开始)
Ruslan: 是的,关于这个 MD……不管叫什么,多维 SIMD。我想知道除了你的硬件,还有没有任何其他硬件能与之对应映射?
Mark: 所以说,思考这个问题也许会成为我未来工作的一部分。我认为,比如 Clang 具有矩阵类型(matrix type)扩展。所以存在着一个我们在编译器中可以用来进行实验的框架。而且,比如那些做图形编程的人,他们通常都在使用 4 x 4 或 3 x 3 的小矩阵工作。所以我可以想象他们是有硬件支持的,或者是有把那些操作映射到 SIMD 的途径的。其实 SIMD 最初始的动机就是图形学,对吧?而图形学里充满了矩阵,所以……
Ruslan: 是的,我的意思是……好吧,GPU。也许不仅仅是你们的 GPU。但我觉得它不能立即映射到 CPU 上。
Mark: 如果你的多维 SIMD 的大小是 2 的幂,那么你就可以把它分解成多个局部的、一维的 SIMD 操作块来执行。
Ruslan: 是的,但这并不是直接在多维寄存器上的映射。我是这么认为的。
Mark: 但即使是我们在 GPU 上的映射,也不一定是绝对“直接”的。因为我们同样必须去管理我们的局部内存。所以,这确实是个好问题。是的。这是个很有意思的开放性问题。
好了。还有什么问题吗?看来大家都迫不及待想去喝咖啡了。我完全理解。好吧。如果稍后有任何问题,欢迎来找我交流。谢谢。
(掌声与道谢)
太棒了,Jim。谢谢。谢谢。非常感谢。太棒了。