研究文章

开放式访问

矩阵转置的基于缓存的Hilbert曲线分块方案

作者:

文章编号：37，页数1-28

https://doi.org/10.1145/3555353

出版:2022年12月19日出版历史

所有格式 PDF格式

摘要

本文提出了一种快速SIMD-Hilbert空间填充曲线生成器，该生成器支持一种新的缓存遗忘阻塞方案技术，该技术应用于一般矩阵的错位换位。高性能计算库中的矩阵操作通常基于主机微处理器规范进行参数化，以最小化不同级别内存层次中的数据移动。缓存可发布算法的性能不依赖于此类参数化。这种算法提供了一种优雅且可移植的解决方案，以解决现代处理器中缺乏标准化的问题。我们的解决方案包含一种迭代阻塞方案，该方案利用希尔伯特空间填充曲线的局部保持特性，最大限度地减少任何内存层次中的数据移动。该方案遍历输入矩阵O（纳米）时间和空间，改进了固有的内存局部性较差的矩阵算法的行为。与最先进的方法相比，将此技术应用于非现成矩阵转置问题取得了具有竞争力的结果。在采用标准软件预取技术后，我们的解决方案的性能超过了“英特尔MKL”版本。

1引言

在计算机科学的大多数领域中，人们会发现矩阵运算是不可避免的原始结构之一，其运算允许计算其他更复杂的算法。处理大型矩阵时最严重的瓶颈是将矩阵数据从主内存提取到缓存所带来的开销。为了实现高效的实现，程序员应该使用缓存友好的解决方案，利用缓存的时间和空间局部性。几个科学计算库以近乎最佳的运行时实现了矩阵运算的缓存友好解决方案[12,18,24]. 这些库通过微调一组参数来有效利用缓存，这些参数的值取决于缓存规范，例如缓存线大小和整体内存层次结构。这些解决方案不可移植，因为缓存规范没有标准化。每个机器架构都需要自己的解决方案。

1.1缓存可用算法

然而，有些高效缓存算法的性能与缓存规范无关。这些算法被称为cache-obliovious[11]. 这组算法通常通过递归的分治方法获得高缓存命中率。在给定的时间点，每个子问题都将适合缓存。在现代CPU上，由于其递归特性，这些解决方案很少比缓存软件更有效。获得缓存遗忘解决方案的一种方法是利用空间填充曲线（如Hilbert、Peano和Z阶曲线）的固有局部保持特性。应用填充曲线优化密集矩阵运算是一个广泛研究的领域[2,三,16]. 两个希尔伯特[13]和Z顺序[20]曲线被定义为能够穿越任何\（2^l\乘以2^l\）网格，其中\（l\in\mathbf{N} _0（0） \)，每个条目只访问一次。这确保了双射函数的存在\（\mathcal{H}^{-1}（H）\）和\（\mathcal{Z}^{-1}（Z）\）映射曲线点\（小时\）或\（z）变成二维点\（（i，j）\）希尔伯特曲线和Z阶曲线。值得注意的是，用于实现这些曲线的大多数方法也是递归的，与递归的分治方法一样效率低下。

在最近的一项工作中，Böhm等人。[5]提出了FUR-Hilbert，一种基于Hilbert空间填充曲线的网格遍历方案。他们的方法使用小的预处理遍历路径，称为纳米投影图，然后根据希尔伯特曲线的给定迭代进行细分。他们的解决方案为每个循环迭代提供了一个恒定的开销，并且是完全迭代的，从而产生了一个高效的策略，可以很容易地由编译器进行优化。此循环的位置保持特性取决于宿主算法中使用的原始矩阵遍历方案。由于希尔伯特填充曲线的性质，此循环只能用于不存在任何类型的数据依赖性的算法，如矩阵乘法。

1.2缓存式异地换位

与其他密集矩阵运算相比，矩形矩阵的缓存式异地转置似乎没有受到科学界的太多关注。Chatterjee等人。方阵同位转置的不同计算方法比较[7]. 这些方法包括基于分治或填充空间曲线的缓存策略。提供最佳执行时间的方法实际上是一种利用Z阶曲线的局部保持特性的策略。虽然他们的解决方案不能直接转化为这个问题，但这表明可以通过使用类似的策略来改进矩形矩阵的非现成转置。

矩形的非局部换位问题\（n倍m）矩阵\（M\）通常定义为写入\（米）作为\（m\次n\）矩阵\（T\）即。，\（对于所有_{i，j}\ T_{j，i}：=M_{i，j}\），其中\（0\le i\lt n\）和\（0\le j \lt m）。此问题通常使用基于以下嵌套的RAM模型方法来解决：

这种类型的解决方案提供了一种更简单的嵌套操作，它引入了最少的算术指令，但忽略了内存层次结构的行为。此循环意味着\（M\）或\（T\）都是在row/column-major顺序布局中分配的，将使用column/row-major顺序遍历进行扫描。对于足够大的矩阵，这会导致不必要的容量和冲突缓存丢失，甚至页面错误。矩阵的转置是一个内存绑定问题。对于这种类型的问题，通常有利于改进内存访问模式，而不是减少算术指令的数量。通过将此方法与简单缓存模型中的不同访问模式进行比较，可以直观地了解标准的非现成转置算法的性能降低效果。假设我们的缓存具有以下特征：

•

它由四条线组成，每条线可以容纳两个单精度值；

•

它是完全关联的，意味着不会发生冲突遗漏；

•

它采用最近租赁的更换政策；

•

它还采用了写入分配策略，即先读后写。

让\（M\）和\（T\）是\（4乘以4）使用row-mahor布局在内存中分配单精度矩阵。\（M\）对应于问题的输入矩阵，而\（T\）对应于输出矩阵。用于计算\（M\）由以下嵌套式表示：

在前面描述的缓存上运行此策略将导致总共24个缓存未命中，如图所示1（a）和1（d）。有人可能会认为，在\（T\）是缓存较小的直接结果。事实并非如此。任何以行主顺序存储并随后使用列主顺序扫描的矩阵都将导致较差的空间局部性保存。每当矩阵列数大于缓存线大小且矩阵行数大于缓存行数时，对输出矩阵的内存访问很可能会导致缓存丢失。

图4。

为了避免遍历\（T\）使用列主顺序，可以使用图1（b）或1（c）所示的方案，该方案利用希尔伯特或Z阶逆函数，表示为\（\mathcal{H}^{-1}\）或\（\mathcal{Z}^{-1}\），在输入矩阵上施加基于Hilbert或Z阶空间填充曲线的遍历\（M\），作为：

两种遍历方案都确保输出矩阵\（T\）将通过之前的希尔伯特或Z阶空间填充曲线的反射进行扫描，如图所示1（e）和1（f）分别在上面定义的缓存上模拟这些策略只会导致16个缓存未命中。通常，从主内存读取数据所需的时间约为70 ns，而从缓存加载数据只需要2 ns，如Arroz等人所述。[1]. 如果不考虑从缓存（cache-hit）和主内存（cache-miss）加载时间以外的任何成本，这种简单的方法会带来以下开销\（70乘以24+2乘以8=1696）ns，而其他基于空间填充曲线的方法只需要\（70乘以16+2乘以16=1152）纳秒。如果忽略这些操作的算术代价，则可以得出结论，后一种策略的性能加速接近\( 33\% \).

1.3贡献和文件组织

本文件提出了两个主要贡献：

(1)

在节中三，我们引入了一个线性希尔伯特填充空间曲线生成器，它不仅是迭代的，而且没有数据依赖性，它允许实现一个主要使用SIMD指令的解决方案，并且可以由大多数编译器进行有效优化。我们的生成器要求将曲线存储在内存中，这意味着只要需要两代或更多代相同的曲线，就可以摊销构建此曲线的成本。该算法的正确性证明基于第节中定义的Lindenmayer系统2.

(2)

在节中4在FUR-Hilbert的基础上，我们开发了一种替代分块方案，该方案能够利用矩形网格上Hilbert和Z阶空间填充曲线的局部保持特性。该策略的总体思想是将输入矩阵划分为\（2^l乘以2^l）较小的子矩阵。然后，这些子矩阵的遍历顺序由一个或多个希尔伯特或Z阶空间填充曲线定义，如图所示2（a）以及2（b） ●●●●。与原始方法相反，我们的分块方案允许使用标准预取技术，这些技术会积极影响非现成转置问题的性能。

两个模块的方案示例\（9乘以17）基于（a）希尔伯特填充曲线和（b）Z阶曲线的网格。每个块，子矩阵，使用行主顺序进行扫描。

在节中5，我们将SIMD Hilbert空间填充曲线生成器与已知计算线性时间内Hilbert曲线的其他解决方案进行了比较，以了解曲线元素的数量。在本节中，我们将在矩形矩阵的非现成转置的背景下对我们的块方案进行广泛的分析，并将我们的缓存可用策略与当前的技术水平（包括缓存软件解决方案）进行比较。在C编程语言中，我们的缓存发布方法的实现可以在以下位置找到：https://github.com/JoaoAlves95/HPC-Cache-Oblivious-Transposition网站.

计算希尔伯特填充曲线的两种已知方法

有几种计算迭代的经典方法\（l）希尔伯特填充曲线。然而，此类算法通常会带来以下时间成本\（O（\log（N））\）每个条目[6,8,15]. 与传统的嵌套for-loop（枚举\（O（N）\）即。，\（O（1）\）每个条目。

林登迈耶体系[22]是一个基于无上下文语法的并行重写系统[9]，能够在线性时间内生成希尔伯特填充曲线的任何迭代\（O（N）\）。该系统接收所需的扩展级别作为输入\（l）并返回中的符号序列\（\lbrace\rightarrow，\\uparrow、\\leftarrow和\\downarrow\rbrace\）描述通过迭代执行的遍历\（l）给定的希尔伯特曲线。每个符号都描述了沿水平或垂直坐标轴的单位步，其中

\（\begin{gather}\rightarrow\text{对应于}j:=j+1，\end{gather{\）

(1)

\（\begin{gather}\leftarrow\text{对应}j:=j-1，\end{gather{\）

(2)

\（\begin{gather}\\uparrow\\text{对应}i:=i+1，\end{gather{）

(3)

\（\begin{gather}\\downarrow\\text{对应于}i:=i-1。\结束{聚集}\）

(4)

为了保持本文件的完整性，我们对贝德[2]定义该系统，如定义中所示1此定义将在第节中有用三以证明我们战略的正确性。

定义1。

希尔伯特空间填充曲线的林登迈耶系统。

设计用于生成希尔伯特空间填充曲线的标准Lindenmayer系统由以下部分组成：

(1)

一组非终结符：\（轨道H、轨道A、轨道B、轨道C）每个非终结符表示图中所示的四种基本模式之一三。出于我们的目的，我们将假设公理，即起始符号，是\（H）或\（A\）.

(2)

一组端子符号：\（\lbrace\rightarrow、\\uparrow和\\leftarrow，\\downarrow以及\\varepsilon\rbrace\）。每个箭头符号代表其所描绘方向上的一个单位步。航站楼\（\varepsilon\）对应于空符号，它用于从语法生成的最后一个字符串中删除任何非终结符，从而生成只包含终结符的字符串。

(3)

一套生产规则：

\（开始{align*}H&#x0026；\xrightarrow[]{{\it产生}}\varepsilon \|\A\\uparrow\H\rightarror H\\downarrow/B，\\A&#x0025；\xrightarrow[]{\it生成}}\varepsilon \ |\H\ rightarrol A\\leftarrowC，\\B&#x0027；\x右箭头[]{\ it产生}}\varepsilon\|\C\左箭头B\\下箭头\B\右箭头H，\\C&#x0026；\xrightarrow[]{{it生成}}\varepsilon\|\B\\downarrow\C\leftarrowC\\uparrow\ A\end{align*}\）

(4)

附加扩展规则1：要构造与希尔伯特曲线迭代相对应的此语法字符串，必须立即扩展当前扩展中存在的所有非终结符。

(5)

附加扩展规则2：我们只允许将非终结符扩展为\（\varepsilon\）一旦我们达到了这个语法的预期扩展水平。此时，所有非终结符必须替换为\（\varepsilon\）.

图4说明了第二次迭代的构造\（l=2）Hilbert空间填充曲线通过Lindenmayer系统的公理\（H）在这一点上，我们有生成线性时间内Hilbert曲线的任何迭代所需的所有信息，即。，\（O（1）\）每个网格条目的成本，以及的内存成本\（O（\log _4（N））\）在函数堆栈中。这种方法看起来很有前途。然而，正如Böhm等人指出的那样，它的递归性质不允许现代编译器高效地生成优化的代码。[5]. 为了解决这个问题，这些作者通过模仿递归Lindenmayer系统的堆栈行为，开发了第一个非递归Lindenmayer系统。

图6。

图7。

他们的方法非常有效\（O（N）\）时间和恒定的内存复杂性。然而，它提供的数据依赖关系限制了SIMD指令的使用，这意味着如果能够解决或缩小这些依赖关系，则可以实现更快的解决方案。

三SIMD（单指令多数据）希尔伯特填充曲线生成器

为了进一步减少第节算法中的开销2，我们的解决方案必须是完全迭代的，并在线性时间内运行，从而生成可以用SIMD指令轻松扩展或由现代编译器自动矢量化的代码。这种解决方案可以通过基于位向排他或运算的线性和增量希尔伯特填充曲线生成器实现，从而实现无数据依赖性的方法。

为了更好地掌握构建这种生成器的基本概念，了解用于表示Hilbert空间填充曲线的编码以及如何转换当前方向很重要\（d\）到目前为止，希尔伯特填充曲线被描述为一系列方向，其中每个方向都由中的一个终端符号表示\（\lbrace\rightarrow、\uparrow，\leftarrow和\downarrow\rbrace\）.使用此符号，方向序列(\（\texttt{曲线}\）)描述第二次迭代的\（l=2）由公理产生的希尔伯特填充曲线\（H）将表示为

\（\[\texttt{curve}=\\rightarrow\\\\uparrow\\\leftarrow\ \\uparrolrow\\right箭头\\downarrow\\n向上箭头\\rightarrow\\down箭头\\\\downarrow\\n向下箭头\\left箭头\\right箭头\\]\）

当然，这种编码不适合将来对该曲线进行算术操作和分析。为了解决这个问题，Böhm等人。[5]提出了一种替代编码，该编码映射\（k\）中的第个端子符号\（\lbrace\rightarrow、\uparrow，\leftarrow和\downarrow\rbrace\）到\（k\）第th个整数\（\lbrace 0,1,2,3\rbrace）此编码提供了将曲线表示为一个以4为基数的数字的数学方法，其中每个方向由一个数字表示，方向序列从左到右读取。将此编码应用于之前的方向序列将导致

\（\[\texttt{curve}=\underrightarrow{\0\1\1\1\\0\1\0\3\3\2\0}_4，\]\）

哪里\（\右下箭头{}_4\）表示曲线由以4为基数的数字表示，以及处理方向的顺序。如前所述，此编码还提供了转换给定方向的方法\（d\in\lbrace\underrightarrow{\0\}_4，\underright arrow}_1，\undertrightar罗w{2\}_4和\underrightarrow{3\}_4\rbrace\）进入尊重方程式的单元步骤(1)至(4)通过使用余数运算符(\( \% \)¹)，作为：

\（开始{align*}j:=j-（（d-1）\%2）；\\i： =i-（（d-2）\%2）；\结束{align*}\）

3.1计算曲线象限

图5表明迭代的第二象限和第三象限\（l）可以通过沿\（y=x）线，而迭代的第四象限\（l）可以通过将该曲线的第一个象限旋转180°来计算。如图所示，似乎可以将反射和旋转操作编码为位互斥或1的应用程序(\（\texttt{xor}_1 \))和逐位互斥或乘2(\（\texttt{xor}_2 \))该曲线第一象限包含的每个方向：

\（\begin｛align｝\texttt{xor}_1（\右下箭头{\0\1\}_4）&#x0026=\underrightarrow{1\0\3\}_4=\\uparrow\\rightarror\\downarrow\，\end{align}\）

(5)

\（\开始{align}\texttt{异或}_2（\右下箭头{\0\1\}_4）&#x0026=\underrightarrow{\2\3\0\}_4=\\leftarrow\\downarrow\\nrightarror。\结束{align}\）

(6)

图8。

这种直觉是由命题形式化的2和三.

提议2。

应用\（\texttt{xor}_1 \)给定希尔伯特空间填充曲线的第一象限的结果是相同迭代的第二象限和第三象限。

证明。

假设我们的曲线是由以定义为特征的Lindenmayer系统生成的1这个系统的公理要么是非终结性的\（A\）或\（H）. The\（k\）给定希尔伯特曲线的第四象限的特征是\（k\）通过将任何公理扩展一次而获得的字符串中存在的第th个非终端符号，即在级别上获得的字符串\（l=1）.

\（\begin{align}H&#x0026；\xrightarrow[]{{it生成}}A\\uparrow\H\rightarror H\\downarrow\ B\end{alinge}\）

(7)

\（\begin{align}A&#x0026；\xrightarrow[]{{it生成}}H\rightarror A\\uparrow\A\leftarrow C\end{alinge}\）

(8)

这一事实使我们能够正式重申提案2作为\（\texttt）{xor}_1（H_l）=A_l\楔形\texttt{xor}_1（B_l）=C_l\），其中\（l\in\mathbb{N}^+\）表示扩展步骤的数量。提议2现在将通过结构归纳法证明。

归纳假设： \（\texttt）{xor}_1（H_l）=A_l\楔形\texttt{xor}_1（B_l）=C_l\），其中\（H_l\）,\（A_l\）,\（B_l\）和\（C_l\）表示通过扩展此非终结符生成的曲线\（l）时间，对于\（l\in\mathbb{N}^+\）.

基本情况：仅在第一个扩展步骤中出现非终结符时进行扩展，\（l=1），每个公理的结果

\（开始{align*}H_1&#x0026；\xrightarrow[]{{\it产生}}\\uparrow\\rightarror\\downarrow\=\underrightarrow{\1\0\3\}_4，\\A_1&#x0.026；\xrightarrow[]{\it生成}}\\rightarrow\\uparror\\leftarrow\=\underrightarrolrow{\0\2\}_4、\\B_1&#x00026；\x右箭头[]{\\it产生}}\\左箭头\\下箭头\\右箭头\\=\右下箭头{\2\3\0\}_4，\\C_1&#x0026；\xrightarrow[]{{it-Products}}\\downarrow\\leftarrow\\nuparrow\=\underrightarrow{3\2\1\}_4。\结束{align*}\）

因此\（\texttt｛异或｝_1 \)非终端扩展的第一步\（H）和\（B\）方程中所示方向序列的结果(9)和(10)分别是。这是

\（开始{方程式}\texttt{xor}_1（H_1）=\texttt{xor}_1（\underrightarrow{\1\0\3\}_4）=\underrightarrow{\0\2\}_4=A_1，\end{方程式}\）

(9)

和

\（开始{方程式}\texttt{xor}_1（B_1）=\texttt{xor}_1（\underlightarrow｛\2\3\0\｝_4）=\underlightarrow｛\3\2\1\｝_4=C_1。\结束{方程式}\）

(10)

可以观察到基本情况成立。非终结符的一个扩展步骤产生的方向序列\（A_1\）等于\（\texttt{xor}_1（H_1）\），而非终结符的一个展开步骤产生的方向序列\（C_1\）等于\（\texttt{xor}_1（B_1）\）.

感应步进：应用\（\texttt{xor}_1 \)非终端符号产生的符号序列\（H）和\（B\）方程式中的结果

\（\begin｛equation｝\texttt）{xor}_1（H_1）=\texttt{xor}_1（A_{l-1}\向上箭头H_{l-1}\右箭头H_}-1-1}\向下箭头B_{1-1}），\结束{方程式}\）

(11)

和

\（\begin｛equation｝\texttt）{xor}_1（B_{l}）=\texttt{xor}_1（C_{1-1}\左箭头B_{l-1}\下箭头B_}-1}\右箭头H_{l-1}）。\结束{方程式}\）

(12)

请注意\（\texttt{xor}_1 \)可以在序列上逐元素地应用，即，将该函数应用于方向序列或应用于该序列的每个单独元素保持相同的结果。因此，可以进一步简化方程(11)和(12)作为方程式(13)和(14)分别为，

\（开始{方程式}\texttt{xor}_1（H_{l}）=\texttt{xor}_1（A_{l-1}）\\texttt{xor}_1（向上箭头）\\texttt{xor}_1（H_{l-1}）\\texttt{xor}_1（向右箭头）\\texttt{xor}_1（H_{l-1}）\\texttt{xor}_1（向下箭头）\\texttt{xor}_1（B_{l-1}），结束{方程式}

(13)

和

\（开始{方程式}\texttt{xor}_1（B_{l}）=\texttt｛异或｝_1（C_{l-1}）\\texttt{xor}_1（\左箭头）\\texttt{xor}_1（B_｛l-1｝）\\texttt{xor}_1（向下箭头）\\texttt{xor}_1（B_{l-1}）\\texttt{xor}_1（向右箭头）\\texttt{xor}_1（H_{l-1}）。\结束{方程式}\）

(14)

提醒终端符号\（（向右箭头）,\（向上箭头）,\（\左箭头\）,\（\向下箭头）\）编码为\( (0,1,2,3) \)。这允许我们重写\（\texttt{xor}_1 \)到方程式中的端子符号(13)和(14)作为

\（开始{方程式}\texttt{xor}_1（H_{l}）=\texttt{xor}_1（A_｛l-1｝）\\rightarrow\\texttt{xor}_1（H_{l-1}）\\uparrow\\texttt{xor}_1（H_{l-1}）\\leftarrow\\texttt{xor}_1（B_{l-1}），结束{方程式}

(15)

和

\（\begin｛equation｝\texttt）{xor}_1（B_{l}）=\texttt{xor}_1（C_{l-1}）\\向下箭头\\texttt｛异或｝_1（B_{l-1}）\\leftarrow\\texttt{xor}_1（B_{l-1}）\\向上箭头\\texttt{xor}_1（H_{l-1}）。\结束{方程式}\）

(16)

我们的归纳假设意味着\（\texttt{xor}_1（H_｛l-1｝）=A_｛l-1｝\）和\（\texttt{xor}_1（B_{l-1}）=C_{l-1}对于\（l\gt 1）.自\（\texttt{xor}_1 \)是幂等函数，它还暗示\（\texttt{xor}_1（A_{l-1}）=H_{l-1}和\（\texttt{xor}_1（C_{l-1}）=B_{l-1}根据我们的归纳假设，可以再次简化方程(15)和(16)作为

\（开始{方程式}\texttt{xor}_1（H_1）=H_{l-1}\右箭头A_{l-1}\\上箭头\A_{1-1}\左箭头C_{1-1}=A_l，\结束{方程式}\）

(17)

和

\（开始{方程式}\texttt{xor}_1（B_l）=B_{l-1}\\downarrow\C_{l-1}\ leftarrow C_{1-1}\\uparrow\ A_{1-1}=C_l.\end{方程式}\）

(18)

请注意，在\（\texttt）{xor}_1（H_l）\）和\（\texttt{xor}_1（B_l）\）与的第一个扩展步骤中出现的非终端相同\（A_l\）和\（C_l\）分别为。命题证明到此结束2.□

提案3。

应用\（\texttt{xor}_2 \)希尔伯特空间填充曲线的给定迭代的第一象限导致相同迭代的第四象限。

与命题证明类似2，假设我们的曲线是由以定义为特征的林登迈耶系统生成的1而且那个\（k\）给定希尔伯特曲线的第四象限的特征是\（k\）字符串中出现的第个非终结符，该字符串通过将任意公理展开一次而获得：

\（\begin{align}H&#x0026；\xrightarrow[]{{it生成}}A\\uparrow\H\rightarror H\\downarrow/B，\end{aling}\）

(19)

\（\begin{align}A&#x0026；\xrightarrow[]{{it Produces}}H\rightarror A\\uparrow\A\leftarrow C.\end{aling}\）

(20)

命题陈述三等于\（\texttt{xor}_2（A_l）=B_l\楔形\texttt{xor}_2（H_l）=C_l\），其中\（l\in\mathbb{N}^+\）表示扩展步骤的数量。提议三可以通过遵循命题中使用的策略进行结构归纳来证明2; 为了保持这份文件简短，我们选择省略这个证明。

3.2连接曲线象限

至此，我们已经了解了如何生成希尔伯特填充曲线给定迭代的最后三个象限。要完成此曲线的生成，需要连接之前计算的象限。用来连接这些象限的方向可以在这个语法扩展的第一步中找到。来自命题2和三接下来是扩张的下一步，\（l+1），可以根据当前的扩展步骤生成\（H_l\）和\（A_l\）作为

\（开始{方程式}H_{l+1}=\texttt{xor}_1（H_l）\向上箭头H_l\向右箭头H_l \向下箭头\texttt{xor}_2（\texttt{xor}_1（H_1）），结束{方程式}

(21)

和

\（开始{方程式}A_{l+1}=\texttt｛异或｝_1（A_l）\右箭头A_l \上箭头A_l\左箭头\texttt{xor}_2（\texttt{xor}_1（A_l））。\结束{方程式}\）

(22)

虽然此方法生成的方向序列与经典林登迈耶系统生成的方向顺序相同，但此方法产生了不必要的开销，因为必须重新计算连续迭代的第一象限，如图所示6（a） ●●●●。

图9。

或者，可以逐步构建公理生成的曲线的下一个展开步骤\（H）或\（A\），完全基于公理产生的当前扩展步骤\（A\）或\（H）分别为。这种见解也是命题的直接结果2和三，自

\（开始{方程式}H_{l+1}=A_l\\uparrow\\texttt{xor}_1（A_l）\右箭头\texttt{xor}_1（A_l）\\向下箭头\\texttt{xor}_2（A_l），结束{方程式}

(23)

和

\（开始｛方程｝A_｛l+1｝=H_l\rightarrow\texttt{xor}_1（H_l）\\uparrow\\texttt{xor}_1（H_l）\leftarrow\texttt{异或}_2（H_l）。\结束{方程式}\）

(24)

注意，如果我们的曲线是以这种方式构建的，那么迭代的第一象限\（l+1）由迭代中计算的曲线表示\（l），如图所示6（b），从而避免了以前的开销。然而，重要的是要理解Lindenmayer系统使用公理生成的曲线\（小时\）或\（A\）将始终在坐标对处结束其遍历\（（2^l-1,0）\）或\（（0,2^l-1）\）分别为。由于公理的第一个扩展步骤中存在不同的终端符号，我们方法的连续公理交换导致每个新曲线迭代的最后一对坐标不同\（H）和\（A\）。这些端子符号表示在给定迭代中用于连接四个象限的方向。从这一点开始，我们将这种类型的方向称为连接方向，并用变量表示\（_k\），其中\（k\in\l轨道1,2,3\r轨道\）表示这些端子在第一个扩展步骤中的外观顺序。提议4形式化连接方向的计算\（s1\）,\（s2\），以及\（s_3\）用于增量生成的曲线。

提案4。

连接方向变量保持的值\（s1\）,\（s2\），以及\（s_3\）用于迭代\（l+1\）可以通过将位互斥或乘以1应用于迭代的连接方向来计算\（l）.

证明。由于我们方法的增量性质，每当Hilbert填充曲线的奇数迭代由公理生成时\（H，\）最后访问的小区是\（（2^l-1,0）\）。如果我们要生成此曲线的偶数迭代，对于同一公理，要访问的最后一个单元格将是\（（0,2^l-1）\）相反。公理生成的曲线正好相反\（A\）。由于每条曲线都开始在网格单元上进行遍历\( (0,0) \)在进入下一个象限之前，必须完全遍历每个象限，该曲线的最后一个访问单元格为\（（0,2^l-1）\）必须首先访问右下象限，然后访问右上象限，完成对左上象限的遍历。这与连接方向相对应\（（（s1，s2，s3）=（\underrightarrow{\0\}_4，\underrightarrow{1\}_4、\under右箭头{2\}_5）。如果最后访问的单元格是\（（2^l-1,0）\），我们的曲线必须首先访问左上象限，然后访问右上象限并在右下象限完成其遍历。表示连接方向等于\（（（s1，s2，s3）=（\underrightarrow{\1\}_4，\underrightarrow{\0\}_4和\under右箭头{\3\}_5）.

独立于所选公理，在每一代迭代中，连接方向将在\（（\underrightarrow{\0\}_4，\under右箭头{1\}_4和\underrightarrow{2\}_5）和\（（\underrightarrow{\1\}_4、\underrightarrow{\0\}_4和\under右箭头{\3\}_3）.得出结论，可以计算迭代的连接方向\（l+1）通过对迭代的每个连接方向执行位向互斥或1\（l，\）

\（开始{方程式*}（s1，s2，s3）：=\大（\texttt{xor}_1（s_1），\texttt{xor}_1（s_2），\texttt{xor}_1（s_3）\大）。\方块\结束{方程式*}\）

理论上，命题4并不是特别相关。同样的曲线可以通过如果。。。其他的声明。然而，它提供了一个优雅的解决方案来计算\（_k\）变量，减少了该解决方案的分支指令数量，并有助于提高解决方案的整体性能和可读性。

3.3实施细节

提议2,三，以及4允许用C语言构建迭代希尔伯特填充空间曲线生成器，由于在每个迭代的构建步骤中缺乏数据依赖性，因此可以通过SSE指令进行扩展。有一些关于内存和尚未寻址的SSE指令数的关键细节。

如前所述，每个方向符号都可以用2位表示，但大多数类C语言不提供这样大小的基本类型。通过实现存储给定希尔伯特空间填充曲线方向的位阵列，可以避免内存浪费，从而减少SIMD指令的数量。我们的位数组实现为\（\texttt{uint8_t}\），一个8位的基元类型，意味着此数组的每个条目最多可以存储4个方向。给定该迭代\（l）希尔伯特填充曲线由\（4^l-1）方向，我们的比特阵列需要\（4^{l-1}\） \（\texttt{uint8_t}\）变量来存储完整的迭代。下面的方案说明了迭代是如何进行的\（l）Hilbert空间填充曲线的\（\texttt{曲线}\），其中变量\（_k\）代表\（k\）要处理的第个方向：

\（开始{等式}\texttt{uint8_t*curve}:=\underbrace{{\underleftarrow{\d_3\d_2\d_1\d_0\}_4}}{\texttt{curve{[0]}，\underlace{{underleft arrow}\d_7\d_6\d_5\d_4}_4{}}}_{\texttt曲线}[1]}，\\cdots\，\ underbrace{\undertleftarror{\d\d_4^{l}}}}\d_{4^{l} -1个}\d{4^{l} -2个}\d{4^{l} -3个}\}_4}{\texttt{曲线}[4^{l-1}]}。\结束｛方程式｝\）

(25)

方向可能看起来很奇怪\（_k\），在每个\（\texttt{uint8_t}\），以小发动机方式分配。此布局允许我们检索\（\texttt{uint8_t}\）在里面\（O（1）\）用余数计算时间(\( \% \))和位向右移位(\（\gt\gt\）)操作员。从比特阵列中检索方向\（\texttt{曲线}\）由以下嵌套循环描述：

需要注意的是，此布局修改了连接方向的位置\（s1\）,\（s2\），以及\（s_3\）在我们的比特阵列中。带迭代的Hilbert曲线的连接方向\（l）可以在该曲线的第一、第二和第三象限之后找到，如中所示

\（开始{方程式}H_l=A_{l-1}\下大括号{\\uparrow\}_{\显示样式s_1}H_{l-1}\下小括号{\\rightarrow}_{\displaystyles_2}H_}-l-1}\\underbrace{\\downarrow\}_{显示样式s_3}\B_{l-1\}，结束{方程式{）

(26)

和

\（开始{方程式}A_l=H_{l-1}\下大括号{\rightarrow}_{\displaystyle s_1}A{1-1}\上大括号{\\uparrow\}_{\显示样式s_2}A{l-1{\underbrace{\leftarrow{}显示样式s_3}C_{l-1}.\end{方程式{\）

(27)

隐含变量\（s1\）,\（s2\），以及\（s_3\）必须在中进行分配\（d_{4^{l-1}}\）,\（d_{2*4^{l-1}}\），以及\（d_{3*4^{l-1}}\）我们的位数组，可以在条目的两个最重要的位中找到\（\texttt{曲线}[4^{l-2}-1] \),\（\texttt｛curve｝[2*（4^｛l-2｝）-1]\），以及\（\texttt{曲线}[3*（4^{l-2}）-1]\）分别为。如果我们取消设置这些位，则可以通过位或操作在位数组中插入连接方向(\( | \)):

\（\开始{align}\texttt{curve}[4^｛l-2｝-1]&#x0026；：=\文本{曲线}[4^{l-2}-1]\|\s_1；\结束{align}\）

(28)

\（\开始{align}\texttt{curve}[（2*4^{l-2}）-1]&#x0026；：=\文本{曲线}[（2*4^{l-2}）-1]\|\s-2；\结束{align}\）

(29)

\（\开始{align}\texttt{curve}[（3*4^{l-2}）-1]&#x0026；：=\texttt｛curve｝[（3*4^｛l-2｝）-1]\|\s_3；\结束{align}\）

(30)

由于每个连接方向现在存储在\（\texttt{uint8_t}\）变量，则下一次曲线迭代的这些方向的计算与Proposition的计算略有不同4。在这种情况下，我们只需要应用位-wise exclusive-or(^{\（\楔形\）})从1到第二个最重要的连接方向。此操作可以通过使用二进制位掩码来实现\（\texttt{0b01000000}\）:

\（\begin{align}s_1:=s_1{\wedge}\\texttt{0b01000000}；\end{align{\）

(31)

\（\begin{align}s2:=s2{\wedge}\\texttt{0b01000000}；\end{align{\）

(32)

\（\begin{align}s_3:=s_3{\wedge}\\texttt{0b01000000}；\end{align{\）

(33)

算法1说明了我们的希尔伯特空间填充曲线生成器在C语言中的实现，并用SSE内蕴进行了扩展[14]. 为了保持伪代码的可读性，我们选择表示每个\（\texttt{uint8_t}\）以4为底的文字，明确描述此变量包含的方向，例如方向序列\（\右下箭头{\0\1\2\3}_4=\左下箭头{3\2\1\0}_4=228_{10}\）。我们现在将解释我们解决方案最重要的实施细节：

•

内存分配（第1行）：生成器通过为位数组分配内存来启动\（\texttt{曲线}\）基于\（\texttt{level}\）表示曲线的所需扩展级别的输入参数。如前所述，我们的比特阵列需要\（4^{\texttt{级别}-1} \) \（\texttt｛uint8_t｝\）存储Hilbert曲线方向的变量，水平等于\（\texttt{level}\）.

•

位阵列初始化（第2行至第11行）：取决于\（\texttt{公理}\）作为输入传递，我们的生成器初始化位数组\（\texttt{曲线}\）以表示希尔伯特填充曲线第三级展开的方向\（\texttt{公理}\）该曲线迭代按照方程中给出的方案进行分配(25). 请注意，具有此扩展级别的曲线完全填充SSE向量（128位），这允许我们避免与希尔伯特曲线的较小迭代相关的角点情况。

•

连接方向初始化（4至11号线）：在此阶段，连接方向\（s1\）,\（s2\），以及\（s_3\）也根据\（\texttt{公理}\）值。这些变量表示为\（\texttt{uint8_t}\）变量，如前所述。

•

专享或面具（第12和13行）：初始化两个SSE整数向量\（（\texttt{__m128i}）\）包含应用位向互斥或1所需的位掩码(\（\texttt{mask_xor1}\）)，或2(\（\texttt｛mask_xor2｝\）)，指向另一个SSE整数向量包含的所有方向。\（\texttt{mask_xor1}\）对应于具有八个\（\texttt{uint8_t}\）全部等于\（\左下箭头{%1\1\}_4\），同时\（\texttt{mask_xor2}\）对应于由八个\（\texttt{uint8_t}\）值等于\（\左下箭头{\2\2\2\\}_4\）.

•

生成下一个曲线迭代（第14行至第31行）：

–

象限边界的识别（第15至22行）：变量\（\texttt{offset_q2}\）,\（\texttt{offset_q3}\）和\（\texttt{offset_q4}\）存储第一个\（\texttt｛uint8_t｝\）位数组的条目\（\texttt{曲线}\）分别包含象限2、3和4的数据。

–

计算最后三个象限（第18至22行）：如提案所示2和三，可以生成迭代的最后三个象限\（l）通过应用\（\texttt{xor}_1 \)和\（\texttt{xor}_2 \)到这个迭代的第一个象限。由于此曲线是以增量方式构建的，因此迭代的第一象限\（l）与迭代中计算的曲线相同\（l-1\）这些函数可以由SSE内部函数编码\（\texttt{_mm_xor_si128}\），需要两个\（\texttt｛__m128i｝\）向量作为参数，并对这些向量的元素执行位排他或。在我们的应用程序中，输入向量是\（q_1\），对应于此迭代的第一象限所包含的数据块，以及\（\texttt｛mask_xor1｝\）或\（\texttt{mask_xor2}\），这取决于我们是要将exclusive-or by 1还是by 2应用于\（q_1\）。完成此位操作后，生成的向量将存储在相应的\（\texttt{曲线}\）通过使用先前计算的\（\texttt{offset_q2}\）,\（\texttt{offset_q3}\）和\（\texttt{offset_q4}\）最后，需要注意的是，重复此步骤直到所有\（\texttt{曲线}\）处理与第一象限相关的条目。

–

重置连接方向位（第23至25行）：由于功能的应用\（\texttt{xor}_1 \)和\（\texttt{xor}_2 \)，位数组中与连接方向相对应的位\（s1\）,\（s2\），以及\（s_3\）将不清楚。为了确保这些方向在我们的位数组中正确分配，我们必须首先重置相应的位。此过程在包含连接方向和位掩码的条目之间实现为位和（&）\（\texttt{0b00111111}\（\左下箭头{0\1\1}_4）\）.

–

分配连接方向（第26至28行）：一旦将分配连接方向的位重置，就可以通过位或(\( | \))如方程式所示(29)至(30).

–

设置下一次迭代（第29行到第31行）：此时，第一个\（4^{l+1}-1\）我们的比特阵列的比特对曲线包含表征迭代的方向\（l+1\）希尔伯特填充曲线。通过从第16行重复该过程，可以计算该曲线的进一步迭代。然而，由于发电机的增量特性，我们需要更新连接方向的值\（s1\）,\（s2\），以及\（s_3\）以下命题25和方程式(31)至(33).

3.4时间和记忆分析

提案5。

算法1时间复杂性为\（O（N）\），其中\（否）是网格条目数。

证明。

让\（l\in\mathbb{N}^+\）表示给定希尔伯特填充曲线的递归水平或迭代次数，以及\（a，b\in\mathbb{N}^+\）算法的外循环（第14行）和最内环（第18行）中的恒定时间操作数1分别为。计算水平曲线的操作总数\（l）对应于以下总和：

\（开始{聚集*}S_{l}=b（l-4）+\sum_{i=3}^{l-1}a（4^{i} -1个). \结束{聚集*}\）

总和\（S^{\素数}_{l}\）是一个几何级数和的上界\（_l\）

\（开始{聚集*}S^{素数}_l=b（l-4）+\sum_{i=3}^{l-1}a^{i} -1个)=S_{l}。\结束{聚集*}\）

自\（l=log _4（N）），这个几何级数可以进一步简化为

\（\begin{gather*}S^{prime}_l=b（\log_4（N）-4）+\frac{a（N-1）}{3}。\结束｛聚集*｝\）

自\（N\gt\log_4（N）），我们得出结论\（O（N）中的S^{\prime}_l\）。这意味着\（O（N）中的SL），鉴于此\（S^{\prime}_l\）是的上界\（_l\）.□

提案6。

算法1呈现的内存复杂性为\（O（N）\），其中\（否）是网格条目数。

证明。

算法1是迭代的，并且仅在第1行中的位阵列的初始化步骤期间分配动态内存。此位数组正好包含\（2N\）位，导致内存开销为\（O（N）\）.□

4希尔伯特和Z阶空间填充曲线的分块方案

为了将希尔伯特曲线推广到矩形网格，Böhm等人。[5]实现了一个解决方案，该解决方案使用了称为nano-programs的小型预处理曲线。每个纳米线代表一系列方向，形成一条曲线，能够填充边长从1到4的网格。然后通过细分，在希尔伯特曲线的任意两个连续方向之间应用正确的纳米线，完成该路径。纳米投影图的使用是一种巧妙的方法，可以填充矩形网格，同时保留子网格内的局部性。然而，我们怀疑有可能在不影响局部性的情况下，将纳米program简化为行主顺序扫描，如图所示7如果子网格足够小，并且宿主算法不存在空间局部性，则适用于非局部矩阵转置问题。提议7确保可以阻止任何\（n倍m）网格，用于\（\lfloor\log_2（n）\rfloor=\floor\log_2（m）\rffloor\）,\（n\gt 1）和\（m\gt 1），变成一个正方形\（2^{\lfloor\log_2（m）\rfloor-1}\乘以2^{\ lfloor\ log_2（m\rfloor-1}\）分区网格。此网格的尺寸允许按照希尔伯特或Z阶曲线施加的顺序遍历其块。

图11。

提案7。

可以拆分\（n倍m）网格，其中\（n\gt 1）,\（m\gt 1）和\（\lfloor\log_2（n）\rfloor=\floor\log_2（m）\rffloor\），到\（2^{\lfloor\log_2（n）\rfloor-1}\乘以2^{\floor\log_2（m）\rffloor-1}）边长从2到4的方块。

证明。

这个命题可以通过以下方程组来证明，其中\（\mathbf{length}\）是原始网格边的大小\（n\）或\（米\）.变量\（a，b\in\mathbb{N} _0（0） \)表示大小为2和3或3和4的段数，

\（\begin{方程式}2^{\lfloor\log_2（\mathbf{length}）\rfloor-1}=a+b\end{方程式{）

(34)

和

\（\begin{方程式}{\left\lbrace\begin{array}{ll}\mathbf{length}=2a+3b\quad\text{if}\；{\mathbf{length}\le\frac{3}{2}2^{\lfloor\log_2（\mathbf1{length-}）\rfloor}，}\\mathbf}length{3a+4b\quae\text{否则。}\end{arrary}\right.}\结束{方程式}\）

(35)

方程式(34)确保尺寸一侧的分段数\（\mathbf{length}\）等于\（2^{\lfloor\log_2（\mathbf{length}）\rfloor-1}\），是2的幂。While等式(35)确保所有线段的大小之和等于原始边长。对该系统应用几个替换步骤会导致

\（\begin{align}&#x0026；{\left\lbrace\begin{array}{ll}a=-\mathbf{length}+\frac{3}{2}\2^{\lfloor\log_2（\mathbf{length}）\rfloor}\\b=\mathbf-{length-}-2^{\lploor\log_2（\mathbf{longth}）\floor}\end{arrary}\right.}\文本{if}\；{\mathbf{length}\le\frac{3}{2}2^{\lfloor\log_2（\mathbf{length}）\rfloor}，}\end{align}\）

(36)

\（开始{align}&#x0026；{\left\lbrace\begin{array}{ll}a=-\mathbf{length}+2^{\lfloor\log_2（\mathbf{length}）\rfloor+1}\\b=\mathbf1{length-}-\frac{3}{2}2^{\lfloor\log_2（\mathbf{length}）\rfloor}\end{array}\right.}\文本{否则。}\结束{align}\）

(37)

变量\（a）和\（b），在方程式中(36)和(37)，为\（\mathbf{length}\gt 1），以及\（2^{\lfloor\log_2（\mathbf{length}）\rfloor}\）可被2整除，确保此系统始终具有以下整数解\（a）和\（b）.收尾侧\（n\）和\（米\），可以分为\（2^｛\lfloor\log_2（n）\lfloor-1｝\）和\（2^{\lfloor\log_2（m）\rfloor-1}\）尺寸介于2和4之间的段。□

算法2描述了我们将大小为2的平方填充曲线（例如希尔伯特和Z阶填充曲线）推广到\（n\次m\）网格，其中\（\lfloor\log_2（n）\rfloor=\floor\log_2（m）\rffloor\）。此阻塞策略的主要优点是它允许及时使用标准预取技术，这大大提高了异地换位问题的性能。我们解决方案最重要的细节是：

•

侧面拼接（第2行至第3行）：分块程序的第一步确定网格边是否拼接成长度为2和3或3和4的线段。如果电网侧尊重\（\mathbf{length}\le\frac{3}{2}2^{\lfloor\log_2（\mathbf{lendth}）\rfloor}\），然后将此侧拼接成长度为2和3的线段。如果这个不等式不成立，那么这一侧将拼接成长度为3和4的线段。每侧使用的分区类型由变量编码\（\texttt{n_type}\）和\（\texttt{m_type}\）。一旦确定了分区类型，就可以获得较小段的数量(\（\texttt{n0}\）和\（\texttt{m0}\）)和更大的部分(\（\texttt｛n_1｝\）和\（\texttt｛m_1｝\）)通过求解方程组(36)或(37)命题的7，取决于该维度的分区类型。

•

块坐标（第11行）：坐标对\（（\texttt{I}，\texttt}）标识分区网格中的单个块。这对坐标对应于\（小时）给定填充曲线的第步，其计算公式为\（\texttt{获取块坐标}\）该函数可以借助我们的SIMD Hilbert空间填充生成器或Böhm等人提出的非递归Lindenmayer系统来计算。[5]，用于基于希尔伯特曲线的遍历。为了创建基于Z顺序曲线的遍历，我们实现了libmorton中提供的内核[4]利用BMI2\（\texttt{_pext_u32}\）有效计算的指令\（（\texttt{I}，\texttt}）在恒定的时间和记忆中：

\（\开始{聚集*}\texttt{I}:=\texttt}_pext_u32}（h，\textt{0xAAAAAAAA}）；\\\文本{J}：=文本{_pext_u32}（h，文本{0x5555555}）；\结束{聚集*}\）

•

块边界（第12至13行）：变量\（\texttt{分钟}\）,\（\texttt{i_max}\）,\（\texttt{j_min}\），以及\（\texttt{j_max}\）表示给定网格块的边界。这些边界的值由以下公式计算\（\texttt{compute_boundaring}\）返回此块的边界，基于\（（\texttt{I}，\texttt}）以及分区网格的形状。通过以下内核，可以在恒定的时间和内存中获得给定块的边界：

\（开始{align*}\texttt{j_max}&#x0026；：=（\texttt{m_0}*\texttt}min}（\texttt{j}+1，\texttt}m_0}））+\文本{i_max}&#x0026；：=（\texttt{n_0}*\texttt}min}（\texttt{I}+1，\texttt}n_0}）+（\textt{n_1}*\text{max}（0，\texttt_I}+1-\text{n_0}））；\\\texttt｛j_min｝&#x0026；：=（\texttt{m_0}*\texttt}最小}（\texttt{J}，\texttt}m_0}）+（\textt{m_1}*\ttexttt{max}（0，\texttt[J}-\texttt（m_0}））；\\\文本{i_min}&#x0026；：=（\texttt{n_0}*\texttt}min}（\texttt{I}，\texttt}n_0}））+（\textt{n_1}*\ttexttt{max}（0，\texttte{I}-\texttt（n_0}））；\结束{align*}\）

•

预取（第14行）：在处理当前块的条目之前，我们的解决方案确定要处理的下一个块的边界，并预取与每个块行的第一个条目对应的缓存线。此步骤由编码\（\texttt{prefetch_next_block}\），雇佣GCC\（\texttt{__builtin_prefetch}（\texttt{addr}，\texttt}rw}，texttt{locality}）[21]. 参数\（\texttt{addr}\）表示要带到内存中的内存地址。\（\texttt{rw}\）是一个编译时间常数，指示预取数据是用于读取（0）还是写入（1）操作。\（\texttt｛locality｝\）也是一个编译时间常数，范围从0到3，其中0表示预取数据没有时间局部性，此缓存线在访问后将被逐出，而3表示预取的数据具有最高程度的时间局部性并将在所有缓存级别中持久存在。

•

块处理（第15行）：预取步骤结束后，我们的解决方案将根据先前分配的变量开始处理当前块的条目\（\texttt{i_min_prev}\）,\（\texttt{i_max_prev}\）,\（\texttt{j_min_prev}\），以及\（\texttt{j_max_prev}\）。然后按行主要顺序处理此块。

4.1严重不对称网格

提议7对于严重不对称的网格，未进行验证，即。，\（\lfloor\log _2（n）\rfloor\ne\lfloor/log _2（m）\rffloor\）。用于穿过这些尺寸的网格的策略与Böhm等人所采用的策略类似。[5]. 我们的目标是拼接原始\（n倍m）将网格划分为多个\（n^{\prime}\乘以m^{\prime}\）子网格，例如\（\lfloor\log_2（n^{\prime}）\rfloor=\lfloor \log_2.主张7为每个人保留\（n^{\prime}\乘以m^{\prime}\）子网格，意味着可以通过迭代遍历所有子网格\（\texttt｛min｝（n，m）-1\）希尔伯特曲线或Z阶曲线。提议8确保在以下情况下可以进行网格拆分\（\lfloor\log _2（n）\rfloor\ne\lfloor/log _2（m）\rffloor\）,\（n \ gt 1 \）和\（m\gt 1）.

提案8。

任何\（n倍m）网格，其中\（\lfloor\log _2（n）\rfloor\ne\lfloor/log _2（m）\rffloor\），可以拆分为\（{\big\lfloor\frac{\texttt{max}（n，m）}{2^{\lfloor \log_2（\texttt{min}（n，m））\rfloor}}\big\floor}\）子网格，其尺寸符合\（\lfloor\log_2（\texttt{min（n，m））}\rfloor=\lfloor \log_2.

证明。

假设初始值\（n倍m）网格，其中\（n\gt 1）,\（m\gt 1）和\（日志2（n）。然后可以重写\（米\）作为

\（开始{方程式}m=\bigg（\bigg\lfloor\frac{m}{2^{\floor\log_2（n）\rfloor}}\bigg\floor-1-bigg）2^{\lfloor \log_2。\结束{方程式}\）

(38)

方程式(38)建议有可能将这个水平严重不对称的网格拆分为一个子网格\（n\） \（\次\） \（（2^{\lfloor\log_2（n）\rfloor}+m\bmod 2^{\ lfloor\ log_2（n）\rffloor}）\）子电网，加上\（{地板m/2^{地板\log_2（n）地板}地板-1}）带尺寸标注的子网格\（n\） \（\次\） \（2^{\lfloor\log_2（n）\rfloor}\）.自\（\lfloor\log_2（n）\rfloor=\lfloor \log_2，命题7保留所有子网格。□

提议8意味着可以通过使用算法将我们的阻塞过程扩展到严重不对称的网格2到每个\（n^{\prime}\乘以m^{\prime}\）网格，同时更新变量值\（\texttt{分钟}\）和\（\texttt{j_min}\）相应地。

4.2电网间位置

如第节所述3.2，公理生成的希尔伯特填充曲线访问的最后一个单元格\（H）将永远是\（（2^l-1,0）\），而公理生成的曲线访问的最后一个单元格\（A\）是\（（0,2^l-1）\）。这一事实使我们能够利用曲线创建过程来改善网格间的局部性，以防我们想要处理严重不对称的网格。为了保持网格间局部性，我们的策略必须确保给定子网格的最后两个象限与下一个子网格的前两个象素相邻。满足此限制的第一步是确定网格是否水平，\（\log _2（n）\） \（\ lt\） \（\log _2（m）\），或垂直，\（\日志2（n）\） \（\gt\） \（\log _2（m）\），严重不对称。如果我们的网格是水平不对称的，那么每个子网格的最后一个访问块将是\（（k（2^l-1），0）\），其中\（k\in\mathbb{N}\）代表当前子电网。如果我们的电网垂直严重不对称，\（（0，k（2^l-1））\）将是最后访问的子网格块。根据命题的观察4，如果我们的曲线在水平方向上严重不对称，那么我们的曲线必须由公理生成\（H）或\（A\）对于奇数或偶数值\（l）分别为。然而，如果我们的曲线是垂直不对称的，那么它必须由公理生成\（A\）或\（H）对于奇数或偶数值\（l）最后，需要注意的是，基于Z阶曲线的遍历总是在与其起点对角的单元格中结束，这意味着不可能利用该曲线的生成来提高网格间的局部性。

4.3时间和记忆分析

提案9。

我们的分块方案、算法2，个礼物\（O（nm）\）时间复杂性，其中\（n\）和\（米\）对应于其输入网格的行数和列数。

证明。

原件\（n倍m）网格按功能划分为非重叠块拆分_侧面，哪些成本\（O（1）\）时间，并且对于生成的每个希尔伯特空间填充曲线调用一次。每个块的边界由以下公式计算获取块坐标，这也需要花费\（O（1）\）时间。最后，函数预取新块发出与给定块的行数成比例的软件预取指令数。因此，算法的复杂性2与输入网格大小成比例\（O（纳米）\）.□

提案10。

我们的分块方案、算法2，个礼物\（O（纳米）\）SIMD Hilbert空间填充曲线生成器非递归Lindenmayer系统生成的曲线的内存复杂性[5]，或\（\texttt{_pext_u32}\）功能。

证明。

独立于所选的生成器，算法2假设存在\（n倍m）内存中的网格。如果我们的阻塞方案使用的曲线是由非递归Lindenmayer系统生成的，或者由\（\texttt{_pext_u32}\），则不会动态分配更多内存。因此，算法的内存复杂性2是\（O（纳米）\）如果该算法使用的曲线是由我们的SIMD Hilbert空间填充曲线生成器生成的，那么算法2增加了每个块2位的开销。由于任何网格的块数都小于条目总数，因此该算法的内存复杂度仍然等于\（O（纳米）\）.□

5实验评估

以下实验在Intel Xeon Gold 6130（Skylake）CPU中进行，具有2.1 GHz固定时钟频率、32 KiB（8路）专用一级数据缓存、1 MiB（16路）专用二级缓存和22 MiB（11路）共享三级缓存。所有缓存级别都采用写回策略。该机器包含6个内存通道，最高内存速度为2.6 GHz。此计算机中安装的操作系统是CentOS Linux 7，版本3.10.0-1127.19.1.el7.x86_64。所有实验都是用C编写的，并使用GNU GCC Red Hat 9.4.0进行编译。使用性能版本3.10.0评估内存指标。

5.1希尔伯特曲线生成器

SIMD Hilbert空间填充曲线生成器的三个版本，算法1，以评估使用SSE或AVX512指令与编译器自动矢量化实现之间的性能差异。将这些解与线性时间内枚举给定Hilbert空间填充曲线所有点的已知方法进行了比较，例如非递归和递归Lindenmayer系统和FUR-Hilbert。

数字8（a）以及8（b）描述了每种方法每秒处理的希尔伯特填充曲线条目数的总共15次运行的平均值，用\（\texttt{-O1-march=native}\）和\（\texttt{-O3-march=native}\）分别为。值得注意的是，每次运行都会计算出完整的希尔伯特曲线。此分析没有利用算法的可能摊销1此外，每个测试都会遍历所有\（（i，j）\）对应于希尔伯特空间填充曲线给定迭代的点的坐标线。

图12。

5.1.1编译器优化-O1公司.

优化级别的最佳方法\（\texttt{-O1}\）是SIMD希尔伯特的AVX512实现，它提供了\（6.70\cdot 10^8\）每秒条数，平均值为\（7.32\cdot 10^8）从10级开始每秒输入个数。通过求助于SSE指令来实现此解决方案的时间很短，平均处理时间很短\( 6.47\% \)每秒条目数少于AVX512。我们的解决方案的自动矢量化版本是\( 9.91\% \)比最后五个级别的SSE版本慢。从6级开始，所有SIMD-Hilbert实现都优于FUR-Hilbert。我们的AVX512和SSE实施平均水平\（1.53倍）和\（1.41\次）在最后五个级别上比FUR Hilbert更快。

5.1.2编译器优化-臭氧.

优化标志\（\texttt{-O3}\）显著减少了不同SIMD Hilbert Space-filling Generator版本的运行时间。从编译器优化中受益最多的方法是我们方法的自动矢量化实现\（5.70\次）加速，及其SSE对应项，其表现出类似的加速\（5.37\次\）。编译器优化的AVX512版本的生成器提供了较低的性能增益，加速比为\（3.58次）。编译器无法有效地优化其余方法。FUR-Hilbert的加速比\（1.04倍），而Lindenmayer系统的递归和非递归版本的速度提高了\（1.02倍）和\（1.27倍）此优化级别最具性能的方法是我们的SIMD Hilbert Space-filling Curve Generator的SSE和自动矢量化实现。从级别10开始，SSE版本计算的平均值为\（3.68\cdot 10^9）每秒曲线点数，而自动矢量化版本的平均值为\（3.51\cdot 10^9）每秒点数。平均而言，SSE实现计算\( 4.55\% \)每秒比自动矢量化版本的曲线点数更多。AVX512版本的发电机是性能第三高的方法\( 71.26\% \)和\( 74.51\% \)SSE和自动矢量化实现分别提供的吞吐量。算法的不同实现1与FUR-Hilbert和两个Lindenmayer-System版本相比，经过编译器优化后，性能有了很大提高。最值得注意的是，我们的SIMD Hilbert空间填充曲线生成器的SSE实现，它是\（6.90\次）比FUR-Hilbert更快，\（27.63\次）和\（14.32\次）比非递归和递归Lindenmayer-Systems更快。

5.2错位换位

本节评估了五种不同方法的性能，这些方法用于计算矩形单精度矩阵的外置转置。被评估的方法包括一种朴素的换位策略[19]; FUR-Hilbert的缓存自适应[5];\（\texttt{somat_copy}\），“英特尔MKL”提供的缓存软件手动调整换位策略[18]; 以及我们的两个版本的缓存缓冲块方案，该方案使用基于Hilbert或Z阶曲线的遍历。除非明确说明，否则我们的希尔伯特阻塞方案对应于我们的阻塞方案的一个版本，其曲线由我们的SIMD希尔伯特填充空间曲线生成器的SSE版本生成。为了进一步优化这些实现，编译标志\（\texttt{-O3}\）和\（\texttt{-march=native}\）已应用。

5.2.1 SIMD并行性。

值得注意的是，“英特尔MKL”提供的解决方案利用了SIMD并行性。为了获得竞争性结果，使用SSE内禀转置核对其余方法进行了扩展[14]执行\（4乘以4）子矩阵。此内核的使用意味着每个访问的网格条目都对应于此大小的子矩阵。实际上，这些实现只能计算大小是4的倍数的矩形矩阵的转置。需要使用内存填充来处理不同大小的矩阵。

5.2.2软件预取。

通过我们的分块方案实现软件预取是通过GCC构建完成的\（\texttt{__builtin_prefetch}\），位置提示等于3。我们没有注意到不同位置值的执行时间之间有任何显著差异。尽管“英特尔MKL”的实现是封闭源代码的，但我们希望此库能够利用软件预取技术。

5.2.3时间基准。

通过使用由\（\texttt｛time.h｝\）库。该测量不考虑初始化输入和输出矩阵所需的时间，而是五次独立运行的平均值。每种方法的性能都表示为吞吐量（GiB/s），计算如下

\（\[\texttt{吞吐量}=\frac{3\cdot\texttt{matrix_elements}\cdot\t{sizeof（float）}}{2^{30}\cdot \texttt}time}}.\]\）

溪流[17]通过计算执行期间加载和存储的内存量来计算基准的吞吐量。在不合时宜的转置上下文中，我们对每个转置的条目进行2次读取和1次存储，因此公式中的常量乘数。此公式可以直接与STREAM基准套件获得的结果进行比较，更具体地说\（\texttt{COPY}\）基准。\（\texttt｛COPY｝\）为大型阵列的现成复制过程提供了给定系统的实际峰值内存吞吐量。该基准不发出算术指令，提供矩阵转置过程可实现的最大吞吐量的上限。这个\（\texttt{COPY}\）在这个实验设置上，基准测试返回10.6 GiB/s的吞吐量。该系统上给定方法的效率计算如下

\（\[\frac{\texttt{approach_throughput}}{\texttt{COPY_throught}}.\]\）

5.2.4记忆基准。

可以通过评估PMU计数器来描述全局内存层次结构访问行为\（\texttt{cycle_activity.stalls_mem_any:u}\），统计内存加载指令挂起时的执行暂停数，以及\（\texttt{resource_stalls.sb:u}\），它统计存储缓冲区已满时暂停的周期数。后缀\（\texttt{：u}\）表示这些事件是在用户空间中测量的。这些事件可能同时发生，这意味着两个事件可能会计算单个暂停。尽管如此，这些计数器的总和应该能够很好地估计与内存相关的暂停。这些值用于计算与朴素换位方法相比，缓存友好方法的内存相关暂停减少量。该值计算如下

\（\[\frac{\texttt{naive_stalls}}{\texttt{approach_stalls{}.\]\）

5.2.5矩阵尺寸。

我们借用了弗拉迪米洛夫[23]将性能分析分为好矩阵和丑矩阵。良好的矩阵预计会在冲突未命中次数较少时发生，这意味着大多数内存暂停都是冷内存和容量未命中的结果。丑陋的矩阵会引发大量冲突未命中，将与矩阵数据对应的内存地址映射到数量减少的缓存集，这会导致很大一部分缓存未使用。丑陋的矩阵对应于边是关键步幅倍数的矩阵。这一步幅，解释如下雾[10]，可以计算为

\（\[\texttt{critical_stride}=\frac{\texttt{cache_size}}{N\texttt{-ways}\cdot\texttt}sizeof（TYPE）}.\]\）

此Intel Xeon Gold 6130上的单精度浮点大小为4 Bytes，这导致缓存L1d的关键步幅为1024个元素，缓存L2的关键步长为16384个元素，最后一级缓存的关键步速为524288个元素。自然，对于与两个或更多缓存关键步幅的倍数相对应的大小，冲突错过的数量会进一步增加。为以下实验选择的矩阵由单精度浮点数组成，当前边长在8192和16384之间。将高非方矩阵和宽非方矩阵的实验平均起来，以保持评估简洁。非矩形矩阵的最大维数在所有测试中保持不变。最大的维度包含16100和16834个元素，用于表示良好和丑陋的矩阵。请注意，测试的矩阵大小都不适合缓存。

5.2.6评估：良好的平方矩阵。

数字9（a）以及9（b）描述了用于大小良好的正方形矩阵的测试方法的吞吐量和与内存相关的停滞的减少。正如预期的那样，与原始的cache-ignorant解决方案相比，高效缓存方法提供了显著的性能提升。令人惊讶的是，与“英特尔MKL”缓存软件解决方案相比，FUR-Hilbert和我们的Hilbert曲线阻塞方案都取得了显著进步。我们的希尔伯特曲线分块方案通过标准预取技术进行扩展，表现出最佳性能，平均吞吐量为6.75 GiB/s（效率为64%）\（3.17倍）与内存相关的暂停减少。将标准预取技术应用于该方法，平均加速率为31%。

图13。

我们的Z阶分块方案也使用相同的预取技术进行了扩展，对于维数小于12的矩阵，它表现出了第二好的性能，\（900次12900次），平均吞吐量为6.15 GiB/s。对于较大的矩阵，该方法的性能稳步下降，对于维数大于12的矩阵，其性能优于FUR-Hilbert，\（100乘以12100）对于尺寸大于14的矩阵，使用MKL，\（500乘以14500）.

将软件预取应用于该方法，平均加速率为24%。总的来说，性能第二好的方法是FUR-Hilbert，平均吞吐量稳定在5.85 GiB/s。最后，值得一提的是，“英特尔MKL”的执行时间与暂停减少似乎没有关联。“英特尔MKL”的封闭源代码特性意味着无法理解此方法是如何实现的，以及随后应分析哪些PMU计数器。

5.2.7评估：丑陋的方阵。

数字10（a）以及10（b）描述了针对丑陋方阵的测试方法的吞吐量和与内存相关的暂停减少。此分析中出现的每个矩阵大小都是缓存L1关键步幅的倍数。值得注意的是，边长等于16384的矩阵也是缓存L2关键步幅的倍数，这进一步增加了冲突未命中的数量。虽然可以观察到每种方法的性能下降，但与第节中的结果相比5.2.6对于步幅是L1和L2关键步幅的倍数的矩阵，“英特尔MKL”提供的解决方案是唯一一个性能严重下降的解决方案。对于丑陋矩阵，性能最好的方法是我们的希尔伯特曲线阻塞方案，该方案再次使用标准预取技术进行了扩展，平均吞吐量为6.15GiB/s（效率为58%），失速减少为\（1.75\次\），对于维数在9之间的矩阵，\（216\乘以9216\）和15，\（360次15360次）。这种方法紧随其后的是我们的Z顺序阻塞方案，该方案的平均吞吐量为5.96 GiB/s，失速减少了\（1.78\次\）应用软件预取丑陋平方矩阵，与希尔伯特和Z阶分块方案的非重取版本相比，速度分别提高了28%和26%。我们认为，由于纳米投影图的使用，FUR-Hilbert将优于我们的两种方法，这将导致较小的垂直扫描序列，从而减少冲突遗漏的数量。然而，FUR-Hilbert的平均吞吐量为5.31 GiB/s，低于其余的缓存策略。MKL表现出最低的性能增益，平均吞吐量为5.07 GiB/s，最低的失速减少平均值为\（1.45\次\）.

图14。

5.2.8评估：良好的矩形矩阵。

数字11（a）以及11（b）描述了针对大小良好的非正方形矩形矩阵的测试方法的吞吐量和与内存相关的停滞减少。这类矩阵的结果与大尺寸方阵的结果几乎相同。可以得出这样的结论：这些方法的吞吐量对输入矩阵形状几乎没有影响。我们的Hilbert阻塞方案通过标准预取技术进行扩展，具有最佳吞吐量，平均6.47 GiB/s（效率为61%），最大失速减少\（3.08\次\）将软件预取应用于我们的Hilbert阻塞方案，与非预取方案相比，平均吞吐量增加了19%。我们的方法紧随其后的是FUR-Hilbert，平均吞吐量为5.98 GiB/s，平均失速减少了\（2.22倍）。其余方法均未与基于希尔伯特的阻塞方案或FUR-Hilbert的阻塞方案产生竞争性结果。

图15。

5.2.9评估：丑陋的矩形矩阵。

数字12（a）以及12（b）描述了针对丑陋的非方形矩形矩阵的测试方法的吞吐量和与内存相关的暂停减少。与评估大小合适的矩阵相比，大小丑陋的矩阵在平方输入矩阵和非平方输入矩阵之间的吞吐量发生了巨大变化。对于这种类型的矩阵，MKL始终优于所有其他方法，平均吞吐量为4.91 GiB/s（效率为46%）。不幸的是，这种方法所带来的失速减少似乎与其吞吐量无关。FUR-Hilbert是唯一一种不受不同矩阵形状严重影响的缓存缓冲方法，与“英特尔MKL”解决方案相比具有竞争力。该策略的平均吞吐量为4.83 GiB/s，平均失速减少为\（2.02\次\）。我们的两种方法都表现出较差的性能，以至于应用标准预取技术导致性能下降。我们的希尔伯特（Hilbert）和Z阶分块方案通过软件预取进行了扩展，其开销为\( 9\% \)和\( 2\% \)与它们的非反射对应物相比。对于丑陋的矩形矩阵，软件预取的使用似乎会导致有用数据被驱逐，转而使用当前不需要的数据。目前尚不清楚为什么这种现象只发生在丑陋的非方阵上。

图16。

5.2.10 SIMD Hilbert生成器与非递归Lindenmayer系统。

数字13（a）以及13（b）描述了我们的SIMD Hilbert曲线生成器生成的Hilbert图的分块方案和böhm等人提出的非递归Lindenmayer系统的吞吐量，包括好的和不好的平方矩阵。[5]. 这两种方法都通过软件预取指令进行了扩展。这些方法对这两种类型的矩阵表现出相同的行为。对于良好和丑陋的矩阵，非递归Lindenmayer-System分别比使用我们的生成器的相同分块方案慢2%和不到1%。

图17。

很明显，用于生成希尔伯特填充曲线的方法不会显著影响我们的分块方案的性能。回想一下，我们的SIMD Hilbert空间填充曲线生成器的平均值\（27.63\次）使用编译器标志进行优化时，比非递归Lindenmayer-System更快\（\texttt{-O3}\）.

值得注意的是，尽管我们的生成器需要一个辅助缓冲区，因此需要更多的内存移动，但它仍然为这个内存绑定问题提供了最佳性能。然而，如果存在内存限制，那么我们的阻塞方案的非递归Lindenmayer系统版本可能是最合适的方法。

6最后备注

在这项工作中，我们介绍了一种用于矩形矩阵外置换位的高性能缓存转发算法。我们的阻塞方案利用了空间填充曲线的局部保持特性，能够填充\（2^l乘以2^l）方形网格和软件预取技术，以便携的方式改进此问题的执行时间。我们的实验表明，对于不合时宜的矩形换位问题，优化的缓存可发布解决方案能够优于当前手工调整的缓存软件最新解决方案。

这些结果是由于我们的算法具有优越的缓存性能才得以实现的。不幸的是，第节的内存度量所显示的结果5.2.4并不总是反映在不同方法的行为中。这个问题源于现代CPU中执行暂停的定义不明确。当代处理器能够在每个周期内发出多个微指令，因此不清楚停滞的周期是指没有发出微指令的周期，还是指发送的微指令数量低于处理单元最大能力的周期。软件预取技术的使用进一步增加了评估不同方法的复杂性，因为非临时预取指令绕过缓存，并且不被缓存相关PMU计数。

我们的块方案的加速表明，我们的算法在共享内存并行环境中表现良好。在这种环境下，与曲线生成相关的成本可能可以通过计算每个处理单元重新使用的较小曲线来摊销。最后，值得注意的是，我们的生成器和阻塞方案为其他高性能缓存可利用解决方案打开了大门，例如一般的稠密线性代数算法，其固有的内存长度较差，并且可以从软件预取中受益。

尽管我们的SIMD希尔伯特填充发生器在矩阵转置的情况下没有表现出显著的性能提升，但与其他已知的枚举希尔伯特填充曲线点的方法相比，它仍然表现出了显著的加速。

致谢

我们要感谢努诺·罗玛教授和佩德罗·托马斯教授，他们是INESC-ID的HPCAS研究小组的成员，他们借用了本研究初步阶段使用的机器。

脚注

在C语言中，如果左侧项为正，则余数运算符等价于模运算。如果左手端项为负，余数运算符的结果等效于计算模运算的对称性。

工具书类

[1]

吉尔赫梅·阿罗兹、何塞·蒙泰罗和阿林多·奥利维拉。2018计算机体系结构《世界科学》。内政部：

交叉参考

谷歌学者

[2]

迈克尔·巴德。2012空间填充曲线：科学计算应用导论第9卷。施普林格科技与商业媒体。

谷歌学者

[3]

迈克尔·巴德（Michael Bader）和克里斯托夫·曾格（Christoph Zenger）。2006.使用基于Peano曲线的元素排序缓存不经意的矩阵乘法。林·阿尔盖布。申请。417, 2–3 (2006), 301–313.

交叉参考

谷歌学者

[4]

J.Baert。2018.Libmorton:C++Morton编码/解码库。检索自https://github.com/Forceflow/libmorton.

谷歌学者

[5]

克里斯蒂安·博姆（Christian Bohm）、马丁·佩达切尔（Martin Perdacher）和克劳迪娅·普兰特（Claudia Plant）。2016年，基于新型填充曲线的Cache-obliovious循环。在IEEE国际大数据会议（大数据）IEEE标准。内政部：

交叉参考

谷歌学者

[6]

A.R.布茨。1971.希尔伯特填充曲线的替代算法。IEEE传输。计算。C-20（1971），424–426。

数字图书馆

谷歌学者

[7]

S.Chatterjee和S.Sen.2000。高效缓存矩阵转置。在第六届高性能计算机体系结构国际研讨会. 195–205.内政部：

交叉参考

谷歌学者

[8]

陈宁涛，王能超，石宝昌。2007.一种新的希尔伯特顺序编码和解码算法。软件：实际。专家。37, 8 (2007), 897–908.内政部：

交叉参考

谷歌学者

[9]

N.乔姆斯基。1956.语言描述的三种模式。IRE事务处理。信息Theor。2, 3 (1956), 113–124.

交叉参考

谷歌学者

[10]

阿格纳·福格。2020年。在C++中优化软件。检索自https://www.agner.org/optimizing/optimiting_cpp.pdf.

谷歌学者

[11]

马特奥·弗里戈（Matteo Frigo）、查尔斯·雷瑟森（Charles E.Leiserson）、哈拉尔德·普罗科普（Harald Prokop）和斯里达尔·拉马钱德兰（Sridhar Ramachandran）。2012.缓存可用算法。ACM事务处理。阿尔戈。2012年1月8日、1日、1–22日。内政部：

数字图书馆

谷歌学者

[12]

Kazushige Goto公司。2007年，转到BLAS。检索自http://www.tacc.utexas.edu/resources/software/.

谷歌学者

[13]

大卫·希尔伯特。1935年，在弗拉琴斯克，你将成为阿比尔顿·埃内尔·利尼。在漂移带：分析·Grundlagen der Mathematik·Physik Verschiedens.弹簧，1–2。

谷歌学者

[14]

英特尔&#xAE；内部函数指南。英特尔&#xAE；内部函数指南。检索自https://software.intel.com/sites/landingpage/Intrinsic指南/.

谷歌学者

[15]

金国华和约翰·梅勒·克鲁米。2005.填充曲线生成：基于表格的方法。算法数学和计算机科学国际会议. 40–46.

谷歌学者

[16]

K.Patrick Lorton和David S.Wise。2007.分析Morton-order和Morton-hybrid矩阵中的块局部性。ACM SIGARCH计算。阿基特。新闻35, 4 (2007), 6–12.

数字图书馆

谷歌学者

[17]

约翰·麦卡宾（John D.McCalpin）。1991-2007.STREAM：高性能计算机中的可持续内存带宽技术报告。弗吉尼亚大学，夏洛茨维尔，弗吉尼亚州。检索自网址：http://www.cs.virginia.edu/stream/.

谷歌学者

[18]

MKL 2020。英特尔&#xAE；数学内核库。检索自https://software.intel.com/content/www/us/en/develop/tools/math-kernel-library.html.

谷歌学者

[19]

Michael Mol.2007年。矩阵转置。检索自https://rosettacode.org/wiki/Matrix_transposition#C.

谷歌学者

[20]

G.M.莫顿。1966.面向计算机的大地测量数据库和文件排序新技术（加拿大渥太华：IBM）。(1966).

谷歌学者

[21]

天鹅座。2013年GCC提供的其他内置功能。检索自https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html.

谷歌学者

[22]

Przemyslaw Prusinkiewicz和James Hanan。1989林登迈耶系统、分形和植物纽约施普林格。内政部：

交叉参考

谷歌学者

[23]

安德烈·弗拉迪米洛夫。2013.使用Intel Xeon处理器和Intel Xeon-Phi协处理器的通用代码对平方矩阵进行多线程换位。https://www.researchgate.net/publication/258048110_Multithreaded_Transposition_of_Square_Matrices_with_Common_Code_for_Intel_Xeon_Processors_and_Intel_Seon_Phi_Coprocessors.

谷歌学者

[24]

张贤义、王倩、扎赫尔·乔提亚。2012年，OpenBLAS。检索自http://xianyi.github.io/OpenBLAS.

谷歌学者

引用人

查看全部

斯沃曼S瓦班斯库A皮门特尔A萨尔茨堡A克拉什纳霍尔凯A巴尔萨莫S纽顿带W阿巴德C尚伟(2024)用进化算法求阵列的无缓存广义Morton布局第十五届ACM/SPEC性能工程国际会议记录10.1145/3629526.3645034(83-94)在线发布日期：2024年5月7日
https://dl.acm.org/doi/10.1145/3629526.3645034

索引术语

基于高速缓存不经意Hilbert曲线的矩阵换位分块方案
1. 计算数学
  1. 数学软件
    1. 数学软件性能
2. 计算理论
  1. 算法的设计和分析
    1. 算法设计技术
  2. 计算模型
    1. 抽象机器

建议

精确TU分解的Cache-oblievious矩阵乘法
PASCO 2017：并行符号计算国际研讨会会议记录

基于莫顿混合空间填充曲线表示，我们提出了一种矩阵乘法的缓存自适应方法，用于有限域上矩形矩阵的并行TU分解。为了实现这一点，我们。。。
阅读更多信息
低深度缓存可用算法
SPAA’10：第二十二届年度ACM算法和架构并行性研讨会会议记录

在本文中，我们探索了一种简单而通用的方法来开发并行算法，从而在具有专用或共享缓存的并行机上实现良好的缓存复杂性。方法是设计深度较低（……）的嵌套并行算法。。。
阅读更多信息
缓存感知程序和缓存感知程序的实验比较
2007年SPAA：第十九届ACM并行算法和架构年会会议记录

为了利用现代微处理器的内存层次结构，可以绕过优化应用程序的一些困难，已经提出了缓存缓冲算法。这些算法基于分治。。。
阅读更多信息

信息和贡献者

问询处

发布于

ACM数学软件汇刊第48卷第4期

2022年12月

339页

ISSN公司：0098-3500

EISSN公司：1557-7295

内政部：10.1145/3572845

编辑：
白昭君
美国加州大学戴维斯分校
,
沃尔夫冈·班格尔
美国科罗拉多州立大学

期刊目录

本作品根据Creative Commons Attribution International 4.0许可证授权。

出版商

计算机协会

美国纽约州纽约市

出版历史

出版：2022年12月19日

在线AM：2022年8月9日

认可的：2022年6月27日

修订过的：2022年5月25日

收到：2021年1月7日

发表于TOMS体积48,问题4

权限

请求对此文章的权限。

请求权限

检查更新

作者标记

限定符

研究文章
推荐

资金来源

技术基金会（FCT）
奥地利科学基金

贡献者

其他指标

查看文章指标

文献计量学和引文

文献计量学

文章指标

1
引文总数
查看引文
1,364
总下载次数

下载量（最近12个月）807
下载次数（最近6周）36

其他指标

查看作者指标

引文

引用人

查看全部

斯沃特曼S瓦班斯库A皮门特尔A萨尔茨堡A克拉什纳霍尔凯ABalsamo S公司纽顿带W阿巴德C尚伟(2024)使用进化算法为阵列寻找缓存友好的广义Morton布局第十五届ACM/SPEC性能工程国际会议记录10.1145/3629526.3645034(83-94)在线发布日期：2024年5月7日
https://dl.acm.org/doi/10.1145/3629526.3645034

视图选项

查看选项

PDF格式

以PDF文件查看或下载。

PDF格式

电子阅读器

使用联机查看电子阅读器.

电子阅读器

HTML格式格式

在中查看本文HTML格式格式。

HTML格式

获取访问权限

登录选项

请检查您是否可以通过登录凭据或您的机构访问此文章以获得完全访问权限。

摘要

1引言

1.1缓存可用算法

1.2缓存式异地换位

1.3贡献和文件组织

计算希尔伯特填充曲线的两种已知方法

三SIMD（单指令多数据）希尔伯特填充曲线生成器

3.1计算曲线象限

3.2连接曲线象限

3.3实施细节

3.4时间和记忆分析

4希尔伯特和Z阶空间填充曲线的分块方案

4.1严重不对称网格

4.2电网间位置

4.3时间和记忆分析

5实验评估

5.1希尔伯特曲线生成器

5.1.1编译器优化-O1公司.

5.1.2编译器优化-臭氧.

5.2错位换位

5.2.1 SIMD并行性。

5.2.2软件预取。

5.2.3时间基准。

5.2.4记忆基准。

5.2.5矩阵尺寸。

5.2.6评估：良好的平方矩阵。

5.2.7评估：丑陋的方阵。

5.2.8评估：良好的矩形矩阵。

5.2.9评估：丑陋的矩形矩阵。

5.2.10 SIMD Hilbert生成器与非递归Lindenmayer系统。

6最后备注

致谢

脚注

工具书类

引用人

索引术语

建议

精确TU分解的Cache-oblievious矩阵乘法

低深度缓存可用算法

缓存感知程序和缓存感知程序的实验比较

评论

问询处

发布于

出版商

出版历史

权限

检查更新

作者标记

限定符

资金来源

贡献者

其他指标

文献计量学

文章指标

其他指标

引文

引用人

查看选项

PDF格式

电子阅读器

HTML格式格式

获取访问权限

登录选项

完全访问权限

数字

其他

分享

共享此出版物链接

在社交媒体上分享

附属公司