我正在尝试将一个循环转换为SSE内在函数。我似乎取得了相当不错的进步，这就是说它的方向正确，但是我似乎在某处做了一些翻译错误，因为我没有从非SSE代码得到相同的“正确”答案。
我的原始循环展开了4倍，看起来像这样：
int unroll_n = (N/4)*4;
for (int j = 0; j < unroll_n; j++) {
for (int i = 0; i < unroll_n; i+=4) {
float rx = x[j] - x[i];
float ry = y[j] - y[i];
float rz = z[j] - z[i];
float r2 = rx*rx + ry*ry + rz*rz + eps;
float r2inv = 1.0f / sqrt(r2);
float r6inv = r2inv * r2inv * r2inv;
float s = m[j] * r6inv;
ax[i] += s * rx;
ay[i] += s * ry;
az[i] += s * rz;
//u
rx = x[j] - x[i+1];
ry = y[j] - y[i+1];
rz = z[j] - z[i+1];
r2 = rx*rx + ry*ry + rz*rz + eps;
r2inv = 1.0f / sqrt(r2);
r6inv = r2inv * r2inv * r2inv;
s = m[j] * r6inv;
ax[i+1] += s * rx;
ay[i+1] += s * ry;
az[i+1] += s * rz;
//unroll i 3
rx = x[j] - x[i+2];
ry = y[j] - y[i+2];
rz = z[j] - z[i+2];
r2 = rx*rx + ry*ry + rz*rz + eps;
r2inv = 1.0f / sqrt(r2);
r6inv = r2inv * r2inv * r2inv;
s = m[j] * r6inv;
ax[i+2] += s * rx;
ay[i+2] += s * ry;
az[i+2] += s * rz;
//unroll i 4
rx = x[j] - x[i+3];
ry = y[j] - y[i+3];
rz = z[j] - z[i+3];
r2 = rx*rx + ry*ry + rz*rz + eps;
r2inv = 1.0f / sqrt(r2);
r6inv = r2inv * r2inv * r2inv;
s = m[j] * r6inv;
ax[i+3] += s * rx;
ay[i+3] += s * ry;
az[i+3] += s * rz;
}
}
然后，我基本上逐行转到顶部，并将其转换为SSE内部函数。代码如下。我不确定是否需要前三行，但是我知道我的数据需要16位对齐才能正常且最佳地工作。
float *x = malloc(sizeof(float) * N);
float *y = malloc(sizeof(float) * N);
float *z = malloc(sizeof(float) * N);
for (int j = 0; j < N; j++) {
for (int i = 0; i < N; i+=4) {
__m128 xj_v = _mm_set1_ps(x[j]);
__m128 xi_v = _mm_load_ps(&x[i]);
__m128 rx_v = _mm_sub_ps(xj_v, xi_v);
__m128 yj_v = _mm_set1_ps(y[j]);
__m128 yi_v = _mm_load_ps(&y[i]);
__m128 ry_v = _mm_sub_ps(yj_v, yi_v);
__m128 zj_v = _mm_set1_ps(z[j]);
__m128 zi_v = _mm_load_ps(&z[i]);
__m128 rz_v = _mm_sub_ps(zj_v, zi_v);
__m128 r2_v = _mm_mul_ps(rx_v, rx_v) + _mm_mul_ps(ry_v, ry_v) + _mm_mul_ps(rz_v, rz_v) + _mm_set1_ps(eps);
__m128 r2inv_v = _mm_div_ps(_mm_set1_ps(1.0f),_mm_sqrt_ps(r2_v));
__m128 r6inv_1v = _mm_mul_ps(r2inv_v, r2inv_v);
__m128 r6inv_v = _mm_mul_ps(r6inv_1v, r2inv_v);
__m128 mj_v = _mm_set1_ps(m[j]);
__m128 s_v = _mm_mul_ps(mj_v, r6inv_v);
__m128 axi_v = _mm_load_ps(&ax[i]);
__m128 ayi_v = _mm_load_ps(&ay[i]);
__m128 azi_v = _mm_load_ps(&az[i]);
__m128 srx_v = _mm_mul_ps(s_v, rx_v);
__m128 sry_v = _mm_mul_ps(s_v, ry_v);
__m128 srz_v = _mm_mul_ps(s_v, rz_v);
axi_v = _mm_add_ps(axi_v, srx_v);
ayi_v = _mm_add_ps(ayi_v, srx_v);
azi_v = _mm_add_ps(azi_v, srx_v);
_mm_store_ps(ax, axi_v);
_mm_store_ps(ay, ayi_v);
_mm_store_ps(az, azi_v);
}
}
我认为主要思想是正确的，但是由于结果答案不正确，所以在某处存在一个/一些错误。
我认为您唯一的错误是简单的拼写错误，而不是逻辑错误，请参见下文。
您不能只使用clang的自动矢量化吗？还是需要为此代码使用gcc？自动向量化功能使您无需修改即可从同一来源制作SSE，AVX和（以后）AVX512版本。不幸的是，本征不能扩展到不同的向量大小。
根据您开始矢量化的过程，我制作了一个优化版本。您应该尝试一下，我很想知道它是否比带错误修正的版本或clang的自动矢量化版本更快。 ：）
_mm_store_ps(ax, axi_v);
_mm_store_ps(ay, ayi_v);
_mm_store_ps(az, azi_v);
您从ax[i]加载，但现在要存储到ax[0]。
此外，clang's unused-variable warning发现了此错误：
axi_v = _mm_add_ps(axi_v, srx_v);
ayi_v = _mm_add_ps(ayi_v, srx_v); // should be sry_v
azi_v = _mm_add_ps(azi_v, srx_v); // should be srz_v
[就像我在上一个问题的回答中说的那样，您可能应该互换循环，所以相同的ax [i + 0..3]，ay [i + 0..3]和az [i + 0。 .3]，避免加载/存储。
[此外，如果您不打算使用rsqrtps + Newton-Raphson，则应使用我在上一个答案中指出的转换：将rsqrtps除以m[j]。用sqrt(k2) ^3将1.0f除以某物，然后再仅相乘一次是没有意义的。
divps might实际上不是一个胜利，因为总的uop吞吐量可能比div / sqrt吞吐量或延迟更是一个瓶颈。 rsqrt比sqrtps + divps多得多。 three multiplies + FMA + rsqrtps对于AVX 256b向量更有用，因为在SnB到Broadwell上，divide / sqrt单位不是全角。 Skylake具有12c的延迟rsqrt，与xmm相同，但是xmm的吞吐量仍然更好（每3c 1而不是6c 1）。
sqrtps ymm。 （当然，只有使用打包版本的clang。）
如果不交换循环，则应手动将仅依赖于clang and gcc were both using rsqrtps / rsqrtss when compiling your code with -ffast-math的所有内容从内部循环中吊起。编译器通常擅长于此，但是使源代码反映您希望编译器能够完成的工作似乎仍然是一个好主意。这有助于“了解”循环的实际作用。
这里是一个对您的原作进行了一些优化的版本：
为了使gcc / clang将mul / adds融合到FMA中，我使用了rsqrtps。这无需使用rsqrtss就可获取高吞吐量的FMA指令。 （三个独立的累加器有很多并行性，因此，与-ffast-math相比，FMA的增加的延迟根本不会受到损害。我希望这里的port0 / 1吞吐量是限制因素。）我j，但似乎没有-ffp-contract=fast就不会出现在这里。
注意v 3/2 = sqrt（v）3 = sqrt（v）* v。这具有更低的延迟和更少的指令。
交换了循环，并在内部循环中使用广播负载来提高局部性（使用AVX将带宽需求减少4或8）。内部循环的每次迭代仅从四个源阵列中的每一个读取4B新数据。 （x，y，z和m）。因此，趁热时，它会充分利用每个缓存行。
在内部循环中使用广播负载还意味着我们并行累加-ffast-math，避免了需要addps。 （有关代码thought gcc did this automatically，请参见此答案的先前版本。）
它很好地编译-ffast-math。 （不过，与clang的自动矢量化256b版本不同，矢量大小仍然仅为128b）。内部循环只有20条指令，其中只有13条是FPU ALU指令，需要Intel SnB系列的端口0/1。即使使用基线SSE2，它也可以编写出体面的代码：无FMA，并且需要shufps的广播负载，但是那些不竞争带有add / mul的执行单元。
ax[i + 0..3]我也尝试过horizontal sum, which takes extra code，但IDK是否会更快。请注意，对于with the loops interchanged, but that used vector loads in the inner loop, with stride = 16B，gcc和clang会将除法转换为for Haswell with gcc, using FMA +牛顿迭代。 （即使出于独立原因，它们也出于某些原因不会将#include <immintrin.h>
void ffunc_sse128(float *restrict ax, float *restrict ay, float *restrict az,
const float *x, const float *y, const float *z,
int N, float eps, const float *restrict m)
{
for (int i = 0; i < N; i+=4) {
__m128 xi_v = _mm_load_ps(&x[i]);
__m128 yi_v = _mm_load_ps(&y[i]);
__m128 zi_v = _mm_load_ps(&z[i]);
// vector accumulators for ax[i + 0..3] etc.
__m128 axi_v = _mm_setzero_ps();
__m128 ayi_v = _mm_setzero_ps();
__m128 azi_v = _mm_setzero_ps();
// AVX broadcast-loads are as cheap as normal loads,
// and data-reuse meant that stand-alone load instructions were used anyway.
// so we're not even losing out on folding loads into other insns
// An inner-loop stride of only 4B is a huge win if memory / cache bandwidth is a bottleneck
// even without AVX, the shufps instructions are cheap,
// and don't compete with add/mul for execution units on Intel
for (int j = 0; j < N; j++) {
__m128 xj_v = _mm_set1_ps(x[j]);
__m128 rx_v = _mm_sub_ps(xj_v, xi_v);
__m128 yj_v = _mm_set1_ps(y[j]);
__m128 ry_v = _mm_sub_ps(yj_v, yi_v);
__m128 zj_v = _mm_set1_ps(z[j]);
__m128 rz_v = _mm_sub_ps(zj_v, zi_v);
__m128 mj_v = _mm_set1_ps(m[j]); // do the load early
// sum of squared differences
__m128 r2_v = _mm_set1_ps(eps) + rx_v*rx_v + ry_v*ry_v + rz_v*rz_v; // GNU extension
/* __m128 r2_v = _mm_add_ps(_mm_set1_ps(eps), _mm_mul_ps(rx_v, rx_v));
r2_v = _mm_add_ps(r2_v, _mm_mul_ps(ry_v, ry_v));
r2_v = _mm_add_ps(r2_v, _mm_mul_ps(rz_v, rz_v));
*/
// rsqrt and a Newton-Raphson iteration might have lower latency
// but there's enough surrounding instructions and cross-iteration parallelism
// that the single-uop sqrtps and divps instructions prob. aren't be a bottleneck
__m128 r2sqrt = _mm_sqrt_ps(r2_v);
__m128 r6sqrt = _mm_mul_ps(r2_v, r2sqrt); // v^(3/2) = sqrt(v)^3 = sqrt(v)*v
__m128 s_v = _mm_div_ps(mj_v, r6sqrt);
__m128 srx_v = _mm_mul_ps(s_v, rx_v);
__m128 sry_v = _mm_mul_ps(s_v, ry_v);
__m128 srz_v = _mm_mul_ps(s_v, rz_v);
axi_v = _mm_add_ps(axi_v, srx_v);
ayi_v = _mm_add_ps(ayi_v, sry_v);
azi_v = _mm_add_ps(azi_v, srz_v);
}
_mm_store_ps(&ax[i], axi_v);
_mm_store_ps(&ay[i], ayi_v);
_mm_store_ps(&az[i], azi_v);
}
}
转换为a version with rcpps + Newton）。 clang做的更好，将FMA用于迭代步骤。
-ffast-math