我目前正在尝试使用内在函数和循环展开来优化矩阵运算。存在我无法弄清楚的段错误。这是我更改的代码:
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()
。我是否以错误的方式索引到矩阵,以致无法从正确的地址加载?还是内存对齐的问题?
最佳答案
至少在一个地方,您正在执行 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
,但在自动矢量化调整时使用它)。
相关:
另请参阅我对 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 .
关于c - 矩阵优化 - 使用内在函数和循环展开时出现段错误,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/47402422/