c - 矩阵优化 - 使用内在函数和循环展开时出现段错误

Question

我目前正在尝试使用内在函数和循环展开来优化矩阵运算。存在我无法弄清楚的段错误。这是我更改的代码:

const int UNROLL = 4;

void outer_product(matrix *vec1, matrix *vec2, matrix *dst) {
    assert(vec1->dim.cols == 1 && vec2->dim.cols == 1 && vec1->dim.rows == dst->dim.rows && vec2->dim.rows == dst->dim.cols);
    __m256 tmp[4];
    for (int x = 0; x < UNROLL; x++) {
        tmp[x] = _mm256_setzero_ps();
    } 
    for (int i = 0; i < vec1->dim.rows; i+=UNROLL*8) {
        for (int j = 0; j < vec2->dim.rows; j++) {     
            __m256 row2 = _mm256_broadcast_ss(&vec2->data[j][0]);
            for (int x = 0; x<UNROLL; x++) {
                tmp[x] = _mm256_mul_ps(_mm256_load_ps(&vec1->data[i+x*8][0]), row2); 
                _mm256_store_ps(&dst->data[i+x*8][j], tmp[x]);
            } 
        }
    }
}

void matrix_multiply(matrix *mat1, matrix *mat2, matrix *dst) {
    assert (mat1->dim.cols == mat2->dim.rows && dst->dim.rows == mat1->dim.rows && dst->dim.cols == mat2->dim.cols);
    for (int i = 0; i < mat1->dim.rows; i+=UNROLL*8) {
        for (int j = 0; j < mat2->dim.cols; j++) {
        __m256 tmp[4];
            for (int x = 0; x < UNROLL; x++) {
                tmp[x] = _mm256_setzero_ps();
            } 
            for (int k = 0; k < mat1->dim.cols; k++) {
                __m256 mat2_s = _mm256_broadcast_ss(&mat2->data[k][j]);
                for (int x = 0; x < UNROLL; x++) {
                    tmp[x] = _mm256_add_ps(tmp[x], _mm256_mul_ps(_mm256_load_ps(&mat1->data[i+x*8][k]), mat2_s));
                }
            }
            for (int x = 0; x < UNROLL; x++) {
                _mm256_store_ps(&dst->data[i+x*8][j], tmp[x]);
            }    
        }
    }    
}

编辑:这是矩阵的结构。我没有修改它。

typedef struct shape {
    int rows;
    int cols;
} shape;

typedef struct matrix {
    shape dim;
    float** data;
} matrix;

编辑:我尝试 gdb 找出哪一行导致了段错误，看起来像是 _mm256_load_ps()。我是否以错误的方式索引到矩阵，以致无法从正确的地址加载？还是内存对齐的问题？

Answer 1

至少在一个地方，您正在执行 32 字节对齐所需的加载，步幅仅为 4 个字节。不过，我认为这不是您真正想要做的:

for (int k = 0; k < mat1->dim.cols; k++) {
    for (int x = 0; x < UNROLL; x++) {
        ...
        _mm256_load_ps(&mat1->data[i+x*8][k])
     }
 }

_mm256_load_ps加载 8 个连续的 float s，即加载 data[i+x*8][k]至data[i+x*8][k+7] 。 我想你想要data[i+x][k*8] ， 并循环 k在最内层循环中。

如果您需要未对齐的加载/存储，请使用 _mm256_loadu_ps/_mm256_storeu_ps 。但更喜欢将数据对齐到 32B，并填充矩阵的存储布局，以便行步长是 32 字节的倍数。 (数组的实际逻辑维度不必与步长匹配；可以将每行末尾的填充保留为 16 或 32 字节的倍数。这使得循环更容易编写。)

<小时/>

您甚至没有使用二维数组(您使用的是指向 float 数组的指针数组)，但语法看起来与 float A[100][100] 相同。 ，尽管 asm 中的含义非常不同。无论如何，在 Fortran 2D 数组中，索引采用相反的方式，递增最左边的索引会将您带到内存中的下一个位置。但在 C 语言中，将左侧索引改变 1 会将您带到一个全新的行。 (由 float **data 的不同元素指向，或者在适当的二维数组中，一行跨过。)当然，由于这种混合与使用 x*8 相结合，您跨了 8 行。 .

说到 asm，这段代码的结果非常糟糕，尤其是使用 gcc，它为每个 vector 重新加载 4 个东西，我认为是因为不确定 vector 存储不会给指针数据别名。将事物分配给局部变量以确保编译器可以将它们提升出循环。 (例如 const float *mat1dat = mat1->data; 。)Clang 做得稍好一些，但源中的访问模式本质上是不好的，并且需要对每个内循环迭代进行指针追踪才能到达新行，因为您循环了 x而不是k 。我把它放在 Godbolt compiler explorer .

<小时/>

但实际上，您应该首先优化内存布局，然后再尝试手动矢量化。转置其中一个数组可能是值得的，因此您可以在一个矩阵的行和另一个矩阵的列上循环连续内存，同时进行行和列的点积来计算结果的一个元素。或者这可能是值得做的c[Arow,Bcol] += a_value_from_A * b[Arow,Bcol]在内循环内部而不是预先转置(但这会占用大量内存流量)。但无论您做什么，请确保您不会跨过对内循环中的矩阵之一的非连续访问。

您还需要放弃指针数组并进行手动 2D 索引( data[row * row_stride + col] ，这样您的数据就全部在一个连续的 block 中，而不是单独分配每一行。在您之前先进行此更改花任何时间手动矢量化，似乎最有意义。

gcc 或 clang 与 -O3应该完成自动矢量化标量 C 的工作，特别是当您使用 -ffast-math 进行编译时。 (完成手动矢量化后，您可以删除 -ffast-math，但在自动矢量化调整时使用它)。

相关:

• How does BLAS get such extreme performance?

• 另请参阅我对 Poor maths performance in C vs Python/numpy 的评论另一个糟糕的内存布局问题。

• how to optimize matrix multiplication (matmul) code to run fast on a single processor core

您可以在查看缓存阻塞之前或之后手动进行矢量化，但是当您这样做时，请参阅 Matrix Multiplication with blocks .