如何帮助gcc矢量化C代码
我有以下C代码。 第一部分只是从标准中读入一个复数的矩阵，称为M矩阵。 有趣的部分是第二部分。
#include #include #include #include #include int main() { int n, m, c, d; float re, im; scanf("%d %d", &n, &m); assert(n==m); complex float M[n][n]; for(c=0; c<n; c++) { for(d=0; d<n; d++) { scanf("%f%fi", &re, &im); M[c][d] = re + im * I; } } for(c=0; c<n; c++) { for(d=0; d<n; d++) { printf("%.2f%+.2fi ", creal(M[c][d]), cimag(M[c][d])); } printf("\n"); } /* Example:input 2 3 1+2i 2+3i 74-4i 3+4i 4+5i -7-8i */ /* Part 2. M is now an n by n matrix of complex numbers */ int s=1, i, j; int *f = malloc(n * sizeof *f); complex float *delta = malloc(n * sizeof *delta); complex float *v = malloc(n * sizeof *v); complex float p = 1, prod; for (i = 0; i < n; i++) { v[i] = 0; for (j = 0; j <n; j++) { v[i] += M[j][i]; } p *= v[i]; f[i] = i; delta[i] = 1; } j = 0; while (j < n-1) { prod = 1.; for (i = 0; i < n; i++) { v[i] -= 2.*delta[j]*M[j][i]; prod *= v[i]; } delta[j] = -delta[j]; s = -s; p += s*prod; f[0] = 0; f[j] = f[j+1]; f[j+1] = j+1; j = f[0]; } free(delta); free(f); free(v); printf("%f + i%f\n", creal(p/pow(2.,(n-1))), cimag(p/pow(2.,(n-1)))); return 0; }
我使用gcc -fopt-info-vec-all -O3 -ffast-math -march=bdver2 permanent-in-cc -lm 。 这向我解释了为什么几乎没有循环被矢量化。
性能最重要的部分是47-50行，它们是：
for (i = 0; i < n; i++) { v[i] -= 2.*delta[j]*M[j][i]; prod *= v[i]; }
gcc告诉我：
permanent-in-cc:47:7: note: reduction used in loop. permanent-in-cc:47:7: note: Unknown def-use cycle pattern. permanent-in-cc:47:7: note: reduction used in loop. permanent-in-cc:47:7: note: Unknown def-use cycle pattern. permanent-in-cc:47:7: note: Unsupported pattern. permanent-in-cc:47:7: note: not vectorized: unsupported use in stmt. permanent-in-cc:47:7: note: unexpected pattern. [...] permanent-in-cc:48:26: note: SLP: step doesn't divide the vector-size. permanent-in-cc:48:26: note: Unknown alignment for access: IMAGPART_EXPR permanent-in-cc:48:26: note: SLP: step doesn't divide the vector-size. permanent-in-cc:48:26: note: Unknown alignment for access: REALPART_EXPR [...] permanent-in-cc:48:26: note: Build SLP failed: unrolling required in basic block SLP permanent-in-cc:48:26: note: Failed to SLP the basic block. permanent-in-cc:48:26: note: not vectorized: failed to find SLP opportunities in basic block.
如何解决阻止此部分被矢量化的问题？
奇怪的是这部分是矢量化的，但我不确定为什么：
for (j = 0; j <n; j++) { v[i] += M[j][i];
gcc -fopt-info-vec-all -O3 -ffast-math -march = bdver2 permanent-in-cc -lm的完整输出位于https://bpaste.net/show/18ebc3d66a53 。
我想我可能已经弄明白了。 经过大量的试验/错误后，很明显gcc内置的矢量化优化是一种硬编码，它不能正确理解复数。 我在代码中做了一些更改，并通过gcc输出确认了内部性能敏感循环的向量化（虽然我不确定所需的结果是否与您想要的计算等效）。 虽然我的理解仅限于您希望代码执行的操作，但我们的结论是，如果您分别计算真实和图像，它将会正常工作。 看一看：
float t_r = 0.0, t_im = 0.0; // two new temporaries while (j < n-1) { prod = 1.; for (i = 0; i < n; i++) { // fill the temps after subtraction from V to avoid stmt error t_r = creal (v[i]) - (2. * creal(delta[j]) * creal (M[j][i])); t_im = cimag(v[i]) - (2. * cimag(delta[j]) * cimag (M[j][i])) * I; //v[i] = 2.*delta[j]*M[j][i]; v[i] = t_r + t_im; // sum of real and img prod *= v[i]; } delta[j] = -delta[j]; s = -s; p += s*prod; f[0] = 0; f[j] = f[j+1]; f[j+1] = j+1; j = f[0]; }
我们先来详细检查一下代码。 我们有
complex float M[rows][cols]; complex float v[cols]; float delta[rows]; complex float p = 1.0f; float s = 1.0f;
虽然delta[]在OP代码中是complex float类型，但它只包含-1.0f或+1.0f 。 （此外，如果计算结果为-2.0f或+2.0f ，则可以简化计算。）因此，我定义的是真实的而不是复杂的。
类似地，OP将s定义为int ，但仅将其有效地用作-1.0f和+1.0f （在计算中）。 这就是我将它明确定义为float 。
我省略了f数组，因为有一种简单的方法可以完全避免它。
代码中有趣部分的第一个循环，
for (i = 0; i < n; i++) { v[i] = 0; for (j = 0; j
执行几件事。 它将delta[]数组中的所有元素初始化为1; 它可以（也可能应该）分成一个单独的循环。
由于外环在i增加， p将是v元素的乘积; 它也可以拆分成一个单独的循环。
因为内循环将列i所有元素与v[i]相加，所以外循环和内循环简单地将每一行（作为向量）与向量v相加。
因此，我们可以用伪代码重写上面的内容
Copy first row of matrix M to vector v For r = 1 .. rows-1: Add complex values in row r of matrix M to vector v p = product of complex elements in vector v delta = 1.0f, 1.0f, 1.0f, .., 1.0f, 1.0f
让我们接下来看看第二个嵌套循环：
j = 0; while (j < n-1) { prod = 1.; for (i = 0; i < n; i++) { v[i] -= 2.*delta[j]*M[j][i]; prod *= v[i]; } delta[j] = -delta[j]; s = -s; p += s*prod; f[0] = 0; f[j] = f[j+1]; f[j+1] = j+1; j = f[0]; }
除非你在循环过程中检查j的值，否则很难看到，但是外循环体中的最后4行实现了j （0,1,0,2,0,1中的OEIS A007814整数序列， 0,3,0,1,0,2,0,1,0,4，......）。 此循环中的迭代次数为2 行 -1 -1。 序列的这一部分是对称的，并实现了高度为行-1的二叉树：
4 3 3 2 2 2 2 (Read horizontally) 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
事实certificate，如果我们遍历i = 1 ... 2 行-1 ，则r是i的零低位数。 GCC提供了一个__builtin_ctz()内置函数，它可以完全计算出来。 （注意__builtin_ctz(0)产生一个未定义的值;所以不要这样做，即使它碰巧在你的计算机上产生一个特定的值。）
内环减去矩阵的行j上的复数值，从矢量v[]缩放2*delta[j] v[] 。 它还计算向量v[] （减法后）中的复数条目到变量prod的乘积。
在内循环之后， delta[j]被否定，比例因子s 。 由s缩放的变量prod的值被添加到p 。
在循环之后，最终（复杂）结果是p除以2 行-1 。 使用ldexp() C99函数（分别在实部和复杂部分上）可以做得更好。
因此，我们可以用伪代码重写第二个嵌套循环
s = 1.0f For k = 1 .. rows-1, inclusive: r = __builtin_ctz(k), ie number of least significant bits that are zero in k Subtract the complex values on row r of matrix M, scaled by delta[r], from vector v[] prod = the product of (complex) elements in vector v[] Negate scale factor s (changing its sign) Add prod times s to result p
根据我的经验，最好将实部和虚部分成单独的向量和矩阵。 考虑以下定义：
typedef struct { float *real; float *imag; size_t floats_per_row; /* Often 'stride' */ size_t rows; size_t cols; } complex_float_matrix; /* Set an array of floats to a specific value */ void float_set(float *, float, size_t); /* Copy an array of floats */ void float_copy(float *, const float *, size_t); /* malloc() vector-aligned memory; size in floats */ float *float_malloc(size_t); /* Elementwise addition of floats */ void float_add(float *, const float *, size_t); /* Elementwise addition of floats, scaled by a real scale factor */ void float_add_scaled(float *, const float *, float, size_t); /* Complex float product, separate real and imag arrays */ complex float complex_float_product(const float *, const float *, size_t);
只要float_malloc()产生充分对齐的指针（并且编译器通过例如GCC函数属性__attribute__ ((__assume_aligned__ (BYTES_IN_VECTOR))); ），并且矩阵中的floats_per_row成员被告知，所有上述内容都很容易向量化。是向量中浮点数的倍数。
（我不知道GCC是否可以自动矢量化上述函数，但我知道它们可以使用GCC向量扩展“手动”进行矢量化。）
有了上述内容，伪C中代码的整个有趣部分就变成了
complex float permanent(const complex_float_matrix *m) { float *v_real, *v_imag; float *scale; /* OP used 'delta' */ complex float result; /* OP used 'p' */ complex float product; /* OP used 'prod' */ float coeff = 1.0f; /* OP used 's' */ size_t i = 1 << (m->rows - 1); size_t r; if (!m || m->cols < 1 || m->rows < 1 || !i) { /* TODO: No input matrix, or too many rows in input matrix */ } v_real = float_malloc(m->cols); v_imag = float_malloc(m->cols); scale = float_malloc(m->rows - 1); if (!v_real || !v_imag || !scale) { free(scale); free(v_imag); free(v_real); /* TODO: Out of memory error */ } float_set(scale, 2.0f, m->rows - 1); /* Sum matrix rows to v. */ float_copy(v_real, m->real, m->cols); for (r = 1; r < m->rows; r++) float_add(v_real, m->real + r * m->floats_per_row, m->cols); float_copy(v_imag, m->imag, m->cols); for (r = 1; r < m->rows; r++) float_add(v_imag, m->imag + r * m->floats_per_row, m->cols); result = complex_float_product(v_real, v_imag, m->cols); while (--i) { r = __builtin_ctz(i); scale[r] = -scale[r]; float_add_scaled(v_real, m->real + r * m->floats_per_row, m->cols); float_add_scaled(v_imag, m->imag + r * m->floats_per_row, m->cols); product = complex_float_product(v_real, v_imag, m->cols); coeff = -coeff; result += coeff * product; } free(scale); free(v_imag); free(v_real); return result; }
在这一点上，我将亲自实现上面没有矢量化，然后广泛测试，直到我确定它正常工作。
然后，我将检查GCC程序集输出（ -S ）以查看它是否可以充分地向量化各个操作（我之前列出的函数）。
使用GCC向量扩展手动矢量化函数非常简单。 声明浮动向量很简单：
typedef float vec2f __attribute__((vector_size (8), aligned (8))); /* 64 bits; MMX, 3DNow! */ typedef float vec4f __attribute__((vector_size (16), aligned (16))); /* 128 bits; SSE */ typedef float vec8f __attribute__((vector_size (32), aligned (32))); /* 256 bits; AVX, AVX2 */ typedef float vec16f __attribute__((vector_size (64), aligned (64))); /* 512 bits; AVX512F */
每个向量中的各个组件可以使用数组表示法进行寻址（ v[0]和v[1]用于vec2f v; ）。 GCC可以对整个向量进行基本操作; 我们这里只需要加法和乘法。 应避免水平操作（在同一向量中的元素之间应用的操作），而是重新排序元素。
即使在没有这种矢量化的架构上，GCC也会为上述矢量大小生成工作代码，但结果代码可能很慢。 （至少5.4的GCC版本会产生很多不必要的动作，通常是堆叠和返回。）
为矢量分配的存储器需要充分对齐。 malloc()在所有情况下都没有提供足够对齐的内存; 你应该使用posix_memalign()代替。 当在本地或静态分配一个GCC时， aligned属性可用于增加GCC用于向量类型的对齐。 在矩阵中，您需要确保行以足够对齐的边界开始; 这就是我在结构中有floats_per_row变量的原因。
如果向量（或行）中的元素数量很大，但不是向量中浮点数的倍数，则应使用“惰性”值填充向量 - 不影响结果的值，像加法和减法的1.0f ，乘法的1.0f 。
至少在x86和x86-64上，GCC将仅使用指针为循环生成更好的代码。 例如，这个
void float_set(float *array, float value, size_t count) { float *const limit = array + count; while (array < limit) *(array++) = value; }
产生更好的代码
void float_set(float *array, float value, size_t count) { size_t i; for (i = 0; i < count; i++) array[i] = value; }
要么
void float_set(float *array, float value, size_t count) { while (count--) *(array++) = value; }
（往往产生类似的代码）。 就个人而言，我将其实施为
void float_set(float *array, float value, size_t count) { if (!((uintptr_t)array & 7) && !(count & 1)) { uint64_t *const end = (uint64_t *)__builtin_assume_aligned((void *)(array + count), 8); uint64_t *ptr = (uint64_t *)__builtin_assume_aligned((void *)array, 8); uint64_t val; __builtin_memcpy(&val, &value, 4); __builtin_memcpy(4 + (char *)&val, &value, 4); while (ptr < end) *(ptr++) = val; } else { uint32_t *const end = (uint32_t *)__builtin_assume_aligned((void *)(array + count), 4); uint32_t *ptr = (uint32_t *)__builtin_assume_aligned((void *)array, 4); uint32_t val; __builtin_memcpy(&val, &value, 4); while (ptr < end) *(ptr++) = val; } }
和float_copy()一样
void float_copy(float *target, const float *source, size_t count) { if (!((uintptr_t)source & 7) && !((uintptr_t)target & 7) && !(count & 1)) { uint64_t *const end = (uint64_t *)__builtin_assume_aligned((void *)(array + count), 8); uint64_t *ptr = (uint64_t *)__builtin_assume_aligned((void *)target, 8); uint64_t *src = (uint64_t *)__builtin_assume_aligned((void *)source, 8); while (ptr < end) *(ptr++) = *(src++); } else { uint32_t *const end = (uint32_t *)__builtin_assume_aligned((void *)(array + count), 4); uint32_t *ptr = (uint32_t *)__builtin_assume_aligned((void *)array, 4); uint32_t *src = (uint32_t *)__builtin_assume_aligned((void *)source, 4); while (ptr < end) *(ptr++) = *(src++); } }
或接近那个。
矢量化最难的是complex_float_product() 。 如果在最终向量中填入未使用的元素，实部为1.0f ，虚部为0.0f ，则可以轻松计算每个向量的复杂乘积。 记住这一点
（ a + b i）×（ c + d i）=（ a c - b d ）+（ a d + b c ）i
这里的难点在于有效地计算向量中元素的复杂乘积。 幸运的是，该部分对整体性能并不重要（除了非常短的向量或列数很少的矩阵），因此在实践中它应该不重要。
（简而言之，“硬”部分是找到重新排序元素的方法，以最大限度地使用打包向量乘法，而不需要这么多的shuffle / move来减慢它。）
对于float_add_scaled()函数，您应该创建一个填充了比例因子的向量; 如下所示，
void float_add_scaled(float *array, const float *source, float scale, size_t count) { const vec4f coeff = { scale, scale, scale, scale }; vec4f *ptr = (vec4f *)__builtin_assume_aligned((void *)array, 16); vec4f *const end = (vec4f *)__builtin_assume_aligned((void *)(array + count), 16); const vec4f *src = (vec4f *)__builtin_assume_aligned((void *)source, 16); while (ptr < end) *(ptr++) += *(src++) * coeff; }
如果我们忽略对齐和大小检查，以及回退实现。
优化程序日志清楚地表明
访问的未知对齐方式：…
在尝试矢量化时
printf("%.2f%+.2fi ", creal(M[c][d]), cimag(M[c][d])); //24 v[i] += M[j][i]; //38 p *= v[i]; //40 v[i] -= 2.*delta[j]*M[j][i]; //48
看起来真的很麻烦你需要在内存中强制对齐数组M ， delta和v 。
GCC中的自动矢量化
仅处理对齐的内存访问（不要尝试向量化包含未对齐访问的循环）
正如之前的评论中提到的，我建议你为此目的使用posix_memalign 。
complex float * restrict delta; posix_memalign(&delta, 64, n * sizeof *delta); //to adapt
你的目标环境是什么？ （OS，CPU）
请查看数据对齐到辅助向量化