-1

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

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()。我是否以错误的方式索引到矩阵，以至于它无法从正确的地址加载？还是对齐内存的问题？

4

1 回答 1

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），但语法看起来与 for 相同float A[100][100]，尽管 asm 中的含义非常不同。无论如何，在 Fortran 二维数组中，索引是另一种方式，增加最左边的索引会将您带到内存中的下一个位置。但是在 C 中，将左索引改变一个会将您带到一个全新的行。（由 的不同元素指向float **data，或在适当的 2D 数组中，相距一排。）当然，由于这种混合与 using 相结合，您跨出了 8 行x*8。

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

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

您还需要放弃指针数组并进行手动 2D 索引（data[row * row_stride + col]因此您的数据都在一个连续的块中，而不是单独分配每一行。首先进行此更改，然后再花任何时间手动矢量化，似乎是最有意义的。

gcc 或 clang with-O3在自动矢量化标量 C 方面应该做得不错，特别是如果您使用-ffast-math. （您可以在完成手动矢量化后删除-ffast-math，但在使用自动矢量化进行调整时使用它）。

有关的：

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

于 2017-11-21T03:36:26.900 回答