c - SIMD 2D矩阵英特尔指令集

Question

我正在开发基于Intel指令集（AVX，FMA等）的高性能算法。当数据按顺序存储时，我的算法（内核）运行良好。但是，现在我面临一个大问题，但没有找到解决方法或解决方案：
see 2D Matrix

int x, y; x = y = 4096;
float data[x*y]__attribute__((aligned(32)));
float buffer[y]__attribute__((aligned(32)));

/* simple test data */ 
for (i = 0; i < x; i++)
    for (j = 0; j < y; j++)
        data[y*i+j] = y*i+j; // 0,1,2,3...4095, | 4096,4097, ... 8191 |...

/* 1) Extract the columns out of matrix */
__m256i vindex; __m256 vec;
    vindex = _mm256_set_epi32(7*y, 6*y, 5*y, 4*y, 3*y, 2*y, y, 0);


 for(i = 0; i < x; i+=8)
 {
   vec = _mm256_i32gather_ps (&data[i*y], vindex, 4);
         _mm256_store_ps (buffer[i], vec);
 }

/* 2) Perform functions */
fft(buffer, x) ;

/*3) write back buffer into matrix*/
/* strided write??? ...*/

我想找到一种非常有效的方法来执行以下操作：


从矩阵中提取列：col1 = 0，4096，8192，... col2 = 1，4097，8193，...
我尝试了collect_ps，这真的很慢。
在提取的列上执行高效算法...
将列存储回矩阵：


有什么特别的把戏吗？
如何使用Intel指令集读取和写入大步距（例如4096）？

还是有任何内存操作选项可以将列从矩阵中取出？

谢谢！

Answer 1

[对于行主要数据，SIMD对行的访问速度快，而对列的访问速度慢]


是的，这就是x86-64和类似体系结构的本质。访问内存中的连续数据很快，但是访问分散的数据（无论是随机的还是规则的模式）的速度都很慢。这是拥有处理器缓存的结果。

有两种基本方法：将数据复制到有助于更好访问模式的新顺序，或者以允许更好访问模式的顺序进行计算。

不，没有经验法则或千篇一律的技巧可以使一切正常。实际上，即使比较不同的实现也很棘手，因为存在许多复杂的交互（从缓存等待时间到操作交错，再到缓存和内存访问模式），因此结果在很大程度上取决于特定的硬件和手头的数据集。



让我们看一下典型的示例情况，矩阵矩阵乘法。假设我们使用标准的C行主要数据顺序将两个5×5矩阵相乘（c = a×b）：

c00 c01 c02 c03 c04   a00 a01 a02 a03 a04   b00 b01 b02 b03 b04
c05 c06 c07 c08 c09   a05 a06 a07 a08 a09   b05 b06 b07 b08 b09
c10 c11 c12 c13 c14 = a10 a11 a12 a13 a14 × b10 b11 b12 b13 b14
c15 c16 c17 c18 c19   a15 a16 a17 a18 a19   b15 b16 b17 b18 b19
c20 c21 c22 c23 c24   a20 a21 a22 a23 a24   b20 b21 b22 b23 b24

如果将结果写为具有五个分量的垂直SIMD向量寄存器，则

c00   a00   b00   a01   b05   a02   b10   a03   b15   a04   b20
c01   a00   b01   a01   b06   a02   b11   a03   b16   a04   b21
c02 = a00 × b02 + a01 × b07 + a02 × b12 + a03 × b17 + a04 × b22
c03   a00   b03   a01   b08   a02   b13   a03   b18   a04   b23
c04   a00   b04   a01   b09   a02   b14   a03   b19   a04   b24

c05   a05   b00   a06   b05   a07   b10   a08   b15   a09   b20
c06   a05   b01   a06   b06   a07   b11   a08   b16   a09   b21
c07 = a05 × b02 + a06 × b07 + a07 × b12 + a08 × b17 + a09 × b22
c08   a05   b03   a06   b08   a07   b13   a08   b18   a09   b23
c09   a05   b04   a06   b09   a07   b14   a08   b19   a09   b24

等等。换句话说，如果 c与 b具有相同的顺序，则我们可以将具有连续存储内容的SIMD寄存器用于 c和 b，并且仅收集 a。此外， a的SIMD寄存器的所有组件都具有相同的值。

但是请注意， b寄存器对 c的所有五行重复。因此，最好将 c初始化为零，然后对具有相同 b SIMD寄存器的乘积进行加法运算：

c00    a00   b00   c05    a05   b00   c10    a10   b00   c15    a15   b00   c20    a20   b00
c01    a00   b01   c06    a05   b01   c11    a10   b01   c16    a15   b01   c21    a20   b01
c02 += a00 × b02,  c07 += a05 × b02,  c12 += a10 × b02,  c17 += a15 × b02,  c22 += a20 × b02
c03    a00 × b03   c08    a05   b03   c13    a10   b03   c18    a15   b03   c23    a20   b03
c04    a00 × b04   c09    a05   b04   c14    a10   b04   c19    a15   b04   c24    a20   b04

如果我们先转置 a，则 a的SIMD矢量寄存器也将从连续的存储器位置中获取值。实际上，如果 a足够大，则线性化 a的内存访问模式也可以提供足够大的速度提升，从而可以更快地进行转置副本（对浮点使用 uint32_t，对双精度使用 uint64_t；即完全不使用SIMD或浮点进行转置，而只是按转置顺序复制存储）。

注意，具有列主要数据顺序的情况，即与上面相比转置的数据顺序的情况非常相似。这里有很深的对称性。例如，如果 c和 b具有相同的数据顺序，而 a具有相反的数据顺序，则可以有效地对矩阵乘积进行SIMD矢量化，而不必复制任何数据。仅求和是不同的，因为这取决于数据顺序，并且矩阵乘法不是可交换的（a×b！= b×a）。

显然，主要的缺点是SIMD向量寄存器的大小是固定的，因此，不能像上面的示例那样使用完整的行作为寄存器，而只能使用部分行。 （如果结果中的列数不是SIMD寄存器宽度的倍数，则也要担心该部分向量。）

SSE和AVX具有相对大量的寄存器（8、16或32，取决于所使用的扩展集），并且取决于特定的处理器类型，可能能够同时执行至少一些矢量操作，或者至少执行较少的操作如果不相关的向量运算被交错，则延迟。因此，甚至可以选择一次操作块的宽度，以及该块是扩展向量还是块子矩阵，都需要讨论，测试和比较。


那么如何使用SIMD最有效地进行矩阵矩阵乘法呢？


就像我说的那样，这取决于数据集。恐怕没有简单的答案。

（选择最有效的方法）主要参数是被乘数和结果矩阵的大小和存储顺序。

如果您计算两个以上不同大小的矩阵的乘积，它将变得更加有趣。这是因为操作次数取决于产品的顺序。




你为什么这么气our？


我不是，实际上。以上所有这些意味着没有太多的人可以处理这种复杂性并保持理智和高效，因此有许多未发现的方法，并且可以在现实世界中获得很多收益。

即使我们忽略SIMD内部函数编译器提供的（在这种情况下为 <x86intrin.h>），我们在设计内部数据结构时也可以应用上面的逻辑，从而使我们使用的C编译器有最好的机会为我们向量化计算。 （不过，它们还不是很擅长。就像我说的那样，复杂的东西。有些人喜欢Fortran比C强，因为它的表达式和规则使Fortran编译器更容易优化和矢量化它们。）

如果这是简单还是容易的话，那么解决方案将是众所周知的。但事实并非如此，因为事实并非如此。但这并不意味着这是不可能的或我们无法实现的。它的全部意思是，足够聪明的开发人员尚未投入足够的精力来解决这个问题。