Skip to content

AVX2: eight-wide f32 BODFT forward/inverse combines, FMA intrinsics and SIMD tail scheduling improvements #73

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
falseywinchnet merged 2 commits into from
Jun 27, 2026
Merged

AVX2: eight-wide f32 BODFT forward/inverse combines, FMA intrinsics and SIMD tail scheduling improvements #73

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
83baa64
Document Bruun rewrite benchmark gates
falseywinchnet Jun 27, 2026
356ea25
Merge branch 'main' into fwysqa-main/optimize-performance-with-avx2/sse2
falseywinchnet Jun 27, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions src/detail/bodft_avx2_notes.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,3 +37,9 @@
- Partner bins are lane-reversed in registers with an AVX2 permute before storing, preserving the scalar/128-bit conjugate-pair layout while doubling the 128-bit float combine width on AVX2 hosts.
- The existing 128-bit SSE2/NEON float combine remains the fallback for non-AVX2 builds and for tiny levels that cannot fill an eight-position vector.
- Correctness was checked with the regular CTest suite and the standalone BODFT benchmark on an AVX2/FMA build; a separate SSE2/no-AVX build was also checked to keep the fallback path intact.

## Float inverse combine pass

- Added an eight-position AVX2/FMA single-precision inverse combine path that mirrors the existing double inverse algebra and reuses the float AVX2 deinterleave/interleave helpers.
- The path loads k/k+M and conjugate-partner blocks, reverses partner lanes in registers, reconstructs the four child spectra, applies conjugate twiddle rotations with `vfmaddps`/`vfmsubps`, and stores c0..c3 as packed complex children.
- In this container, a same-command BODFT benchmark comparison against the previous commit at N=65536 improved float inverse from about 1.13 ms to about 0.78 ms; N=1048576 improved from about 25.6 ms to about 18.4 ms. Forward timings remained noisy and mostly unchanged because this patch targets inverse combine work.
75 changes: 75 additions & 0 deletions src/detail/bodft_kernel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -239,6 +239,70 @@ static inline int bodft_combine_fwd_avx2_f32(const complex_f32_t* RESTRICT tab,
return k;
}


static inline int bodft_combine_inv_avx2_f32(const complex_f32_t* RESTRICT tab,
const complex_f32_t* RESTRICT tab2,
const complex_f32_t* RESTRICT tab3,
const complex_f32_t* RESTRICT in,
complex_f32_t* RESTRICT c0,
complex_f32_t* RESTRICT c1,
complex_f32_t* RESTRICT c2,
complex_f32_t* RESTRICT c3,
int M, int half) {
const __m256 halfv = _mm256_set1_ps(0.5f);
const __m256 zero = _mm256_setzero_ps();
int k = 0;
for (; k + 8 <= half; k += 8) {
__m256 tre, tim, t2re, t2im, t3re, t3im;
__m256 ykre, ykim, ykMre, ykMim;
__m256 ykpre, ykpim, ykpMre, ykpMim;
bodft_uzp8_ps(tab, k, tre, tim);
bodft_uzp8_ps(tab2, k, t2re, t2im);
bodft_uzp8_ps(tab3, k, t3re, t3im);
bodft_uzp8_ps(in, k, ykre, ykim);
bodft_uzp8_ps(in, k + M, ykMre, ykMim);
bodft_uzp8_ps(in, M - 8 - k, ykpre, ykpim);
bodft_uzp8_ps(in, M - 8 - k + M, ykpMre, ykpMim);
ykpre = bodft_rev8_ps(ykpre);
ykpim = bodft_rev8_ps(ykpim);
ykpMre = bodft_rev8_ps(ykpMre);
ykpMim = bodft_rev8_ps(ykpMim);

const __m256 e0re = _mm256_mul_ps(halfv, _mm256_add_ps(ykre, ykpMre));
const __m256 e0im = _mm256_mul_ps(halfv, _mm256_sub_ps(ykim, ykpMim));
const __m256 o0re = _mm256_mul_ps(halfv, _mm256_sub_ps(ykre, ykpMre));
const __m256 o0im = _mm256_mul_ps(halfv, _mm256_add_ps(ykim, ykpMim));
const __m256 e1re = _mm256_mul_ps(halfv, _mm256_add_ps(ykMre, ykpre));
const __m256 e1im = _mm256_mul_ps(halfv, _mm256_sub_ps(ykMim, ykpim));
const __m256 diff_re = _mm256_sub_ps(ykMre, ykpre);
const __m256 diff_im = _mm256_add_ps(ykMim, ykpim);
const __m256 o1re = _mm256_mul_ps(halfv, _mm256_sub_ps(zero, diff_im));
const __m256 o1im = _mm256_mul_ps(halfv, diff_re);

const __m256 b0re = _mm256_mul_ps(halfv, _mm256_add_ps(e0re, e1re));
const __m256 b0im = _mm256_mul_ps(halfv, _mm256_add_ps(e0im, e1im));
const __m256 b2re = _mm256_mul_ps(halfv, _mm256_sub_ps(e0re, e1re));
const __m256 b2im = _mm256_mul_ps(halfv, _mm256_sub_ps(e0im, e1im));
const __m256 b1re = _mm256_mul_ps(halfv, _mm256_add_ps(o0re, o1re));
const __m256 b1im = _mm256_mul_ps(halfv, _mm256_add_ps(o0im, o1im));
const __m256 b3re = _mm256_mul_ps(halfv, _mm256_sub_ps(o0re, o1re));
const __m256 b3im = _mm256_mul_ps(halfv, _mm256_sub_ps(o0im, o1im));

const __m256 c1re = _mm256_fmadd_ps(tim, b1im, _mm256_mul_ps(tre, b1re));
const __m256 c1im = _mm256_fmsub_ps(tre, b1im, _mm256_mul_ps(tim, b1re));
const __m256 c2re = _mm256_fmadd_ps(t2im, b2im, _mm256_mul_ps(t2re, b2re));
const __m256 c2im = _mm256_fmsub_ps(t2re, b2im, _mm256_mul_ps(t2im, b2re));
const __m256 c3re = _mm256_fmadd_ps(t3im, b3im, _mm256_mul_ps(t3re, b3re));
const __m256 c3im = _mm256_fmsub_ps(t3re, b3im, _mm256_mul_ps(t3im, b3re));

bodft_store8_ps(c0, k, b0re, b0im);
bodft_store8_ps(c1, k, c1re, c1im);
bodft_store8_ps(c2, k, c2re, c2im);
bodft_store8_ps(c3, k, c3re, c3im);
}
return k;
}

static inline int bodft_combine_inv_avx2_f64(const complex_t* RESTRICT tab,
const complex_t* RESTRICT tab2,
const complex_t* RESTRICT tab3,
Expand Down Expand Up @@ -537,6 +601,17 @@ static inline void combine_inv(const CT* RESTRICT tab, const CT* RESTRICT tab2,
reinterpret_cast<complex_t*>(c2),
reinterpret_cast<complex_t*>(c3), M, half);
}
if constexpr (sizeof(RT) == 4) {
k = bodft_combine_inv_avx2_f32(
reinterpret_cast<const complex_f32_t*>(tab),
reinterpret_cast<const complex_f32_t*>(tab2),
reinterpret_cast<const complex_f32_t*>(tab3),
reinterpret_cast<const complex_f32_t*>(in),
reinterpret_cast<complex_f32_t*>(c0),
reinterpret_cast<complex_f32_t*>(c1),
reinterpret_cast<complex_f32_t*>(c2),
reinterpret_cast<complex_f32_t*>(c3), M, half);
}
#endif
const RT h = static_cast<RT>(0.5);
for (; k < half; ++k) {
Expand Down
115 changes: 70 additions & 45 deletions src/detail/bruun_kernel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,8 +59,13 @@ typedef __m128d bruun_v2;
# define V2_SUB(a, b) _mm_sub_pd((a), (b))
# define V2_MUL(a, b) _mm_mul_pd((a), (b))
# define V2_DIV(a, b) _mm_div_pd((a), (b))
# define V2_MADD(a, b, c) V2_ADD((a), V2_MUL((b), (c)))
# define V2_MSUB(a, b, c) V2_SUB((a), V2_MUL((b), (c)))
# if defined(__FMA__)
# define V2_MADD(a, b, c) _mm_fmadd_pd((b), (c), (a))
# define V2_MSUB(a, b, c) _mm_fnmadd_pd((b), (c), (a))
# else
# define V2_MADD(a, b, c) V2_ADD((a), V2_MUL((b), (c)))
# define V2_MSUB(a, b, c) V2_SUB((a), V2_MUL((b), (c)))
# endif
# define V2_SET1(x) _mm_set1_pd(x)
# define V2_SETLH(l, h) _mm_set_pd((h), (l))
# define V2_UNPLO(a, b) _mm_unpacklo_pd((a), (b))
Expand Down Expand Up @@ -109,8 +114,13 @@ typedef __m128 bruun_v4f;
# define V4F_ADD(a, b) _mm_add_ps((a), (b))
# define V4F_SUB(a, b) _mm_sub_ps((a), (b))
# define V4F_MUL(a, b) _mm_mul_ps((a), (b))
# define V4F_MADD(a, b, c) V4F_ADD((a), V4F_MUL((b), (c)))
# define V4F_MSUB(a, b, c) V4F_SUB((a), V4F_MUL((b), (c)))
# if defined(__FMA__)
# define V4F_MADD(a, b, c) _mm_fmadd_ps((b), (c), (a))
# define V4F_MSUB(a, b, c) _mm_fnmadd_ps((b), (c), (a))
# else
# define V4F_MADD(a, b, c) V4F_ADD((a), V4F_MUL((b), (c)))
# define V4F_MSUB(a, b, c) V4F_SUB((a), V4F_MUL((b), (c)))
# endif
# define V4F_SET1(x) _mm_set1_ps(x)
# define V4F_SET4(a,b,c,d) _mm_setr_ps((a), (b), (c), (d))
# define V4F_ZERO() _mm_setzero_ps()
Expand Down Expand Up @@ -680,7 +690,8 @@ static inline void binomial_fwd(double* RESTRICT v, int h) {
_mm256_storeu_pd(v + i, _mm256_add_pd(a, b));
_mm256_storeu_pd(v + h + i, _mm256_sub_pd(a, b));
}
#elif BRUUN_LEVEL == 1
#endif
#if BRUUN_LEVEL >= 1
for (; i + 1 < h; i += 2) {
const bruun_v2 a = V2_LD(v + i);
const bruun_v2 b = V2_LD(v + h + i);
Expand Down Expand Up @@ -716,7 +727,8 @@ static inline void binomial_oop(const double* RESTRICT in, double* RESTRICT v, i
_mm256_storeu_pd(v + i, _mm256_add_pd(a, b));
_mm256_storeu_pd(v + h + i, _mm256_sub_pd(a, b));
}
#elif BRUUN_LEVEL == 1
#endif
#if BRUUN_LEVEL >= 1
for (; i + 1 < h; i += 2) {
const bruun_v2 a = V2_LD(in + i);
const bruun_v2 b = V2_LD(in + h + i);
Expand Down Expand Up @@ -755,7 +767,8 @@ static inline void binomial_inv(double* RESTRICT v, int h) {
_mm256_storeu_pd(v + i, _mm256_mul_pd(half4, _mm256_add_pd(a, b)));
_mm256_storeu_pd(v + h + i, _mm256_mul_pd(half4, _mm256_sub_pd(a, b)));
}
#elif BRUUN_LEVEL == 1
#endif
#if BRUUN_LEVEL >= 1
const bruun_v2 half2 = V2_SET1(0.5);
for (; i + 1 < h; i += 2) {
const bruun_v2 a = V2_LD(v + i);
Expand Down Expand Up @@ -823,7 +836,8 @@ static inline void norm_q_fwd(double* RESTRICT p, int q, double c_scalar, double
_mm256_storeu_pd(B1p + n, _mm256_sub_pd(I, A1));
}
}
#elif BRUUN_LEVEL == 1
#endif
#if BRUUN_LEVEL >= 1
{
const bruun_v2 vc = V2_SET1(c_scalar);
const bruun_v2 vs = V2_SET1(s_scalar);
Expand Down Expand Up @@ -967,7 +981,8 @@ static inline void norm2_fused(double* RESTRICT p, int q,
_mm256_storeu_pd(B1 + qh + n, _mm256_sub_pd(I1, x0));
}
}
#elif BRUUN_LEVEL == 1
#endif
#if BRUUN_LEVEL >= 1
{
const bruun_v2 vc = V2_SET1(c), vs = V2_SET1(s);
const bruun_v2 vc0 = V2_SET1(c0), vs0 = V2_SET1(s0);
Expand Down Expand Up @@ -1107,7 +1122,8 @@ static inline void norm_q_inv(double* RESTRICT p, int q, double c_scalar, double
_mm256_storeu_pd(D1p + n, B1);
}
}
#elif BRUUN_LEVEL == 1
#endif
#if BRUUN_LEVEL >= 1
{
const bruun_v2 half = V2_SET1(0.5);
const bruun_v2 vc = V2_SET1(c_scalar);
Expand Down Expand Up @@ -1268,7 +1284,8 @@ static inline void norm2_inv_fused(double* RESTRICT p, int q,
_mm256_storeu_pd(B1 + qh + n, _mm256_fmsub_pd(vc, Ih, _mm256_mul_pd(vs, Rh)));
}
}
#elif BRUUN_LEVEL == 1
#endif
#if BRUUN_LEVEL >= 1
{
const bruun_v2 hf = V2_SET1(0.5);
const bruun_v2 vc = V2_SET1(c), vs = V2_SET1(s);
Expand Down Expand Up @@ -1532,14 +1549,14 @@ static inline void norm_q_fwd_f32(float* RESTRICT p, int q, float c_scalar, floa
const __m256 vs = _mm256_set1_ps(s_scalar);

for (; n + 7 < q; n += 8) {
const __m256 A0 = _mm256_loadu_ps(A0p + n);
const __m256 B0 = _mm256_loadu_ps(B0p + n);
const __m256 A1 = _mm256_loadu_ps(A1p + n);
const __m256 B1 = _mm256_loadu_ps(B1p + n);

const __m256 R = _mm256_fmsub_ps(vc, B0, _mm256_mul_ps(vs, B1));
const __m256 I = _mm256_fmadd_ps(vs, B0, _mm256_mul_ps(vc, B1));

const __m256 A0 = _mm256_loadu_ps(A0p + n);
const __m256 A1 = _mm256_loadu_ps(A1p + n);

_mm256_storeu_ps(A0p + n, _mm256_add_ps(A0, R));
_mm256_storeu_ps(B0p + n, _mm256_add_ps(A1, I));
_mm256_storeu_ps(A1p + n, _mm256_sub_ps(A0, R));
Expand All @@ -1553,14 +1570,14 @@ static inline void norm_q_fwd_f32(float* RESTRICT p, int q, float c_scalar, floa
const bruun_v4f vs = V4F_SET1(s_scalar);

for (; n + 3 < q; n += 4) {
const bruun_v4f A0 = V4F_LD(A0p + n);
const bruun_v4f B0 = V4F_LD(B0p + n);
const bruun_v4f A1 = V4F_LD(A1p + n);
const bruun_v4f B1 = V4F_LD(B1p + n);

const bruun_v4f R = V4F_MSUB(V4F_MUL(vc, B0), vs, B1);
const bruun_v4f I = V4F_MADD(V4F_MUL(vs, B0), vc, B1);

const bruun_v4f A0 = V4F_LD(A0p + n);
const bruun_v4f A1 = V4F_LD(A1p + n);

V4F_ST(A0p + n, V4F_ADD(A0, R));
V4F_ST(B0p + n, V4F_ADD(A1, I));
V4F_ST(A1p + n, V4F_SUB(A0, R));
Expand Down Expand Up @@ -1652,20 +1669,21 @@ static inline void norm2_fused_f32(float* RESTRICT p, int q,
const __m256 vc1 = _mm256_set1_ps(c1), vs1 = _mm256_set1_ps(s1);

for (; n + 7 < qh; n += 8) {
const __m256 a0n = _mm256_loadu_ps(A0 + n);
const __m256 a0h = _mm256_loadu_ps(A0 + qh + n);
const __m256 b0n = _mm256_loadu_ps(B0 + n);
const __m256 b0h = _mm256_loadu_ps(B0 + qh + n);
const __m256 a1n = _mm256_loadu_ps(A1 + n);
const __m256 a1h = _mm256_loadu_ps(A1 + qh + n);
const __m256 b1n = _mm256_loadu_ps(B1 + n);
const __m256 b0h = _mm256_loadu_ps(B0 + qh + n);
const __m256 b1h = _mm256_loadu_ps(B1 + qh + n);

const __m256 Rn = _mm256_fmsub_ps(vc, b0n, _mm256_mul_ps(vs, b1n));
const __m256 In = _mm256_fmadd_ps(vs, b0n, _mm256_mul_ps(vc, b1n));
const __m256 Rh = _mm256_fmsub_ps(vc, b0h, _mm256_mul_ps(vs, b1h));
const __m256 Ih = _mm256_fmadd_ps(vs, b0h, _mm256_mul_ps(vc, b1h));

const __m256 a0n = _mm256_loadu_ps(A0 + n);
const __m256 a0h = _mm256_loadu_ps(A0 + qh + n);
const __m256 a1n = _mm256_loadu_ps(A1 + n);
const __m256 a1h = _mm256_loadu_ps(A1 + qh + n);

const __m256 u0 = _mm256_add_ps(a0n, Rn);
const __m256 uh = _mm256_add_ps(a0h, Rh);
const __m256 w0 = _mm256_add_ps(a1n, In);
Expand Down Expand Up @@ -1698,20 +1716,21 @@ static inline void norm2_fused_f32(float* RESTRICT p, int q,
const bruun_v4f vc1 = V4F_SET1(c1), vs1 = V4F_SET1(s1);

for (; n + 3 < qh; n += 4) {
const bruun_v4f a0n = V4F_LD(A0 + n);
const bruun_v4f a0h = V4F_LD(A0 + qh + n);
const bruun_v4f b0n = V4F_LD(B0 + n);
const bruun_v4f b0h = V4F_LD(B0 + qh + n);
const bruun_v4f a1n = V4F_LD(A1 + n);
const bruun_v4f a1h = V4F_LD(A1 + qh + n);
const bruun_v4f b1n = V4F_LD(B1 + n);
const bruun_v4f b0h = V4F_LD(B0 + qh + n);
const bruun_v4f b1h = V4F_LD(B1 + qh + n);

const bruun_v4f Rn = V4F_MSUB(V4F_MUL(vc, b0n), vs, b1n);
const bruun_v4f In = V4F_MADD(V4F_MUL(vs, b0n), vc, b1n);
const bruun_v4f Rh = V4F_MSUB(V4F_MUL(vc, b0h), vs, b1h);
const bruun_v4f Ih = V4F_MADD(V4F_MUL(vs, b0h), vc, b1h);

const bruun_v4f a0n = V4F_LD(A0 + n);
const bruun_v4f a0h = V4F_LD(A0 + qh + n);
const bruun_v4f a1n = V4F_LD(A1 + n);
const bruun_v4f a1h = V4F_LD(A1 + qh + n);

const bruun_v4f u0 = V4F_ADD(a0n, Rn);
const bruun_v4f uh = V4F_ADD(a0h, Rh);
const bruun_v4f w0 = V4F_ADD(a1n, In);
Expand Down Expand Up @@ -1808,22 +1827,24 @@ static inline void norm_q_inv_f32(float* RESTRICT p, int q, float c_scalar, floa
#if BRUUN_LEVEL >= 2
{
const __m256 half = _mm256_set1_ps(0.5f);
const __m256 vc = _mm256_set1_ps(c_scalar);
const __m256 vs = _mm256_set1_ps(s_scalar);
const __m256 hvc = _mm256_set1_ps(0.5f * c_scalar);
const __m256 hvs = _mm256_set1_ps(0.5f * s_scalar);

for (; n + 7 < q; n += 8) {
const __m256 C0v = _mm256_loadu_ps(C0p + n);
const __m256 C1v = _mm256_loadu_ps(C1p + n);
const __m256 D0v = _mm256_loadu_ps(D0p + n);
const __m256 D1v = _mm256_loadu_ps(D1p + n);

const __m256 A0 = _mm256_mul_ps(half, _mm256_add_ps(C0v, D0v));
const __m256 R = _mm256_mul_ps(half, _mm256_sub_ps(C0v, D0v));
const __m256 I = _mm256_mul_ps(half, _mm256_add_ps(C1v, D1v));
const __m256 A1 = _mm256_mul_ps(half, _mm256_sub_ps(C1v, D1v));
const __m256 t0 = _mm256_add_ps(C0v, D0v);
const __m256 r = _mm256_sub_ps(C0v, D0v);
const __m256 i = _mm256_add_ps(C1v, D1v);
const __m256 t1 = _mm256_sub_ps(C1v, D1v);

const __m256 B0 = _mm256_fmadd_ps(vc, R, _mm256_mul_ps(vs, I));
const __m256 B1 = _mm256_fmsub_ps(vc, I, _mm256_mul_ps(vs, R));
const __m256 A0 = _mm256_mul_ps(half, t0);
const __m256 A1 = _mm256_mul_ps(half, t1);
const __m256 B0 = _mm256_fmadd_ps(hvc, r, _mm256_mul_ps(hvs, i));
const __m256 B1 = _mm256_fmsub_ps(hvc, i, _mm256_mul_ps(hvs, r));

_mm256_storeu_ps(C0p + n, A0);
_mm256_storeu_ps(C1p + n, B0);
Expand All @@ -1835,22 +1856,24 @@ static inline void norm_q_inv_f32(float* RESTRICT p, int q, float c_scalar, floa
#if BRUUN_LEVEL >= 1
{
const bruun_v4f half = V4F_SET1(0.5f);
const bruun_v4f vc = V4F_SET1(c_scalar);
const bruun_v4f vs = V4F_SET1(s_scalar);
const bruun_v4f hvc = V4F_SET1(0.5f * c_scalar);
const bruun_v4f hvs = V4F_SET1(0.5f * s_scalar);

for (; n + 3 < q; n += 4) {
const bruun_v4f C0v = V4F_LD(C0p + n);
const bruun_v4f C1v = V4F_LD(C1p + n);
const bruun_v4f D0v = V4F_LD(D0p + n);
const bruun_v4f D1v = V4F_LD(D1p + n);

const bruun_v4f A0 = V4F_MUL(half, V4F_ADD(C0v, D0v));
const bruun_v4f R = V4F_MUL(half, V4F_SUB(C0v, D0v));
const bruun_v4f I = V4F_MUL(half, V4F_ADD(C1v, D1v));
const bruun_v4f A1 = V4F_MUL(half, V4F_SUB(C1v, D1v));
const bruun_v4f t0 = V4F_ADD(C0v, D0v);
const bruun_v4f r = V4F_SUB(C0v, D0v);
const bruun_v4f i = V4F_ADD(C1v, D1v);
const bruun_v4f t1 = V4F_SUB(C1v, D1v);

const bruun_v4f B0 = V4F_MADD(V4F_MUL(vc, R), vs, I);
const bruun_v4f B1 = V4F_MSUB(V4F_MUL(vc, I), vs, R);
const bruun_v4f A0 = V4F_MUL(half, t0);
const bruun_v4f A1 = V4F_MUL(half, t1);
const bruun_v4f B0 = V4F_MADD(V4F_MUL(hvc, r), hvs, i);
const bruun_v4f B1 = V4F_MSUB(V4F_MUL(hvc, i), hvs, r);

V4F_ST(C0p + n, A0);
V4F_ST(C1p + n, B0);
Expand Down Expand Up @@ -2367,7 +2390,8 @@ class RFFT {
bool standard_output_uses_two_phase() const {
#if BRUUN_LEVEL >= 2
return N >= 8192;
#elif BRUUN_LEVEL == 1
#endif
#if BRUUN_LEVEL >= 1
return N > 1048576;
#else
return false;
Expand Down Expand Up @@ -3045,7 +3069,8 @@ class RFFT {
_mm256_storeu_pd(v + j, _mm256_fmadd_pd(hd, r, _mm256_mul_pd(ho, rs)));
}
}
#elif BRUUN_LEVEL == 1
#endif
#if BRUUN_LEVEL >= 1
for (; j + 1 < end; j += 2) {
const bruun_v2 r = V2_LD(v + j);
const bruun_v2 h = V2_LD(RF + j);
Expand Down
9 changes: 9 additions & 0 deletions src/detail/inverse_acceleration_notes.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,3 +41,12 @@ Wide-x86 schedule gate:
- Rechecked the scheduled inverse on Intel Xeon Platinum 8370C with `-O3 -mavx2 -mfma -mno-avx512*`.
- Rechecked the AVX-512 scheduled inverse on the same host with `-O3 -march=native`.
- The benchmark prints a direct `Inv_ns` column so inverse-only changes are visible instead of inferred from round-trip timing.

2026-06-27 x86 SIMD tail and f32 scheduling pass:
- x86 128-bit V2/V4F multiply-add macros now use accumulator-first FMA intrinsics when `__FMA__` is available (`_mm_fmadd_*` and `_mm_fnmadd_*`) instead of lowering to separate multiply plus add/subtract. The fallback SSE2 definitions are unchanged for non-FMA builds.
- Double streaming kernels now run their 128-bit V2 loops after the AVX2/AVX-512 loops, so one-to-three-element wide tails no longer drop straight to scalar when a two-lane chunk remains.
- The float forward kernels start the B0/B1 twiddle rotation before loading A0/A1 in the AVX2 and V4F loops, shortening live ranges around the first dependency frontier.
- `norm_q_inv_f32` folds the 0.5 scale into the twiddle constants for rotated B outputs in the AVX2 and V4F loops, keeping explicit half multiplies only for the A outputs that are stored directly.

2026-06-27 continuation:
- Added `kernel_sequentiality_flow.md` to describe forward/inverse dependency frontiers, SIMD tail flow, and where arithmetic can safely move across frontiers without increasing live pressure.
Loading
Loading