我在C中使用AVX2、FMA制作了矩阵向量乘法程序。我使用GCC ver7和-mfma、-mavx编译。
然而，我得到了错误"释放对象的校验和不正确-对象可能在被释放后被修改了。">
我认为如果矩阵的维数不是4的倍数，就会产生误差。
我知道AVX2使用ymm寄存器，可以使用4个双精度浮点数字。因此，在矩阵是4的倍数的情况下，我可以无误差地使用AVX2。
但是，这是我的问题。如果矩阵不是4的倍数，我如何有效地使用AVX2？？？
这是我的密码。
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#include "time.h"
#include "x86intrin.h"
void mv(double *a,double *b,double *c, int m, int n, int l)
{
__m256d va,vb,vc;
int k;
int i;
for (k = 0; k < l; k++) {
vb = _mm256_broadcast_sd(&b[k]);
for (i = 0; i < m; i+=4) {
va = _mm256_loadu_pd(&a[m*k+i]);
vc = _mm256_loadu_pd(&c[i]);
vc = _mm256_fmadd_pd(vc, va, vb);
_mm256_storeu_pd( &c[i], vc );
}
}
}
int main(int argc, char* argv[]) {
// set variables
int m;
double* a;
double* b;
double* c;
int i;
int temp=0;
struct timespec startTime, endTime;
m=9;
// main program
// set vector or matrix
a=(double *)malloc(sizeof(double) * m*m);
b=(double *)malloc(sizeof(double) * m*1);
c=(double *)malloc(sizeof(double) * m*1);
for (i=0;i<m;i++) {
a[i]=1;
b[i]=1;
c[i]=0.0;
}
for (i=m;i<m*m;i++) {
a[i]=1;
}
// check start time
clock_gettime(CLOCK_REALTIME, &startTime);
mv(a, b, c, m, 1, m);
// check end time
clock_gettime(CLOCK_REALTIME, &endTime);
free(a);
free(b);
free(c);
return 0;
}
您加载并存储4个double的向量，但循环条件仅检查第一个向量元素是否在边界内，因此当m不是4的倍数时，您最多可以写入3x8=24字节的外部对象。
在主循环中需要类似i < (m-3)的东西，以及处理数据的最后一个部分向量的清理策略。SIMD的矢量化非常类似于展开：您必须检查在循环条件下执行多个未来元素是否可以。
标量清理循环运行良好，但我们可以做得更好。例如，在最后一个完整的256位矢量(即，最多1个(之后，在变为标量之前，进行尽可能多的128位矢量。
在许多情况下(例如，只写目的地(，在数组末尾结束的未对齐矢量加载非常好(当m>=4时(。如果m%4 != 0，它可能会与主循环重叠，但这很好，因为输出数组不会与输入重叠，所以作为单个清理的一部分重做元素比避免它的分支更便宜。
但这在这里不起作用，因为您的逻辑是c[i+0..3] += ...，所以重做一个元素会使它出错。
// cleanup using a 128-bit FMA, then scalar if there's an odd element.
// untested
void mv(double *a,double *b,double *c, int m, int n, int l)
{
/* the loop below should actually work for m=1..3, but a separate strategy might be good.
if (m < 4) {
// maybe check m >= 2 and use __m128 vectors?
// or vectorize differently?
}
*/
for (int k = 0; k < l; k++) {
__m256 vb = _mm256_broadcast_sd(&b[k]);
int i;
for (i = 0; i < (m-3); i+=4) {
__m256d va = _mm256_loadu_pd(&a[m*k+i]);
__m256d vc = _mm256_loadu_pd(&c[i]);
vc = _mm256_fmadd_pd(vc, va, vb);
_mm256_storeu_pd( &c[i], vc );
}
if (i<(m-1)) {
__m128d lasta = _mm_loadu_pd(&a[m*k+i]);
__m128d lastc = _mm_loadu_pd(&c[i]);
lastc = _mm_fmadd_pd(lastc, va, _mm256_castpd256_pd128(vb));
_mm_storeu_pd( &c[i], lastc );
// i+=2; // last element only checks m odd/even, doesn't use i
}
// if (i<m)
if (m&1) {
// odd number of elements, do the last non-vector one
c[m-1] += a[m*k + m-1] * _mm256_cvtsd_f64(vb);
}
}
}
我还没有研究过gcc/clang-O3是如何编译的。有时编译器试图在清理代码方面变得过于聪明(例如，试图自动向量化标量清理循环(。
其他策略可能包括使用AVX掩码存储来完成最后4个元素：每个矩阵行的末尾都需要相同的掩码，因此生成一次，然后在每行的末尾使用它可能会很好。请参阅使用未对齐缓冲区矢量化：使用VMASKMOVPS：根据未对齐计数生成掩码？或者根本不使用insn。(为了简化分支，你应该设置它，使主循环只进入i < (m-4)，然后总是运行清理。在m%4 == 0的情况下，掩码是全1，所以你做最后的全向量。(如果你不能安全地读取矩阵的末尾，你可能需要一个掩码加载和掩码存储。
为了提高效率，您还可以考虑对齐行，或者与行的逻辑长度分开的行步长。(即焊盘排到32字节边界(。在行的末尾留下填充可以简化清理，因为您总是可以执行写填充的整个向量。
特殊情况m==2：您希望将2个元素广播到__m256d的两个128位通道中，而不是从b[]广播一个元素，因此一个256位FMA可以同时执行2行。