C++ 在 C 中使用 AVX 实现矩阵运算

发布于05月01日

我正在try 使用AVX实现以下操作:

for (i=0; i<N; i++) {
  for(j=0; j<N; j++) {
    for (k=0; k<K; k++) {
      d[i][j] += 2 * a[i][k] * ( b[k][j]- c[k]);
    }
  }
}

for (int i=0; i<N; i++){
   f+= d[ind[i]][ind[i]]/2;
}

其中d是NxN矩阵，a是NxK，b是KxN，c是长度K的向量.它们都是双精度的.当然，所有的数据都是统一的，我使用#pragma vector aligned帮助编译器(gcc).

我知道如何在一维数组中使用AVX扩展，但使用矩阵对我来说有点棘手.目前，我有以下几点，但我没有得到正确的结果:

    for (int i=0; i< floor (N/4); i++){
        for (int j=0; j< floor (N/4); j++){
            __m256d D, A, B, C;
            D = _mm256_setzero_pd();
            #pragma vector aligned
            for (int k=0; k<K_MAX; k++){
                A = _mm256_load_pd(a[i] + k*4);
                B = _mm256_load_pd(b[k] + j*4);
                C = _mm256_load_pd(c + 4*k);
                B = _mm256_sub_pd(B, C);
                A = _mm256_mul_pd(A, B);
                D = _mm256_add_pd(_mm256_set1_pd(2.0), A);
                _mm256_store_pd(d[i] + j*4, D);
            }

        }
    }


    for (int i=0; i<N; i++){
        f+= d[ind[i]][ind[i]]/2;
    }

我希望有人能告诉我哪里出了错.

提前谢谢.

注意:我不想介绍OpenMP，只是使用SIMD英特尔指令

#include <stdint.h> #include <stdlib.h> #include <immintrin.h> // This double [N][N] C99 VLA syntax isn't portable to C++ even with GNU extensions // restrict tells the compiler the output doesn't overlap with any of the inputs void matop(size_t N, size_t K, double d[restrict N][N], const double a[restrict N][K], const double b[restrict K][N], const double c[restrict K]) { for( size_t i = 0; i < N; i++ ) { // loop-invariant pointers for this outer iteration //double* restrict rowDi = &d[ i ][ 0 ]; const double* restrict rowAi = &a[ i ][ 0 ]; for( size_t k = 0; k < K; k++ ) { const double* restrict rowBk = &b[ k ][ 0 ]; double* restrict rowDi = &d[ i ][ 0 ]; #if 0 // pure scalar // auto-vectorizes ok; still a lot of extra checking outside outermost loop even with restrict for (size_t j=0 ; j<N ; j++){ rowDi[j] += 2*rowAi[k] * (rowBk[j] - c[k]); } #else // SIMD inner loop with cleanup // *** TODO: unroll over 2 or 3 i values // and maybe also 2 or 3 k values, to reuse each bv a few times while it's loaded. __m256d av = _mm256_broadcast_sd( rowAi + k ); av = _mm256_add_pd( av, av ); // 2*a[ i ][ k ] broadcasted const __m256d cv = _mm256_broadcast_sd( &c[ k ] ); const __m256d minus_ck = _mm256_xor_pd(cv, _mm256_set1_pd(-0.0)); // broadcasted -c[k] //const size_t N_aligned = ( (size_t)N / 4 ) * 4; size_t N_aligned = N & -4; // round down to a multiple of 4 j iterations const double* endBk = rowBk + N_aligned; //for( ; j < N_aligned; j += 4 ) for ( ; rowBk != endBk ; rowBk += 4, rowDi += 4) { // coax GCC into using pointer-increments in the asm, instead of j+=4 // Load the output vector to update __m256d dv = _mm256_loadu_pd( rowDi ); // Update with FMA __m256d bv = _mm256_loadu_pd( rowBk ); __m256d t2 = _mm256_add_pd( minus_ck, bv ); // bv - cv dv = _mm256_fmadd_pd( av, t2, dv ); // Store back to the same address _mm256_storeu_pd( rowDi, dv ); } // rowDi and rowBk point to the double after the last full vector // The remainder, if you can't pad your rows to a multiple of 4 and step on that padding for(int j=0 ; j < (N&3); j++ ) rowDi[ j ] += _mm256_cvtsd_f64( av ) * ( rowBk[ j ] + _mm256_cvtsd_f64( minus_ck ) ); #endif } } }

.L7: # do{ vaddpd ymm0, ymm2, YMMWORD PTR [rax] # -c[k] + rowBk[0..3] add rax, 32 # rowBk += 4 add rdx, 32 # rowDi += 4 vfmadd213pd ymm0, ymm1, YMMWORD PTR [rdx-32] # fma(2aik, Bkj-ck, Dij) vmovupd YMMWORD PTR [rdx-32], ymm0 # store FMA result cmp rcx, rax jne .L7 # }while(p != endp)

size_t j = 0; // used for cleanup loop after for( ; j < N_aligned; j += 4 ) { // Load the output vector to update __m256d dv = _mm256_loadu_pd( rowDi + j ); // Update with FMA __m256d bv = _mm256_loadu_pd( rowBk + j ); __m256d t2 = _mm256_sub_pd( bv, cv ); // bv - cv dv = _mm256_fmadd_pd( av, t2, dv ); // Store back to the same address _mm256_storeu_pd( rowDi + j, dv ); } // The remainder, if you can't pad your rows to a multiple of 4 for( ; j < N; j++ ) rowDi[ j ] += _mm256_cvtsd_f64( av ) * ( rowBk[ j ] - _mm256_cvtsd_f64( cv ) );

智能展开:缓存阻塞的第一步

如果我们在加载每个数据向量时对其进行更多处理，那么上述一些问题可能会消失，这些问题包括最大化加载/存储执行单元吞吐量，以及要求编译器使用非索引寻址模式.

例如，我们可以处理2、3或4，而不是只处理一行d[][].然后，对于d[i+unroll][j + 0..vec]向量，每个(rowBk[j] - c[k])结果可以使用多次(使用不同的2aik).

我们还可以加载几个不同的(rowBk+K*0..unroll)[j+0..3]，每个都有相应的minus_ck0、minus_ck1等等(或者保留一个向量数组；只要它很小，编译器有足够的寄存器，这些元素就不会存在于内存中)

由于寄存器中同时有多个bv-cv和dv个向量，我们可以在不增加FP工作总量的情况下，按负载执行更多的FMA.然而，常数需要更多的寄存器，否则我们可能会通过强制更多的重新加载来达到目的.

d[i][j] += (2*a_ik) * b[k][j] - (2*a_ik*c_k)转换在这里没有用；我们希望将bv-cv和i分开，这样我们就可以将结果作为不同FMA的输入进行重用.

b[k][j]+(-c[k])仍然可以从负载与vaddpd的微融合中获益，因此理想情况下，它仍将使用指针增量，但前端可能不再是瓶颈.

不要做得过火；太多的内存输入流可能是缓存冲突未命中的问题，尤其是对于可能会产生混叠的某些N值，以及对跟踪它们的硬件预取的问题.(虽然英特尔的L2拖缆据说每4k页跟踪1个正向流和1个反向流，但IIRC.)大概4到8次就可以了.但如果L1d中没有d[][]，那么它就不是真正的内存输入流.不过，您不希望b[][]个输入行逐出d个数据，因为您将重复循环2到4行d个数据.

相比之下:Soons的循环——清理频率较低，但内存访问模式更差.

Soonts目前的3个负载和3个ALU操作的循环并不理想，尽管每个FMA操作1个负载在缓存中已经可以了(大多数现代CPU每个时钟可以执行2个，尽管AMD Zen也可以与mul/FMA并行执行2个FP加法).如果额外的ALU操作是一个瓶颈，我们可以将a[][]乘以2一次，只需要O(N*K)个工作，而不是O(N^2*K)个工作就可以完成.但这可能不是瓶颈，因此不值得.

更重要的是，Soons当前答案中的内存访问模式是，对于c[k]和a[i][k]的广播负载，一次向前循环1倍，这很好，但是the 102 of 103 is unfortunately striding down a column.

如果你要像Soons建议的那样展开，不要只为一个dv做两个dep链，至少要做两个向量，d[i][j + 0..3]和4..7，这样你就可以从你touch 的每b[k][j]处使用整个64字节(完整缓存线).或一对缓存线的四个向量.(英特尔CPU至少使用了一个邻线预取器，它喜欢完成128字节aligned对缓存线，因此您可以将b[][]行与128行或至少与64行对齐，并从邻线预取中获得一些好处.

如果b[][]的垂直切片适合某个级别的缓存(以及当前累积到的d[i][]行)，那么下一步下一组列可以受益于预取和局部性.如果没有，那么充分利用你接触到的线就更重要了，这样以后就不必再拉它们了.

因此，对于Soons的矢量化策略，对于L1d缓存中不适用的大型问题，最好确保b的行按64对齐，即使这意味着在每行末尾填充.(存储几何体不必与实际的矩阵维数匹配；分别传递N和row_stride.一个用于索引计算，另一个用于循环边界.)

C++ 在 C 中使用 AVX 实现矩阵运算

推荐答案

智能展开:缓存阻塞的第一步

相比之下:Soons的循环——清理频率较低，但内存访问模式更差.

C++相关问答推荐

初始化char数组-用于初始化数组的字符串是否除了存储数组的位置之外单独存储在内存中

理解没有返回语句的递归C函数的行为

如何正确地索引C中的 struct 指针数组？

为什么在C中进行大量的位移位？

堆栈帧和值指针

如何使fputs功能提示错误输入并要求用户重新输入.程序停止而不是请求新的输入

每个 struct 变量在C中都有自己的命名空间吗？

Char变量如何在不使用方括号或花括号的情况下存储字符串，以及它如何迭代到下一个字符？

CSAPP微型shell 实验室：卡在sigprocmask

==284==错误：AddressSaniizer：堆栈缓冲区下溢

在C语言中，指针指向一个数组

S和查尔有什么不同[1]？

变量值不正确的问题

哪些C++功能可以在外部C块中使用

关于不同C编译器中的attribute支持

clion.我无法理解 Clion 中发生的 scanf 错误

C 中从 Unix 纪元时间转换的损坏

为什么需要struct in_addr

仅使用其内存地址取消引用 C 中的 struct

为什么这里的符号没有解析？