Оптимизация с использованием SSE Intrinsics
Я пытаюсь преобразовать цикл, который у меня есть, в свойства 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);
}
}
Я чувствую, что основная идея верна, однако где-то есть некоторые ошибки, так как полученный ответ неверен.
Заранее спасибо.
1 ответ
Я думаю, что ваши единственные ошибки - это простые опечатки, а не логические ошибки, см. Ниже.
Вы не можете просто использовать 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 обнаружило эту ошибку:
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
Как я уже сказал в своем ответе на ваш предыдущий вопрос, вам, вероятно, следует поменять местами циклы, поэтому те же топоры [i+0..3], ay[i+0..3] и az[i+0..3] используются, избегая этой загрузки / хранения.
Кроме того, если вы не собираетесь использовать rsqrtps + Ньютон-Рафсон, вы должны использовать преобразование, которое я указал в моем последнем ответе: делить m[j] от sqrt(k2) ^3, Там нет смысла делить 1.0f чем-то, используя divps, а затем умножение только один раз.
rsqrt на самом деле не может быть победой, потому что общая пропускная способность UOP может быть более узким местом, чем пропускная способность div / sqrt или задержки. три умножения + FMA + rsqrtps значительно больше мопов, чем sqrtps + divps. rsqrt более полезен для векторов AVX 256b, потому что модуль div / sqrt не является полноширинным на SnB для Broadwell. Skylake имеет задержку 12С sqrtps ymm То же, что и для xmm, но пропускная способность все же лучше для xmm (один на 3c вместо одного на 6c).
Clang и GCC оба использовали rsqrtps / rsqrtss при компиляции вашего кода с -ffast-math, (конечно, только лягушатник, использующий упакованную версию.)
Если вы не меняете местами петли, вы должны вручную поднять все, что зависит только от j из внутреннего цикла. Компиляторы, как правило, хороши в этом, но все же кажется хорошей идеей, чтобы исходный код отражал то, что вы ожидаете от компилятора. Это помогает "видеть", что на самом деле делает цикл.
Вот версия с некоторыми оптимизациями по сравнению с оригиналом:
Чтобы заставить gcc / clang объединить мул / аддит в FMA, я использовал
-ffp-contract=fast, Это получает инструкции FMA для высокой пропускной способности без использования-ffast-math, (Существует много параллелизма с тремя отдельными аккумуляторами, поэтому увеличенная задержка FMA по сравнению сaddpsне должно быть больно. Я ожидаю, что пропускная способность port0/1 является ограничивающим фактором.) Я думал, что gcc сделал это автоматически, но, похоже, что здесь без него-ffast-math,Обратите внимание, что v 3/2 = sqrt (v) 3 = sqrt (v) * v. Это имеет меньшую задержку и меньше инструкций.
Поменяйте местами циклы и используйте широковещательные нагрузки во внутреннем цикле для улучшения локальности (уменьшите требования к полосе пропускания на 4 или 8 с AVX). Каждая итерация внутреннего цикла считывает только 4B новых данных из каждого из четырех исходных массивов. (x,y,z и m). Таким образом, она активно использует каждую строку кэша, пока она горячая.
Использование широковещательных нагрузок во внутреннем цикле также означает, что мы накапливаем
ax[i + 0..3]параллельно, избегая необходимости горизонтальной суммы, которая требует дополнительного кода. (См. Предыдущую версию этого ответа для кода с замененными циклами, но который использовал векторные нагрузки во внутреннем цикле, с шагом = 16B.)
Он хорошо компилируется для Haswell с gcc, используя FMA. (Тем не менее, размер вектора только 128b, в отличие от версии 256b с автографом вектора clang). Внутренний цикл содержит только 20 инструкций, и только 13 из них являются инструкциями FPU ALU, которым требуется порт 0/1 в семействе Intel SnB. Он делает достойный код даже с базовым SSE2: нет FMA, и ему нужны shufps для широковещательных нагрузок, но они не конкурируют за исполнительные блоки с add/mul.
#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);
}
}
Я также попробовал версию с rcpps, но IDK, если она будет быстрее. Обратите внимание, что с -ffast-math, gcc и clang преобразуют деление в rcpps + итерация Ньютона. (Они почему-то не конвертируют 1.0f/sqrtf(x) в rsqrt + Ньютон даже в автономной функции). Clang делает лучше, используя FMA для шага итерации.
#define USE_RSQRT
#ifndef USE_RSQRT
// even with -mrecip=vec-sqrt after -ffast-math, this still does sqrt(v)*v, then rcpps
__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);
#else
__m128 r2isqrt = rsqrt_float4_single(r2_v);
// can't use the sqrt(v)*v trick, unless we either do normal sqrt first then rcpps
// or rsqrtps and rcpps. Maybe it's possible to do a Netwon Raphson iteration on that product
// instead of refining them both separately?
__m128 r6isqrt = r2isqrt * r2isqrt * r2isqrt;
__m128 s_v = _mm_mul_ps(mj_v, r6isqrt);
#endif