使用代码向量化的矩阵运算

Question

我已经编写了一个函数来进行 4x4 矩阵的转置，但我不知道如何扩展矩阵 m x n 的代码。

在哪里可以找到一些关于 SSE 矩阵运算的示例代码？乘积、转置、逆等？

这是转置4x4的代码：

 void transpose(float* src, int n) {
    __m128  row0,   row1,   row2,   row3;
    __m128 tmp1;
    tmp1=_mm_loadh_pi(_mm_loadl_pi(tmp1, (__m64*)(src)), (__m64*)(src+ 4));
    row1=_mm_loadh_pi(_mm_loadl_pi(row1, (__m64*)(src+8)), (__m64*)(src+12));
    row0=_mm_shuffle_ps(tmp1, row1, 0x88);
    row1=_mm_shuffle_ps(row1, tmp1, 0xDD);

    tmp1=_mm_movelh_ps(tmp1, row1);
    row1=_mm_movehl_ps(tmp1, row1);

    tmp1=_mm_loadh_pi(_mm_loadl_pi(tmp1, (__m64*)(src+ 2)), (__m64*)(src+ 6));
    row3= _mm_loadh_pi(_mm_loadl_pi(row3, (__m64*)(src+10)), (__m64*)(src+14));
    row2=_mm_shuffle_ps(tmp1, row3, 0x88);
    row3=_mm_shuffle_ps(row3, tmp1, 0xDD);

    tmp1=_mm_movelh_ps(tmp1, row3);
    row3=_mm_movehl_ps(tmp1, row3);

    _mm_store_ps(src, row0);
    _mm_store_ps(src+4, row1);
    _mm_store_ps(src+8, row2);
    _mm_store_ps(src+12, row3);
}

Answer 1

这是一种通用方法，可用于使用平铺转置 NxN 矩阵。您甚至可以使用现有的 4x4 转置并使用 4x4 平铺大小：

for each 4x4 block in the matrix with top left indices r, c
    if block is on diagonal (i.e. if r == c)
         get block a = 4x4 block at r, c
         transpose block a
         store block a at r, c
    else if block is above diagonal (i.e. if r < c)
         get block a = 4x4 block at r, c
         get block b = 4x4 block at c, r
         transpose block a
         transpose block b
         store transposed block a at c, r
         store transposed block b at r, c
    else // block is below diagonal
         do nothing
    endif
endfor

显然 N 需要是 4 的倍数才能起作用，否则您将需要做一些额外的整理工作。

正如上面在 cmets 中提到的，一个 MxN in-place 转置很难做到——你需要使用一个额外的临时矩阵（这实际上使它成为一个非就地转置）或使用here 中描述的方法，但这将难以使用 SIMD 进行矢量化。

Answer 2

我不确定如何有效地使用 SIMD 对任意矩阵进行就地转置，但我知道如何就地进行转换。让我描述一下如何做到这两点

原地转置

对于就地转置，您应该查看 Agner Fog 的 Optimizing software in C++ 手册。请参阅第 9.10 节“大型数据结构中的缓存争用”示例 9.5a。对于某些矩阵大小，由于缓存混叠，您会看到性能大幅下降。有关示例，请参见表 9.1 和此Why is transposing a matrix of 512x512 much slower than transposing a matrix of 513x513?。 Agner 在示例 9.5b 中提供了一种使用循环平铺（类似于 Paul R 所描述的）来解决此问题的方法。

不合适的转置

在这里查看我的答案（得票最多的那个）What is the fastest way to transpose a matrix in C++?。我已经很久没有研究过了，但让我在这里重复一下我的代码：

inline void transpose4x4_SSE(float *A, float *B, const int lda, const int ldb) {
    __m128 row1 = _mm_load_ps(&A[0*lda]);
    __m128 row2 = _mm_load_ps(&A[1*lda]);
    __m128 row3 = _mm_load_ps(&A[2*lda]);
    __m128 row4 = _mm_load_ps(&A[3*lda]);
     _MM_TRANSPOSE4_PS(row1, row2, row3, row4);
     _mm_store_ps(&B[0*ldb], row1);
     _mm_store_ps(&B[1*ldb], row2);
     _mm_store_ps(&B[2*ldb], row3);
     _mm_store_ps(&B[3*ldb], row4);
}

inline void transpose_block_SSE4x4(float *A, float *B, const int n, const int m, const int lda, const int ldb ,const int block_size) {
    #pragma omp parallel for
    for(int i=0; i<n; i+=block_size) {
        for(int j=0; j<m; j+=block_size) {
            int max_i2 = i+block_size < n ? i + block_size : n;
            int max_j2 = j+block_size < m ? j + block_size : m;
            for(int i2=i; i2<max_i2; i2+=4) {
                for(int j2=j; j2<max_j2; j2+=4) {
                    transpose4x4_SSE(&A[i2*lda +j2], &B[j2*ldb + i2], lda, ldb);
                }
            }
        }
    }   
}