NumericAL
A header-only C++20 lexer built on REAL
Loading...
Searching...
No Matches
x86.hpp
Go to the documentation of this file.
1
11#ifndef NUMERICAL_X86_HPP
12#define NUMERICAL_X86_HPP
13
14#include <algorithm>
15#include <array>
16#include <cmath>
17#include <complex>
18#include <cstddef>
19#include <cstdint>
20
21#include "simd.hpp"
22
23#if defined(__x86_64__) || defined(_M_X64)
24# define NUMERICAL_X86_COMPILED // the x86-64 tiles are compiled in (tested with #ifdef only)
25# include <immintrin.h>
26# if defined(_MSC_VER) && !defined(__clang__)
27# include <intrin.h>
28# define NUMERICAL_TARGET(isa)
29# define NUMERICAL_UNROLL_TILE
30# define NUMERICAL_UNROLL_DEPTH
31# define NUMERICAL_UNROLL_DEPTH_AVX2
32# else
33# include <cpuid.h>
34# define NUMERICAL_TARGET(isa) __attribute__((target(isa))) // NOLINT(cppcoreguidelines-macro-usage): an attribute, which no function can stand for
35// A tile loop unrolled whole, so each accumulator stays in its register (rolled, GCC keeps them in memory).
36# define NUMERICAL_UNROLL_TILE _Pragma("GCC unroll 24")
37// The depth loop unrolled by two: half the loop's own instructions, which compete with the multiply-adds for
38// issue slots (an AVX-512 tile in the L1 cache: 103 GFLOP/s rolled, 115 unrolled, on one x86-64 core).
39# define NUMERICAL_UNROLL_DEPTH _Pragma("GCC unroll 2")
40// The AVX2 tile's depth loop unrolled by four, as OpenBLAS's Haswell kernel unrolls it (by eight): a product of order
41// 1000 on one AVX2 core, 50.5 → 52.9 GFLOP/s, by two; by eight no better.
42# define NUMERICAL_UNROLL_DEPTH_AVX2 _Pragma("GCC unroll 4")
43# endif
44#endif
45
46namespace numerical::detail {
47
48#ifdef NUMERICAL_X86_COMPILED
50 struct x86_features
51 {
52 bool avx2_fma {false};
53 bool avx512f {false};
54 };
55
57 [[nodiscard]] inline std::array<std::uint32_t, 4> cpuid(std::uint32_t leaf,
58 std::uint32_t sub) noexcept
59 {
60 std::array<std::uint32_t, 4> r {};
61# if defined(_MSC_VER) && !defined(__clang__)
62 std::array<int, 4> regs {};
63 __cpuidex(regs.data(), static_cast<int>(leaf), static_cast<int>(sub));
64 for (std::size_t i {0}; i < 4; ++i) {
65 r[i] = static_cast<std::uint32_t>(regs[i]);
66 }
67# else
68 __cpuid_count(leaf, sub, r[0], r[1], r[2], r[3]);
69# endif
70 return r;
71 }
72
74 [[nodiscard]] NUMERICAL_TARGET("xsave") inline std::uint64_t saved_state(std::uint32_t leaf1_ecx) noexcept
75 {
76 if ((leaf1_ecx & (1U << 27U)) == 0U) { // OSXSAVE
77 return 0;
78 }
79 return _xgetbv(0);
80 }
81
83 [[nodiscard]] constexpr x86_features decode_features(std::uint32_t leaf1_ecx,
84 std::uint32_t leaf7_ebx,
85 std::uint64_t state) noexcept
86 {
87 const bool avx {(leaf1_ecx & (1U << 28U)) != 0U};
88 const bool ymm {avx && (state & 0x6U) == 0x6U}; // the AVX instructions, and SSE and AVX state
89 const bool zmm {avx && (state & 0xe6U) == 0xe6U}; // and the opmask and both halves of the 512-bit registers
90 const bool fma {(leaf1_ecx & (1U << 12U)) != 0U};
91 const bool avx2 {(leaf7_ebx & (1U << 5U)) != 0U};
92 const bool avx512f {(leaf7_ebx & (1U << 16U)) != 0U};
93 return {.avx2_fma = ymm && avx2 && fma, .avx512f = zmm && avx512f};
94 }
95
97 [[nodiscard]] inline x86_features cpu_features() noexcept
98 {
99 static const x86_features features {[] {
100 const std::array<std::uint32_t, 4> leaf1 {cpuid(1, 0)};
101 // Leaf 7 only where the processor has it: past its last leaf, cpuid repeats
102 // another leaf's words.
103 const std::array<std::uint32_t, 4> leaf7 {cpuid(0, 0)[0] >= 7U ? cpuid(7, 0)
104 : std::array<std::uint32_t, 4> {}};
105 return decode_features(leaf1[2], leaf7[1], saved_state(leaf1[2]));
106 }()};
107 return features;
108 }
109
112 struct zmm_register
113 {
114 __m512d v;
115 };
116
118 struct ymm_register
119 {
120 __m256d v;
121 };
122
124 inline constexpr std::size_t avx512_mr {12};
125
127 inline constexpr std::size_t avx512_nr {16};
128
130 inline constexpr std::size_t avx2_mr {6};
131
133 inline constexpr std::size_t avx2_nr {8};
134
140 NUMERICAL_TARGET("avx512f") inline void avx512_full_tile_f64(std::size_t kc,
141 const double* ap,
142 const double* bp,
143 std::size_t ldb,
144 double* tile,
145 std::size_t ld) noexcept
146 {
147 std::array<zmm_register, avx512_mr * 2> acc {};
148 NUMERICAL_UNROLL_TILE
149 for (std::size_t i {0}; i < avx512_mr; ++i) {
150 acc[2 * i].v = _mm512_loadu_pd(tile + (i * ld));
151 acc[(2 * i) + 1].v = _mm512_loadu_pd(tile + (i * ld) + 8);
152 }
153 NUMERICAL_UNROLL_DEPTH
154 for (std::size_t p {0}; p < kc; ++p) {
155 const __m512d b0 {_mm512_loadu_pd(bp + (p * ldb))};
156 const __m512d b1 {_mm512_loadu_pd(bp + (p * ldb) + 8)};
157 NUMERICAL_UNROLL_TILE
158 for (std::size_t i {0}; i < avx512_mr; ++i) {
159 const __m512d a {_mm512_set1_pd(ap[(p * avx512_mr) + i])};
160 acc[2 * i].v = _mm512_fmadd_pd(a, b0, acc[2 * i].v);
161 acc[(2 * i) + 1].v = _mm512_fmadd_pd(a, b1, acc[(2 * i) + 1].v);
162 }
163 }
164 NUMERICAL_UNROLL_TILE
165 for (std::size_t i {0}; i < avx512_mr; ++i) {
166 _mm512_storeu_pd(tile + (i * ld), acc[2 * i].v);
167 _mm512_storeu_pd(tile + (i * ld) + 8, acc[(2 * i) + 1].v);
168 }
169 _mm256_zeroupper(); // the upper halves cleared: code past the tile may be SSE alone
170 }
171
174 inline void avx512_tile_f64(std::size_t kc,
175 const double* ap,
176 const double* bp,
177 std::size_t ldb,
178 double* c,
179 std::size_t ldc,
180 std::size_t rows,
181 std::size_t cols) noexcept
182 {
183 through_edge<avx512_mr, avx512_nr>(c, ldc, rows, cols, [&](double* tile, std::size_t ld) noexcept { avx512_full_tile_f64(kc, ap, bp, ldb, tile, ld); });
184 }
185
187 inline constexpr std::size_t avx2_prefetch_ahead {64};
188
193 NUMERICAL_TARGET("avx2,fma") inline void avx2_full_tile_f64(std::size_t kc,
194 const double* ap,
195 const double* bp,
196 std::size_t ldb,
197 double* tile,
198 std::size_t ld) noexcept
199 {
200 std::array<ymm_register, avx2_mr * 2> acc {};
201 NUMERICAL_UNROLL_TILE
202 for (std::size_t i {0}; i < avx2_mr; ++i) {
203 acc[2 * i].v = _mm256_loadu_pd(tile + (i * ld));
204 acc[(2 * i) + 1].v = _mm256_loadu_pd(tile + (i * ld) + 4);
205 }
206 NUMERICAL_UNROLL_DEPTH_AVX2
207 for (std::size_t p {0}; p < kc; ++p) {
208 // A and B a few steps ahead, into the L1 cache: the hardware's streams start a step late at each sliver (as
209 // OpenBLAS's Haswell kernel fetches both 512 bytes ahead; with the unrolling, order 1000 52.9 → 53.9 GFLOP/s).
210 _mm_prefetch(reinterpret_cast<const char*>(bp + (p * ldb) + avx2_prefetch_ahead), _MM_HINT_T0); // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the intrinsic takes a char pointer
211 _mm_prefetch(reinterpret_cast<const char*>(ap + (p * avx2_mr) + avx2_prefetch_ahead), _MM_HINT_T0); // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the same
212 const __m256d b0 {_mm256_loadu_pd(bp + (p * ldb))};
213 const __m256d b1 {_mm256_loadu_pd(bp + (p * ldb) + 4)};
214 NUMERICAL_UNROLL_TILE
215 for (std::size_t i {0}; i < avx2_mr; ++i) {
216 const __m256d a {_mm256_broadcast_sd(ap + (p * avx2_mr) + i)};
217 acc[2 * i].v = _mm256_fmadd_pd(a, b0, acc[2 * i].v);
218 acc[(2 * i) + 1].v = _mm256_fmadd_pd(a, b1, acc[(2 * i) + 1].v);
219 }
220 }
221 NUMERICAL_UNROLL_TILE
222 for (std::size_t i {0}; i < avx2_mr; ++i) {
223 _mm256_storeu_pd(tile + (i * ld), acc[2 * i].v);
224 _mm256_storeu_pd(tile + (i * ld) + 4, acc[(2 * i) + 1].v);
225 }
226 _mm256_zeroupper(); // the upper halves cleared: code past the tile may be SSE alone
227 }
228
231 inline void avx2_tile_f64(std::size_t kc,
232 const double* ap,
233 const double* bp,
234 std::size_t ldb,
235 double* c,
236 std::size_t ldc,
237 std::size_t rows,
238 std::size_t cols) noexcept
239 {
240 through_edge<avx2_mr, avx2_nr>(c, ldc, rows, cols, [&](double* tile, std::size_t ld) noexcept { avx2_full_tile_f64(kc, ap, bp, ldb, tile, ld); });
241 }
242
244 struct zmm_register_f32
245 {
246 __m512 v;
247 };
248
250 struct ymm_register_f32
251 {
252 __m256 v;
253 };
254
256 inline constexpr std::size_t avx512_nr_f32 {32};
257
259 inline constexpr std::size_t avx2_nr_f32 {16};
260
262 NUMERICAL_TARGET("avx512f") inline void avx512_full_tile_f32(std::size_t kc,
263 const float* ap,
264 const float* bp,
265 std::size_t ldb,
266 float* tile,
267 std::size_t ld) noexcept
268 {
269 std::array<zmm_register_f32, avx512_mr * 2> acc {};
270 NUMERICAL_UNROLL_TILE
271 for (std::size_t i {0}; i < avx512_mr; ++i) {
272 acc[2 * i].v = _mm512_loadu_ps(tile + (i * ld));
273 acc[(2 * i) + 1].v = _mm512_loadu_ps(tile + (i * ld) + 16);
274 }
275 NUMERICAL_UNROLL_DEPTH
276 for (std::size_t p {0}; p < kc; ++p) {
277 const __m512 b0 {_mm512_loadu_ps(bp + (p * ldb))};
278 const __m512 b1 {_mm512_loadu_ps(bp + (p * ldb) + 16)};
279 NUMERICAL_UNROLL_TILE
280 for (std::size_t i {0}; i < avx512_mr; ++i) {
281 const __m512 a {_mm512_set1_ps(ap[(p * avx512_mr) + i])};
282 acc[2 * i].v = _mm512_fmadd_ps(a, b0, acc[2 * i].v);
283 acc[(2 * i) + 1].v = _mm512_fmadd_ps(a, b1, acc[(2 * i) + 1].v);
284 }
285 }
286 NUMERICAL_UNROLL_TILE
287 for (std::size_t i {0}; i < avx512_mr; ++i) {
288 _mm512_storeu_ps(tile + (i * ld), acc[2 * i].v);
289 _mm512_storeu_ps(tile + (i * ld) + 16, acc[(2 * i) + 1].v);
290 }
291 _mm256_zeroupper(); // the upper halves cleared: code past the tile may be SSE alone
292 }
293
295 inline void avx512_tile_f32(std::size_t kc,
296 const float* ap,
297 const float* bp,
298 std::size_t ldb,
299 float* c,
300 std::size_t ldc,
301 std::size_t rows,
302 std::size_t cols) noexcept
303 {
304 through_edge<avx512_mr, avx512_nr_f32>(c, ldc, rows, cols, [&](float* tile, std::size_t ld) noexcept { avx512_full_tile_f32(kc, ap, bp, ldb, tile, ld); });
305 }
306
308 NUMERICAL_TARGET("avx2,fma") inline void avx2_full_tile_f32(std::size_t kc,
309 const float* ap,
310 const float* bp,
311 std::size_t ldb,
312 float* tile,
313 std::size_t ld) noexcept
314 {
315 std::array<ymm_register_f32, avx2_mr * 2> acc {};
316 NUMERICAL_UNROLL_TILE
317 for (std::size_t i {0}; i < avx2_mr; ++i) {
318 acc[2 * i].v = _mm256_loadu_ps(tile + (i * ld));
319 acc[(2 * i) + 1].v = _mm256_loadu_ps(tile + (i * ld) + 8);
320 }
321 NUMERICAL_UNROLL_DEPTH_AVX2
322 for (std::size_t p {0}; p < kc; ++p) {
323 _mm_prefetch(reinterpret_cast<const char*>(bp + (p * ldb) + (2 * avx2_prefetch_ahead)), _MM_HINT_T0); // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the intrinsic takes a char pointer
324 _mm_prefetch(reinterpret_cast<const char*>(ap + (p * avx2_mr) + (2 * avx2_prefetch_ahead)), _MM_HINT_T0); // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the same
325 const __m256 b0 {_mm256_loadu_ps(bp + (p * ldb))};
326 const __m256 b1 {_mm256_loadu_ps(bp + (p * ldb) + 8)};
327 NUMERICAL_UNROLL_TILE
328 for (std::size_t i {0}; i < avx2_mr; ++i) {
329 const __m256 a {_mm256_broadcast_ss(ap + (p * avx2_mr) + i)};
330 acc[2 * i].v = _mm256_fmadd_ps(a, b0, acc[2 * i].v);
331 acc[(2 * i) + 1].v = _mm256_fmadd_ps(a, b1, acc[(2 * i) + 1].v);
332 }
333 }
334 NUMERICAL_UNROLL_TILE
335 for (std::size_t i {0}; i < avx2_mr; ++i) {
336 _mm256_storeu_ps(tile + (i * ld), acc[2 * i].v);
337 _mm256_storeu_ps(tile + (i * ld) + 8, acc[(2 * i) + 1].v);
338 }
339 _mm256_zeroupper(); // the upper halves cleared: code past the tile may be SSE alone
340 }
341
343 inline void avx2_tile_f32(std::size_t kc,
344 const float* ap,
345 const float* bp,
346 std::size_t ldb,
347 float* c,
348 std::size_t ldc,
349 std::size_t rows,
350 std::size_t cols) noexcept
351 {
352 through_edge<avx2_mr, avx2_nr_f32>(c, ldc, rows, cols, [&](float* tile, std::size_t ld) noexcept { avx2_full_tile_f32(kc, ap, bp, ldb, tile, ld); });
353 }
354
356 inline constexpr std::size_t avx2_complex_mr {3};
357
359 inline constexpr std::size_t avx2_complex_nr {4};
360
368 NUMERICAL_TARGET("avx2,fma") inline void avx2_full_tile_c64(std::size_t kc,
369 const std::complex<double>* ap,
370 const std::complex<double>* bp,
371 std::size_t ldb,
372 std::complex<double>* tile,
373 std::size_t ld) noexcept
374 {
375 const auto* const a {reinterpret_cast<const double*>(ap)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
376 const auto* const b {reinterpret_cast<const double*>(bp)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the same
377 auto* const c {reinterpret_cast<double*>(tile)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the same
378 std::array<ymm_register, avx2_complex_mr * 4> acc {};
379 NUMERICAL_UNROLL_DEPTH
380 for (std::size_t p {0}; p < kc; ++p) {
381 const __m256d b0 {_mm256_loadu_pd(b + (2 * p * ldb))};
382 const __m256d b1 {_mm256_loadu_pd(b + (2 * p * ldb) + 4)};
383 NUMERICAL_UNROLL_TILE
384 for (std::size_t i {0}; i < avx2_complex_mr; ++i) {
385 const __m256d re {_mm256_broadcast_sd(a + (2 * ((p * avx2_complex_mr) + i)))};
386 const __m256d im {_mm256_broadcast_sd(a + (2 * ((p * avx2_complex_mr) + i)) + 1)};
387 acc[4 * i].v = _mm256_fmadd_pd(re, b0, acc[4 * i].v);
388 acc[(4 * i) + 1].v = _mm256_fmadd_pd(re, b1, acc[(4 * i) + 1].v);
389 acc[(4 * i) + 2].v = _mm256_fmadd_pd(im, b0, acc[(4 * i) + 2].v);
390 acc[(4 * i) + 3].v = _mm256_fmadd_pd(im, b1, acc[(4 * i) + 3].v);
391 }
392 }
393 NUMERICAL_UNROLL_TILE
394 for (std::size_t i {0}; i < avx2_complex_mr; ++i) {
395 double* const row {c + (2 * i * ld)};
396 const __m256d lo {_mm256_addsub_pd(acc[4 * i].v, _mm256_permute_pd(acc[(4 * i) + 2].v, 0b0101))};
397 const __m256d hi {_mm256_addsub_pd(acc[(4 * i) + 1].v, _mm256_permute_pd(acc[(4 * i) + 3].v, 0b0101))};
398 _mm256_storeu_pd(row, _mm256_add_pd(_mm256_loadu_pd(row), lo));
399 _mm256_storeu_pd(row + 4, _mm256_add_pd(_mm256_loadu_pd(row + 4), hi));
400 }
401 _mm256_zeroupper(); // the upper halves cleared: code past the tile may be SSE alone
402 }
403
406 inline void avx2_tile_c64(std::size_t kc,
407 const std::complex<double>* ap,
408 const std::complex<double>* bp,
409 std::size_t ldb,
410 std::complex<double>* c,
411 std::size_t ldc,
412 std::size_t rows,
413 std::size_t cols) noexcept
414 {
415 through_edge<avx2_complex_mr, avx2_complex_nr>(c, ldc, rows, cols, [&](std::complex<double>* tile, std::size_t ld) noexcept {
416 avx2_full_tile_c64(kc, ap, bp, ldb, tile, ld);
417 });
418 }
419
421 inline constexpr std::size_t avx512_complex_mr {6};
422
424 inline constexpr std::size_t avx512_complex_nr {8};
425
431 NUMERICAL_TARGET("avx512f") inline void avx512_full_tile_c64(std::size_t kc,
432 const std::complex<double>* ap,
433 const std::complex<double>* bp,
434 std::size_t ldb,
435 std::complex<double>* tile,
436 std::size_t ld) noexcept
437 {
438 const auto* const a {reinterpret_cast<const double*>(ap)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee above
439 const auto* const b {reinterpret_cast<const double*>(bp)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the same
440 auto* const c {reinterpret_cast<double*>(tile)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the same
441 std::array<zmm_register, avx512_complex_mr * 4> acc {};
442 NUMERICAL_UNROLL_DEPTH
443 for (std::size_t p {0}; p < kc; ++p) {
444 const __m512d b0 {_mm512_loadu_pd(b + (2 * p * ldb))};
445 const __m512d b1 {_mm512_loadu_pd(b + (2 * p * ldb) + 8)};
446 NUMERICAL_UNROLL_TILE
447 for (std::size_t i {0}; i < avx512_complex_mr; ++i) {
448 const __m512d re {_mm512_set1_pd(a[2 * ((p * avx512_complex_mr) + i)])};
449 const __m512d im {_mm512_set1_pd(a[(2 * ((p * avx512_complex_mr) + i)) + 1])};
450 acc[4 * i].v = _mm512_fmadd_pd(re, b0, acc[4 * i].v);
451 acc[(4 * i) + 1].v = _mm512_fmadd_pd(re, b1, acc[(4 * i) + 1].v);
452 acc[(4 * i) + 2].v = _mm512_fmadd_pd(im, b0, acc[(4 * i) + 2].v);
453 acc[(4 * i) + 3].v = _mm512_fmadd_pd(im, b1, acc[(4 * i) + 3].v);
454 }
455 }
456 const __m512d signs {_mm512_set_pd(1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0)}; // −1 on the real parts
457 NUMERICAL_UNROLL_TILE
458 for (std::size_t i {0}; i < avx512_complex_mr; ++i) {
459 double* const row {c + (2 * i * ld)};
460 const __m512d lo {_mm512_fmadd_pd(_mm512_permute_pd(acc[(4 * i) + 2].v, 0x55), signs, acc[4 * i].v)};
461 const __m512d hi {_mm512_fmadd_pd(_mm512_permute_pd(acc[(4 * i) + 3].v, 0x55), signs, acc[(4 * i) + 1].v)};
462 _mm512_storeu_pd(row, _mm512_add_pd(_mm512_loadu_pd(row), lo));
463 _mm512_storeu_pd(row + 8, _mm512_add_pd(_mm512_loadu_pd(row + 8), hi));
464 }
465 }
466
468 inline void avx512_tile_c64(std::size_t kc,
469 const std::complex<double>* ap,
470 const std::complex<double>* bp,
471 std::size_t ldb,
472 std::complex<double>* c,
473 std::size_t ldc,
474 std::size_t rows,
475 std::size_t cols) noexcept
476 {
477 through_edge<avx512_complex_mr, avx512_complex_nr>(c, ldc, rows, cols, [&](std::complex<double>* tile, std::size_t ld) noexcept {
478 avx512_full_tile_c64(kc, ap, bp, ldb, tile, ld);
479 });
480 }
481
483 inline constexpr std::size_t avx2_complex_nr_f32 {8};
484
489 NUMERICAL_TARGET("avx2,fma") inline void avx2_full_tile_c32(std::size_t kc,
490 const std::complex<float>* ap,
491 const std::complex<float>* bp,
492 std::size_t ldb,
493 std::complex<float>* tile,
494 std::size_t ld) noexcept
495 {
496 const auto* const a {reinterpret_cast<const float*>(ap)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee of [complex.numbers]
497 const auto* const b {reinterpret_cast<const float*>(bp)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the same
498 auto* const c {reinterpret_cast<float*>(tile)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the same
499 std::array<ymm_register_f32, avx2_complex_mr * 4> acc {};
500 NUMERICAL_UNROLL_DEPTH
501 for (std::size_t p {0}; p < kc; ++p) {
502 const __m256 b0 {_mm256_loadu_ps(b + (2 * p * ldb))};
503 const __m256 b1 {_mm256_loadu_ps(b + (2 * p * ldb) + 8)};
504 NUMERICAL_UNROLL_TILE
505 for (std::size_t i {0}; i < avx2_complex_mr; ++i) {
506 const __m256 re {_mm256_broadcast_ss(a + (2 * ((p * avx2_complex_mr) + i)))};
507 const __m256 im {_mm256_broadcast_ss(a + (2 * ((p * avx2_complex_mr) + i)) + 1)};
508 acc[4 * i].v = _mm256_fmadd_ps(re, b0, acc[4 * i].v);
509 acc[(4 * i) + 1].v = _mm256_fmadd_ps(re, b1, acc[(4 * i) + 1].v);
510 acc[(4 * i) + 2].v = _mm256_fmadd_ps(im, b0, acc[(4 * i) + 2].v);
511 acc[(4 * i) + 3].v = _mm256_fmadd_ps(im, b1, acc[(4 * i) + 3].v);
512 }
513 }
514 NUMERICAL_UNROLL_TILE
515 for (std::size_t i {0}; i < avx2_complex_mr; ++i) {
516 float* const row {c + (2 * i * ld)};
517 const __m256 lo {_mm256_addsub_ps(acc[4 * i].v, _mm256_permute_ps(acc[(4 * i) + 2].v, 0xB1))};
518 const __m256 hi {_mm256_addsub_ps(acc[(4 * i) + 1].v, _mm256_permute_ps(acc[(4 * i) + 3].v, 0xB1))};
519 _mm256_storeu_ps(row, _mm256_add_ps(_mm256_loadu_ps(row), lo));
520 _mm256_storeu_ps(row + 8, _mm256_add_ps(_mm256_loadu_ps(row + 8), hi));
521 }
522 _mm256_zeroupper(); // the upper halves cleared: code past the tile may be SSE alone
523 }
524
526 inline void avx2_tile_c32(std::size_t kc,
527 const std::complex<float>* ap,
528 const std::complex<float>* bp,
529 std::size_t ldb,
530 std::complex<float>* c,
531 std::size_t ldc,
532 std::size_t rows,
533 std::size_t cols) noexcept
534 {
535 through_edge<avx2_complex_mr, avx2_complex_nr_f32>(c, ldc, rows, cols, [&](std::complex<float>* tile, std::size_t ld) noexcept {
536 avx2_full_tile_c32(kc, ap, bp, ldb, tile, ld);
537 });
538 }
539
541 inline constexpr std::size_t avx512_complex_nr_f32 {16};
542
545 NUMERICAL_TARGET("avx512f") inline void avx512_full_tile_c32(std::size_t kc,
546 const std::complex<float>* ap,
547 const std::complex<float>* bp,
548 std::size_t ldb,
549 std::complex<float>* tile,
550 std::size_t ld) noexcept
551 {
552 const auto* const a {reinterpret_cast<const float*>(ap)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee of [complex.numbers]
553 const auto* const b {reinterpret_cast<const float*>(bp)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the same
554 auto* const c {reinterpret_cast<float*>(tile)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the same
555 std::array<zmm_register_f32, avx512_complex_mr * 4> acc {};
556 NUMERICAL_UNROLL_DEPTH
557 for (std::size_t p {0}; p < kc; ++p) {
558 const __m512 b0 {_mm512_loadu_ps(b + (2 * p * ldb))};
559 const __m512 b1 {_mm512_loadu_ps(b + (2 * p * ldb) + 16)};
560 NUMERICAL_UNROLL_TILE
561 for (std::size_t i {0}; i < avx512_complex_mr; ++i) {
562 const __m512 re {_mm512_set1_ps(a[2 * ((p * avx512_complex_mr) + i)])};
563 const __m512 im {_mm512_set1_ps(a[(2 * ((p * avx512_complex_mr) + i)) + 1])};
564 acc[4 * i].v = _mm512_fmadd_ps(re, b0, acc[4 * i].v);
565 acc[(4 * i) + 1].v = _mm512_fmadd_ps(re, b1, acc[(4 * i) + 1].v);
566 acc[(4 * i) + 2].v = _mm512_fmadd_ps(im, b0, acc[(4 * i) + 2].v);
567 acc[(4 * i) + 3].v = _mm512_fmadd_ps(im, b1, acc[(4 * i) + 3].v);
568 }
569 }
570 const __m512 signs {_mm512_set_ps(1.0F, -1.0F, 1.0F, -1.0F, 1.0F, -1.0F, 1.0F, -1.0F, 1.0F, -1.0F, 1.0F, -1.0F, 1.0F, -1.0F, 1.0F, -1.0F)}; // −1 on the real parts
571 NUMERICAL_UNROLL_TILE
572 for (std::size_t i {0}; i < avx512_complex_mr; ++i) {
573 float* const row {c + (2 * i * ld)};
574 const __m512 lo {_mm512_fmadd_ps(_mm512_permute_ps(acc[(4 * i) + 2].v, 0xB1), signs, acc[4 * i].v)};
575 const __m512 hi {_mm512_fmadd_ps(_mm512_permute_ps(acc[(4 * i) + 3].v, 0xB1), signs, acc[(4 * i) + 1].v)};
576 _mm512_storeu_ps(row, _mm512_add_ps(_mm512_loadu_ps(row), lo));
577 _mm512_storeu_ps(row + 16, _mm512_add_ps(_mm512_loadu_ps(row + 16), hi));
578 }
579 }
580
582 inline void avx512_tile_c32(std::size_t kc,
583 const std::complex<float>* ap,
584 const std::complex<float>* bp,
585 std::size_t ldb,
586 std::complex<float>* c,
587 std::size_t ldc,
588 std::size_t rows,
589 std::size_t cols) noexcept
590 {
591 through_edge<avx512_complex_mr, avx512_complex_nr_f32>(c, ldc, rows, cols, [&](std::complex<float>* tile, std::size_t ld) noexcept {
592 avx512_full_tile_c32(kc, ap, bp, ldb, tile, ld);
593 });
594 }
595
601 NUMERICAL_TARGET("avx512f") inline double avx512_symmetric_row(const double* row,
602 const double* v,
603 double f,
604 double* out,
605 std::size_t count) noexcept
606 {
607 constexpr std::size_t width {8};
608 constexpr std::size_t lanes {4};
609 const __m512d splat {_mm512_set1_pd(f)};
610 std::array<zmm_register, lanes> sum {};
611 const std::size_t done {count - (count % (width * lanes))};
612 for (std::size_t i {0}; i < done; i += width * lanes) {
613 NUMERICAL_UNROLL_TILE
614 for (std::size_t l {0}; l < lanes; ++l) {
615 const __m512d r {_mm512_loadu_pd(row + i + (l * width))};
616 sum[l].v = _mm512_fmadd_pd(r, _mm512_loadu_pd(v + i + (l * width)), sum[l].v);
617 _mm512_storeu_pd(out + i + (l * width), _mm512_fmadd_pd(r, splat, _mm512_loadu_pd(out + i + (l * width))));
618 }
619 }
620 double total {0.0};
621 for (std::size_t l {0}; l < lanes; ++l) {
622 total += _mm512_reduce_add_pd(sum[l].v);
623 }
624 for (std::size_t i {done}; i < count; ++i) {
625 total += row[i] * v[i];
626 out[i] += row[i] * f;
627 }
628 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
629 return total;
630 }
631
633 NUMERICAL_TARGET("avx2,fma") inline double avx2_symmetric_row(const double* row,
634 const double* v,
635 double f,
636 double* out,
637 std::size_t count) noexcept
638 {
639 constexpr std::size_t width {4};
640 constexpr std::size_t lanes {4};
641 const __m256d splat {_mm256_set1_pd(f)};
642 std::array<ymm_register, lanes> sum {};
643 const std::size_t done {count - (count % (width * lanes))};
644 for (std::size_t i {0}; i < done; i += width * lanes) {
645 NUMERICAL_UNROLL_TILE
646 for (std::size_t l {0}; l < lanes; ++l) {
647 const __m256d r {_mm256_loadu_pd(row + i + (l * width))};
648 sum[l].v = _mm256_fmadd_pd(r, _mm256_loadu_pd(v + i + (l * width)), sum[l].v);
649 _mm256_storeu_pd(out + i + (l * width), _mm256_fmadd_pd(r, splat, _mm256_loadu_pd(out + i + (l * width))));
650 }
651 }
652 double total {0.0};
653 for (std::size_t l {0}; l < lanes; ++l) {
654 alignas(32) std::array<double, width> part {}; // a register stored whole: some compilers store it aligned
655 _mm256_storeu_pd(part.data(), sum[l].v);
656 total += (part[0] + part[1]) + (part[2] + part[3]);
657 }
658 for (std::size_t i {done}; i < count; ++i) {
659 total += row[i] * v[i];
660 out[i] += row[i] * f;
661 }
662 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
663 return total;
664 }
665
670 NUMERICAL_TARGET("avx512f") inline std::array<double, symmetric_row_group> avx512_symmetric_rows(const double* rows,
671 std::size_t n,
672 const double* v,
673 double* out,
674 std::size_t i) noexcept
675 {
676 constexpr std::size_t width {8};
677 std::array<zmm_register, symmetric_row_group> sum {};
678 std::array<zmm_register, symmetric_row_group> splat {};
679 for (std::size_t r {0}; r < symmetric_row_group; ++r) {
680 splat[r].v = _mm512_set1_pd(v[i + r]);
681 }
682 const std::size_t done {i - (i % width)};
683 for (std::size_t j {0}; j < done; j += width) {
684 const __m512d vv {_mm512_loadu_pd(v + j)};
685 __m512d ov {_mm512_loadu_pd(out + j)};
686 NUMERICAL_UNROLL_TILE
687 for (std::size_t r {0}; r < symmetric_row_group; ++r) {
688 const __m512d rv {_mm512_loadu_pd(rows + (r * n) + j)};
689 sum[r].v = _mm512_fmadd_pd(rv, vv, sum[r].v);
690 ov = _mm512_fmadd_pd(rv, splat[r].v, ov);
691 }
692 _mm512_storeu_pd(out + j, ov);
693 }
694 std::array<double, symmetric_row_group> sums {};
695 for (std::size_t r {0}; r < symmetric_row_group; ++r) {
696 sums[r] = _mm512_reduce_add_pd(sum[r].v);
697 }
698 _mm256_zeroupper(); // the upper halves cleared: the tail may be SSE alone
699 return symmetric_rows_tail(rows, n, v, out, i, done, sums);
700 }
701
703 NUMERICAL_TARGET("avx2,fma") inline std::array<double, symmetric_row_group> avx2_symmetric_rows(const double* rows,
704 std::size_t n,
705 const double* v,
706 double* out,
707 std::size_t i) noexcept
708 {
709 constexpr std::size_t width {4};
710 std::array<ymm_register, symmetric_row_group> sum {};
711 std::array<ymm_register, symmetric_row_group> splat {};
712 for (std::size_t r {0}; r < symmetric_row_group; ++r) {
713 splat[r].v = _mm256_set1_pd(v[i + r]);
714 }
715 const std::size_t done {i - (i % width)};
716 for (std::size_t j {0}; j < done; j += width) {
717 const __m256d vv {_mm256_loadu_pd(v + j)};
718 __m256d ov {_mm256_loadu_pd(out + j)};
719 NUMERICAL_UNROLL_TILE
720 for (std::size_t r {0}; r < symmetric_row_group; ++r) {
721 const __m256d rv {_mm256_loadu_pd(rows + (r * n) + j)};
722 sum[r].v = _mm256_fmadd_pd(rv, vv, sum[r].v);
723 ov = _mm256_fmadd_pd(rv, splat[r].v, ov);
724 }
725 _mm256_storeu_pd(out + j, ov);
726 }
727 std::array<double, symmetric_row_group> sums {};
728 for (std::size_t r {0}; r < symmetric_row_group; ++r) {
729 alignas(32) std::array<double, width> part {}; // a register stored whole: some compilers store it aligned
730 _mm256_storeu_pd(part.data(), sum[r].v);
731 sums[r] = (part[0] + part[1]) + (part[2] + part[3]);
732 }
733 _mm256_zeroupper(); // the upper halves cleared: the tail may be SSE alone
734 return symmetric_rows_tail(rows, n, v, out, i, done, sums);
735 }
736
739 NUMERICAL_TARGET("avx512f") inline void avx512_scaled_add(double* y,
740 const double* x,
741 double f,
742 std::size_t count) noexcept
743 {
744 const __m512d splat {_mm512_set1_pd(f)};
745 const std::size_t done {count - (count % 8)};
746 for (std::size_t i {0}; i < done; i += 8) {
747 _mm512_storeu_pd(y + i, _mm512_fmadd_pd(splat, _mm512_loadu_pd(x + i), _mm512_loadu_pd(y + i)));
748 }
749 for (std::size_t i {done}; i < count; ++i) {
750 y[i] = std::fma(f, x[i], y[i]);
751 }
752 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
753 }
754
760 NUMERICAL_TARGET("avx512f") inline void avx512_subtract_products(double* y,
761 const double* f,
762 const double* x,
763 std::size_t ldx,
764 std::size_t terms,
765 std::size_t count) noexcept
766 {
767 std::size_t i {0};
768 for (; i + 32 <= count; i += 32) {
769 __m512d y0 {_mm512_loadu_pd(y + i)};
770 __m512d y1 {_mm512_loadu_pd(y + i + 8)};
771 __m512d y2 {_mm512_loadu_pd(y + i + 16)};
772 __m512d y3 {_mm512_loadu_pd(y + i + 24)};
773 for (std::size_t t {0}; t < terms; ++t) {
774 const __m512d minus {_mm512_set1_pd(-f[t])};
775 const double* const row {x + (t * ldx) + i};
776 y0 = _mm512_fmadd_pd(minus, _mm512_loadu_pd(row), y0);
777 y1 = _mm512_fmadd_pd(minus, _mm512_loadu_pd(row + 8), y1);
778 y2 = _mm512_fmadd_pd(minus, _mm512_loadu_pd(row + 16), y2);
779 y3 = _mm512_fmadd_pd(minus, _mm512_loadu_pd(row + 24), y3);
780 }
781 _mm512_storeu_pd(y + i, y0);
782 _mm512_storeu_pd(y + i + 8, y1);
783 _mm512_storeu_pd(y + i + 16, y2);
784 _mm512_storeu_pd(y + i + 24, y3);
785 }
786 for (; i + 8 <= count; i += 8) {
787 __m512d sum {_mm512_loadu_pd(y + i)};
788 for (std::size_t t {0}; t < terms; ++t) {
789 sum = _mm512_fmadd_pd(_mm512_set1_pd(-f[t]), _mm512_loadu_pd(x + (t * ldx) + i), sum);
790 }
791 _mm512_storeu_pd(y + i, sum);
792 }
793 for (; i < count; ++i) {
794 for (std::size_t t {0}; t < terms; ++t) {
795 y[i] = std::fma(-f[t], x[(t * ldx) + i], y[i]);
796 }
797 }
798 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
799 }
800
806 NUMERICAL_TARGET("avx512f") inline void avx512_subtract_rows(double* y,
807 std::size_t ldy,
808 const double* f,
809 std::size_t ldf,
810 const double* x,
811 std::size_t ldx,
812 std::size_t terms,
813 std::size_t count) noexcept
814 {
815 struct lane // a register in a struct: an array of the register type itself loses its alignment attribute
816 {
817 __m512d value;
818 };
819 for (std::size_t i {0}; i < count; i += 8) {
820 const auto mask {static_cast<__mmask8>(count - i >= 8 ? 0xFFU : (1U << (count - i)) - 1U)};
821 std::array<lane, 8> held {};
822 for (std::size_t r {0}; r < 8; ++r) {
823 held[r].value = _mm512_maskz_loadu_pd(mask, y + (r * ldy) + i);
824 }
825 for (std::size_t t {0}; t < terms; ++t) {
826 const __m512d row {_mm512_maskz_loadu_pd(mask, x + (t * ldx) + i)};
827 for (std::size_t r {0}; r < 8; ++r) {
828 held[r].value = _mm512_fmadd_pd(_mm512_set1_pd(-f[(r * ldf) + t]), row, held[r].value);
829 }
830 }
831 for (std::size_t r {0}; r < 8; ++r) {
832 _mm512_mask_storeu_pd(y + (r * ldy) + i, mask, held[r].value);
833 }
834 }
835 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
836 }
837
839 NUMERICAL_TARGET("avx2,fma") inline void avx2_subtract_rows(double* y,
840 std::size_t ldy,
841 const double* f,
842 std::size_t ldf,
843 const double* x,
844 std::size_t ldx,
845 std::size_t terms,
846 std::size_t count) noexcept
847 {
848 struct lane // as avx512_subtract_rows's
849 {
850 __m256d value;
851 };
852 std::size_t i {0};
853 for (; i + 4 <= count; i += 4) {
854 std::array<lane, 8> held {};
855 for (std::size_t r {0}; r < 8; ++r) {
856 held[r].value = _mm256_loadu_pd(y + (r * ldy) + i);
857 }
858 for (std::size_t t {0}; t < terms; ++t) {
859 const __m256d row {_mm256_loadu_pd(x + (t * ldx) + i)};
860 for (std::size_t r {0}; r < 8; ++r) {
861 held[r].value = _mm256_fmadd_pd(_mm256_set1_pd(-f[(r * ldf) + t]), row, held[r].value);
862 }
863 }
864 for (std::size_t r {0}; r < 8; ++r) {
865 _mm256_storeu_pd(y + (r * ldy) + i, held[r].value);
866 }
867 }
868 for (; i < count; ++i) {
869 std::array<double, 8> held {};
870 for (std::size_t r {0}; r < 8; ++r) {
871 held[r] = y[(r * ldy) + i];
872 }
873 for (std::size_t t {0}; t < terms; ++t) {
874 const double value {x[(t * ldx) + i]};
875 for (std::size_t r {0}; r < 8; ++r) {
876 held[r] = std::fma(-f[(r * ldf) + t], value, held[r]);
877 }
878 }
879 for (std::size_t r {0}; r < 8; ++r) {
880 y[(r * ldy) + i] = held[r];
881 }
882 }
883 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
884 }
885
887 NUMERICAL_TARGET("avx2,fma") inline void avx2_subtract_products(double* y,
888 const double* f,
889 const double* x,
890 std::size_t ldx,
891 std::size_t terms,
892 std::size_t count) noexcept
893 {
894 std::size_t i {0};
895 for (; i + 16 <= count; i += 16) {
896 __m256d y0 {_mm256_loadu_pd(y + i)};
897 __m256d y1 {_mm256_loadu_pd(y + i + 4)};
898 __m256d y2 {_mm256_loadu_pd(y + i + 8)};
899 __m256d y3 {_mm256_loadu_pd(y + i + 12)};
900 for (std::size_t t {0}; t < terms; ++t) {
901 const __m256d minus {_mm256_set1_pd(-f[t])};
902 const double* const row {x + (t * ldx) + i};
903 y0 = _mm256_fmadd_pd(minus, _mm256_loadu_pd(row), y0);
904 y1 = _mm256_fmadd_pd(minus, _mm256_loadu_pd(row + 4), y1);
905 y2 = _mm256_fmadd_pd(minus, _mm256_loadu_pd(row + 8), y2);
906 y3 = _mm256_fmadd_pd(minus, _mm256_loadu_pd(row + 12), y3);
907 }
908 _mm256_storeu_pd(y + i, y0);
909 _mm256_storeu_pd(y + i + 4, y1);
910 _mm256_storeu_pd(y + i + 8, y2);
911 _mm256_storeu_pd(y + i + 12, y3);
912 }
913 for (; i + 4 <= count; i += 4) {
914 __m256d sum {_mm256_loadu_pd(y + i)};
915 for (std::size_t t {0}; t < terms; ++t) {
916 sum = _mm256_fmadd_pd(_mm256_set1_pd(-f[t]), _mm256_loadu_pd(x + (t * ldx) + i), sum);
917 }
918 _mm256_storeu_pd(y + i, sum);
919 }
920 for (; i < count; ++i) {
921 for (std::size_t t {0}; t < terms; ++t) {
922 y[i] = std::fma(-f[t], x[(t * ldx) + i], y[i]);
923 }
924 }
925 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
926 }
927
934 NUMERICAL_TARGET("avx512f") inline void avx512_subtract_complex_products(double* y,
935 const double* f,
936 const double* x,
937 std::size_t ldx,
938 std::size_t terms,
939 std::size_t count) noexcept
940 {
941 const __m512d sign {_mm512_set_pd(1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0)}; // (−, +) for each number
942 std::size_t i {0};
943 for (; i + 16 <= count; i += 16) {
944 std::array<zmm_register, 4> acc {};
945 NUMERICAL_UNROLL_TILE
946 for (std::size_t r {0}; r < 4; ++r) {
947 acc[r].v = _mm512_loadu_pd(y + (2 * i) + (8 * r));
948 }
949 for (std::size_t t {0}; t < terms; ++t) {
950 const __m512d real {_mm512_set1_pd(f[2 * t])};
951 const __m512d imag {_mm512_mul_pd(_mm512_set1_pd(f[(2 * t) + 1]), sign)};
952 const double* const row {x + (2 * ((t * ldx) + i))};
953 NUMERICAL_UNROLL_TILE
954 for (std::size_t r {0}; r < 4; ++r) {
955 const __m512d value {_mm512_loadu_pd(row + (8 * r))};
956 acc[r].v = _mm512_fnmadd_pd(real, value, acc[r].v);
957 acc[r].v = _mm512_fnmadd_pd(imag, _mm512_permute_pd(value, 0x55), acc[r].v);
958 }
959 }
960 NUMERICAL_UNROLL_TILE
961 for (std::size_t r {0}; r < 4; ++r) {
962 _mm512_storeu_pd(y + (2 * i) + (8 * r), acc[r].v);
963 }
964 }
965 for (; i + 4 <= count; i += 4) {
966 __m512d sum {_mm512_loadu_pd(y + (2 * i))};
967 for (std::size_t t {0}; t < terms; ++t) {
968 const __m512d value {_mm512_loadu_pd(x + (2 * ((t * ldx) + i)))};
969 sum = _mm512_fnmadd_pd(_mm512_set1_pd(f[2 * t]), value, sum);
970 sum = _mm512_fnmadd_pd(_mm512_mul_pd(_mm512_set1_pd(f[(2 * t) + 1]), sign), _mm512_permute_pd(value, 0x55), sum);
971 }
972 _mm512_storeu_pd(y + (2 * i), sum);
973 }
974 for (; i < count; ++i) {
975 double real {y[2 * i]};
976 double imag {y[(2 * i) + 1]};
977 for (std::size_t t {0}; t < terms; ++t) {
978 const double xr {x[2 * ((t * ldx) + i)]};
979 const double xi {x[(2 * ((t * ldx) + i)) + 1]};
980 real = std::fma(-f[2 * t], xr, real);
981 imag = std::fma(-f[2 * t], xi, imag);
982 real = std::fma(f[(2 * t) + 1], xi, real);
983 imag = std::fma(-f[(2 * t) + 1], xr, imag);
984 }
985 y[2 * i] = real;
986 y[(2 * i) + 1] = imag;
987 }
988 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
989 }
990
992 NUMERICAL_TARGET("avx2,fma") inline void avx2_subtract_complex_products(double* y,
993 const double* f,
994 const double* x,
995 std::size_t ldx,
996 std::size_t terms,
997 std::size_t count) noexcept
998 {
999 const __m256d sign {_mm256_set_pd(1.0, -1.0, 1.0, -1.0)};
1000 std::size_t i {0};
1001 for (; i + 8 <= count; i += 8) {
1002 std::array<ymm_register, 4> acc {};
1003 NUMERICAL_UNROLL_TILE
1004 for (std::size_t r {0}; r < 4; ++r) {
1005 acc[r].v = _mm256_loadu_pd(y + (2 * i) + (4 * r));
1006 }
1007 for (std::size_t t {0}; t < terms; ++t) {
1008 const __m256d real {_mm256_set1_pd(f[2 * t])};
1009 const __m256d imag {_mm256_mul_pd(_mm256_set1_pd(f[(2 * t) + 1]), sign)};
1010 const double* const row {x + (2 * ((t * ldx) + i))};
1011 NUMERICAL_UNROLL_TILE
1012 for (std::size_t r {0}; r < 4; ++r) {
1013 const __m256d value {_mm256_loadu_pd(row + (4 * r))};
1014 acc[r].v = _mm256_fnmadd_pd(real, value, acc[r].v);
1015 acc[r].v = _mm256_fnmadd_pd(imag, _mm256_permute_pd(value, 0x5), acc[r].v);
1016 }
1017 }
1018 NUMERICAL_UNROLL_TILE
1019 for (std::size_t r {0}; r < 4; ++r) {
1020 _mm256_storeu_pd(y + (2 * i) + (4 * r), acc[r].v);
1021 }
1022 }
1023 for (; i + 2 <= count; i += 2) {
1024 __m256d sum {_mm256_loadu_pd(y + (2 * i))};
1025 for (std::size_t t {0}; t < terms; ++t) {
1026 const __m256d value {_mm256_loadu_pd(x + (2 * ((t * ldx) + i)))};
1027 sum = _mm256_fnmadd_pd(_mm256_set1_pd(f[2 * t]), value, sum);
1028 sum = _mm256_fnmadd_pd(_mm256_mul_pd(_mm256_set1_pd(f[(2 * t) + 1]), sign), _mm256_permute_pd(value, 0x5), sum);
1029 }
1030 _mm256_storeu_pd(y + (2 * i), sum);
1031 }
1032 for (; i < count; ++i) {
1033 double real {y[2 * i]};
1034 double imag {y[(2 * i) + 1]};
1035 for (std::size_t t {0}; t < terms; ++t) {
1036 const double xr {x[2 * ((t * ldx) + i)]};
1037 const double xi {x[(2 * ((t * ldx) + i)) + 1]};
1038 real = std::fma(-f[2 * t], xr, real);
1039 imag = std::fma(-f[2 * t], xi, imag);
1040 real = std::fma(f[(2 * t) + 1], xi, real);
1041 imag = std::fma(-f[(2 * t) + 1], xr, imag);
1042 }
1043 y[2 * i] = real;
1044 y[(2 * i) + 1] = imag;
1045 }
1046 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
1047 }
1048
1055 NUMERICAL_TARGET("avx512f") inline std::complex<double> avx512_hermitian_row(const double* row,
1056 const double* v,
1057 double vr,
1058 double vi,
1059 double* out,
1060 std::size_t count) noexcept
1061 {
1062 const __m512d real {_mm512_set_pd(-vr, vr, -vr, vr, -vr, vr, -vr, vr)}; // (vr, −vr) for each number
1063 const __m512d imag {_mm512_set1_pd(vi)};
1064 __m512d same {_mm512_setzero_pd()}; // rr·vr, ri·vi
1065 __m512d across {_mm512_setzero_pd()}; // rr·vi, ri·vr
1066 std::size_t j {0};
1067 for (; j + 4 <= count; j += 4) {
1068 const __m512d r {_mm512_loadu_pd(row + (2 * j))};
1069 const __m512d x {_mm512_loadu_pd(v + (2 * j))};
1070 const __m512d s {_mm512_permute_pd(r, 0x55)};
1071 __m512d o {_mm512_loadu_pd(out + (2 * j))};
1072 o = _mm512_fmadd_pd(real, r, o);
1073 o = _mm512_fmadd_pd(imag, s, o);
1074 same = _mm512_fmadd_pd(r, x, same);
1075 across = _mm512_fmadd_pd(s, x, across);
1076 _mm512_storeu_pd(out + (2 * j), o);
1077 }
1078 // Registers stored whole: GCC 14 stores the second aligned (in C++20, under the address sanitizer), so both are.
1079 alignas(64) std::array<double, 8> lanes_same {};
1080 alignas(64) std::array<double, 8> lanes_across {};
1081 _mm512_storeu_pd(lanes_same.data(), same);
1082 _mm512_storeu_pd(lanes_across.data(), across);
1083 double re {0.0};
1084 double im {0.0};
1085 for (std::size_t l {0}; l < 8; l += 2) {
1086 re += lanes_same[l] - lanes_same[l + 1];
1087 im += lanes_across[l] + lanes_across[l + 1];
1088 }
1089 for (; j < count; ++j) {
1090 const double rr {row[2 * j]};
1091 const double ri {row[(2 * j) + 1]};
1092 out[2 * j] = std::fma(vi, ri, std::fma(vr, rr, out[2 * j]));
1093 out[(2 * j) + 1] = std::fma(vi, rr, std::fma(-vr, ri, out[(2 * j) + 1]));
1094 re += (rr * v[2 * j]) - (ri * v[(2 * j) + 1]);
1095 im += (ri * v[2 * j]) + (rr * v[(2 * j) + 1]);
1096 }
1097 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
1098 return std::complex<double> {re, im};
1099 }
1100
1102 NUMERICAL_TARGET("avx2,fma") inline std::complex<double> avx2_hermitian_row(const double* row,
1103 const double* v,
1104 double vr,
1105 double vi,
1106 double* out,
1107 std::size_t count) noexcept
1108 {
1109 const __m256d real {_mm256_set_pd(-vr, vr, -vr, vr)};
1110 const __m256d imag {_mm256_set1_pd(vi)};
1111 __m256d same {_mm256_setzero_pd()};
1112 __m256d across {_mm256_setzero_pd()};
1113 std::size_t j {0};
1114 for (; j + 2 <= count; j += 2) {
1115 const __m256d r {_mm256_loadu_pd(row + (2 * j))};
1116 const __m256d x {_mm256_loadu_pd(v + (2 * j))};
1117 const __m256d s {_mm256_permute_pd(r, 0x5)};
1118 __m256d o {_mm256_loadu_pd(out + (2 * j))};
1119 o = _mm256_fmadd_pd(real, r, o);
1120 o = _mm256_fmadd_pd(imag, s, o);
1121 same = _mm256_fmadd_pd(r, x, same);
1122 across = _mm256_fmadd_pd(s, x, across);
1123 _mm256_storeu_pd(out + (2 * j), o);
1124 }
1125 alignas(32) std::array<double, 4> lanes_same {}; // registers stored whole, aligned (as above)
1126 alignas(32) std::array<double, 4> lanes_across {};
1127 _mm256_storeu_pd(lanes_same.data(), same);
1128 _mm256_storeu_pd(lanes_across.data(), across);
1129 double re {(lanes_same[0] - lanes_same[1]) + (lanes_same[2] - lanes_same[3])};
1130 double im {(lanes_across[0] + lanes_across[1]) + (lanes_across[2] + lanes_across[3])};
1131 for (; j < count; ++j) {
1132 const double rr {row[2 * j]};
1133 const double ri {row[(2 * j) + 1]};
1134 out[2 * j] = std::fma(vi, ri, std::fma(vr, rr, out[2 * j]));
1135 out[(2 * j) + 1] = std::fma(vi, rr, std::fma(-vr, ri, out[(2 * j) + 1]));
1136 re += (rr * v[2 * j]) - (ri * v[(2 * j) + 1]);
1137 im += (ri * v[2 * j]) + (rr * v[(2 * j) + 1]);
1138 }
1139 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
1140 return std::complex<double> {re, im};
1141 }
1142
1144 NUMERICAL_TARGET("avx2,fma") inline void avx2_scaled_add(double* y,
1145 const double* x,
1146 double f,
1147 std::size_t count) noexcept
1148 {
1149 const __m256d splat {_mm256_set1_pd(f)};
1150 const std::size_t done {count - (count % 4)};
1151 for (std::size_t i {0}; i < done; i += 4) {
1152 _mm256_storeu_pd(y + i, _mm256_fmadd_pd(splat, _mm256_loadu_pd(x + i), _mm256_loadu_pd(y + i)));
1153 }
1154 for (std::size_t i {done}; i < count; ++i) {
1155 y[i] = std::fma(f, x[i], y[i]);
1156 }
1157 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
1158 }
1159
1166 NUMERICAL_TARGET("avx2") inline std::size_t avx2_scaled_add_c32(std::complex<float>* y,
1167 const std::complex<float>* x,
1168 std::complex<float> f,
1169 std::size_t count) noexcept
1170 {
1171 auto* const out {reinterpret_cast<float*>(y)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee of [complex.numbers]
1172 const auto* const in {reinterpret_cast<const float*>(x)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the same
1173 const __m256 re {_mm256_set1_ps(f.real())};
1174 const __m256 im {_mm256_set1_ps(f.imag())};
1175 const std::size_t done {count - (count % 4)};
1176 for (std::size_t i {0}; i < 2 * done; i += 8) {
1177 const __m256 xv {_mm256_loadu_ps(in + i)};
1178 const __m256 product {_mm256_addsub_ps(_mm256_mul_ps(re, xv), _mm256_mul_ps(im, _mm256_permute_ps(xv, 0xB1)))};
1179 _mm256_storeu_ps(out + i, _mm256_add_ps(_mm256_loadu_ps(out + i), product));
1180 }
1181 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
1182 return done;
1183 }
1184
1186 NUMERICAL_TARGET("avx2") inline std::size_t avx2_scaled_add_c64(std::complex<double>* y,
1187 const std::complex<double>* x,
1188 std::complex<double> f,
1189 std::size_t count) noexcept
1190 {
1191 auto* const out {reinterpret_cast<double*>(y)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the layout guarantee of [complex.numbers]
1192 const auto* const in {reinterpret_cast<const double*>(x)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): the same
1193 const __m256d re {_mm256_set1_pd(f.real())};
1194 const __m256d im {_mm256_set1_pd(f.imag())};
1195 const std::size_t done {count - (count % 2)};
1196 for (std::size_t i {0}; i < 2 * done; i += 4) {
1197 const __m256d xv {_mm256_loadu_pd(in + i)};
1198 const __m256d product {_mm256_addsub_pd(_mm256_mul_pd(re, xv), _mm256_mul_pd(im, _mm256_permute_pd(xv, 0b0101)))};
1199 _mm256_storeu_pd(out + i, _mm256_add_pd(_mm256_loadu_pd(out + i), product));
1200 }
1201 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
1202 return done;
1203 }
1204
1211 NUMERICAL_TARGET("avx2,fma") inline void avx2_transposed_columns(const double* a,
1212 const double* bt,
1213 std::size_t ldbt,
1214 std::size_t inner,
1215 std::size_t first,
1216 std::size_t last,
1217 double* y) noexcept
1218 {
1219 constexpr std::size_t group {8};
1220 std::size_t j {first};
1221 for (; j + group <= last; j += group) {
1222 std::array<double, group> sums {};
1223 for (std::size_t k {0}; k < inner; ++k) {
1224 const double f {a[k]};
1225 for (std::size_t c {0}; c < group; ++c) {
1226 sums[c] = std::fma(f, bt[((j + c) * ldbt) + k], sums[c]);
1227 }
1228 }
1229 std::copy(sums.begin(), sums.end(), y + j);
1230 }
1231 for (; j < last; ++j) {
1232 double sum {};
1233 for (std::size_t k {0}; k < inner; ++k) {
1234 sum = std::fma(a[k], bt[(j * ldbt) + k], sum);
1235 }
1236 y[j] = sum;
1237 }
1238 }
1239
1242 NUMERICAL_TARGET("avx512f") inline double avx512_dot(const double* x,
1243 const double* y,
1244 std::size_t count) noexcept
1245 {
1246 constexpr std::size_t width {8};
1247 constexpr std::size_t lanes {4};
1248 std::array<zmm_register, lanes> sum {};
1249 const std::size_t done {count - (count % (width * lanes))};
1250 for (std::size_t i {0}; i < done; i += width * lanes) {
1251 NUMERICAL_UNROLL_TILE
1252 for (std::size_t l {0}; l < lanes; ++l) {
1253 sum[l].v = _mm512_fmadd_pd(_mm512_loadu_pd(x + i + (l * width)), _mm512_loadu_pd(y + i + (l * width)), sum[l].v);
1254 }
1255 }
1256 double total {0.0};
1257 for (std::size_t l {0}; l < lanes; ++l) {
1258 total += _mm512_reduce_add_pd(sum[l].v);
1259 }
1260 for (std::size_t i {done}; i < count; ++i) {
1261 total += x[i] * y[i];
1262 }
1263 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
1264 return total;
1265 }
1266
1268 NUMERICAL_TARGET("avx2,fma") inline double avx2_dot(const double* x,
1269 const double* y,
1270 std::size_t count) noexcept
1271 {
1272 constexpr std::size_t width {4};
1273 constexpr std::size_t lanes {4};
1274 std::array<ymm_register, lanes> sum {};
1275 const std::size_t done {count - (count % (width * lanes))};
1276 for (std::size_t i {0}; i < done; i += width * lanes) {
1277 NUMERICAL_UNROLL_TILE
1278 for (std::size_t l {0}; l < lanes; ++l) {
1279 sum[l].v = _mm256_fmadd_pd(_mm256_loadu_pd(x + i + (l * width)), _mm256_loadu_pd(y + i + (l * width)), sum[l].v);
1280 }
1281 }
1282 double total {0.0};
1283 for (std::size_t l {0}; l < lanes; ++l) {
1284 alignas(32) std::array<double, width> part {}; // a register stored whole: some compilers store it aligned
1285 _mm256_storeu_pd(part.data(), sum[l].v);
1286 total += (part[0] + part[1]) + (part[2] + part[3]);
1287 }
1288 for (std::size_t i {done}; i < count; ++i) {
1289 total += x[i] * y[i];
1290 }
1291 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
1292 return total;
1293 }
1294
1302 NUMERICAL_TARGET("avx2,fma") inline std::complex<double> avx2_complex_dot(const std::complex<double>* x,
1303 const std::complex<double>* y,
1304 std::size_t count) noexcept
1305 {
1306 constexpr std::size_t width {4}; // parts in a register: two complex numbers
1307 constexpr std::size_t lanes {4};
1308 const auto* const a {reinterpret_cast<const double*>(x)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): a complex number is an array of its two parts
1309 const auto* const b {reinterpret_cast<const double*>(y)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): as above
1310 const std::size_t parts {2 * count};
1311 std::array<ymm_register, lanes> same {};
1312 std::array<ymm_register, lanes> cross {};
1313 const std::size_t done {parts - (parts % (width * lanes))};
1314 for (std::size_t i {0}; i < done; i += width * lanes) {
1315 NUMERICAL_UNROLL_TILE
1316 for (std::size_t l {0}; l < lanes; ++l) {
1317 const __m256d av {_mm256_loadu_pd(a + i + (l * width))};
1318 const __m256d bv {_mm256_loadu_pd(b + i + (l * width))};
1319 same[l].v = _mm256_fmadd_pd(av, bv, same[l].v);
1320 cross[l].v = _mm256_fmadd_pd(av, _mm256_permute_pd(bv, 0x5), cross[l].v); // each pair's parts swapped
1321 }
1322 }
1323 double rr {0.0};
1324 double ii {0.0};
1325 double ri {0.0};
1326 double ir {0.0};
1327 for (std::size_t l {0}; l < lanes; ++l) {
1328 alignas(32) std::array<double, width> s {}; // a register stored whole: some compilers store it aligned
1329 alignas(32) std::array<double, width> c {};
1330 _mm256_storeu_pd(s.data(), same[l].v);
1331 _mm256_storeu_pd(c.data(), cross[l].v);
1332 rr += s[0] + s[2];
1333 ii += s[1] + s[3];
1334 ri += c[0] + c[2];
1335 ir += c[1] + c[3];
1336 }
1337 for (std::size_t i {done}; i < parts; i += 2) {
1338 rr += a[i] * b[i];
1339 ii += a[i + 1] * b[i + 1];
1340 ri += a[i] * b[i + 1];
1341 ir += a[i + 1] * b[i];
1342 }
1343 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
1344 return {rr - ii, ri + ir};
1345 }
1346
1348 NUMERICAL_TARGET("avx2,fma") inline std::complex<float> avx2_complex_dot(const std::complex<float>* x,
1349 const std::complex<float>* y,
1350 std::size_t count) noexcept
1351 {
1352 constexpr std::size_t width {8};
1353 constexpr std::size_t lanes {4};
1354 const auto* const a {reinterpret_cast<const float*>(x)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): a complex number is an array of its two parts
1355 const auto* const b {reinterpret_cast<const float*>(y)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): as above
1356 const std::size_t parts {2 * count};
1357 std::array<ymm_register_f32, lanes> same {};
1358 std::array<ymm_register_f32, lanes> cross {};
1359 const std::size_t done {parts - (parts % (width * lanes))};
1360 for (std::size_t i {0}; i < done; i += width * lanes) {
1361 NUMERICAL_UNROLL_TILE
1362 for (std::size_t l {0}; l < lanes; ++l) {
1363 const __m256 av {_mm256_loadu_ps(a + i + (l * width))};
1364 const __m256 bv {_mm256_loadu_ps(b + i + (l * width))};
1365 same[l].v = _mm256_fmadd_ps(av, bv, same[l].v);
1366 cross[l].v = _mm256_fmadd_ps(av, _mm256_permute_ps(bv, 0xB1), cross[l].v); // each pair's parts swapped
1367 }
1368 }
1369 float rr {0.0F};
1370 float ii {0.0F};
1371 float ri {0.0F};
1372 float ir {0.0F};
1373 for (std::size_t l {0}; l < lanes; ++l) {
1374 alignas(32) std::array<float, width> s {};
1375 alignas(32) std::array<float, width> c {};
1376 _mm256_storeu_ps(s.data(), same[l].v);
1377 _mm256_storeu_ps(c.data(), cross[l].v);
1378 rr += (s[0] + s[2]) + (s[4] + s[6]);
1379 ii += (s[1] + s[3]) + (s[5] + s[7]);
1380 ri += (c[0] + c[2]) + (c[4] + c[6]);
1381 ir += (c[1] + c[3]) + (c[5] + c[7]);
1382 }
1383 for (std::size_t i {done}; i < parts; i += 2) {
1384 rr += a[i] * b[i];
1385 ii += a[i + 1] * b[i + 1];
1386 ri += a[i] * b[i + 1];
1387 ir += a[i + 1] * b[i];
1388 }
1389 _mm256_zeroupper(); // the upper halves cleared: code past the kernel may be SSE alone
1390 return {rr - ii, ri + ir};
1391 }
1392
1393#endif
1394} // namespace numerical::detail
1395
1396#endif // NUMERICAL_X86_HPP
numerical::detail::symmetric_row_group
constexpr std::size_t symmetric_row_group
The rows of a lower triangle symmetric_rows takes at once.
Definition simd.hpp:34
numerical::detail::symmetric_rows_tail
std::array< T, symmetric_row_group > symmetric_rows_tail(const T *rows, std::size_t n, const T *v, T *out, std::size_t i, std::size_t done, std::array< T, symmetric_row_group > sums)
The columns done…i+r of symmetric_rows, one by one: for each row r in order, out[k] += row_r[k] · v[i...
Definition simd.hpp:41
numerical::sum
T sum(const tensor< T > &t)
The sum of every element, by pairwise halving (0 for an empty tensor).
Definition tensor_ops.hpp:759
simd.hpp
The native vector of float and double the kernels hold their lanes in (GNU vector extensions,...
Generated by doxygen 1.9.8