NumericAL
A header-only C++20 lexer built on REAL
Loading...
Searching...
No Matches
simd.hpp
Go to the documentation of this file.
1
6#ifndef NUMERICAL_SIMD_HPP
7#define NUMERICAL_SIMD_HPP
8
9#include <algorithm>
10#include <array>
11#include <bit>
12#include <cmath>
13#include <complex>
14#include <cstddef>
15#include <cstring>
16#include <type_traits>
17
18// A function inlined wherever it is called, past the compiler's budget: a call there costs more than its body.
19#if defined(_MSC_VER) && !defined(__clang__)
20# define NUMERICAL_ALWAYS_INLINE __forceinline
21#else
22# define NUMERICAL_ALWAYS_INLINE [[gnu::always_inline]] inline // NOLINT(cppcoreguidelines-macro-usage): an attribute, which no function can stand for
23#endif
24
25namespace numerical::detail {
26
31 inline constexpr std::size_t transpose_tile {32};
32
34 inline constexpr std::size_t symmetric_row_group {4};
35
40 template <typename T>
41 [[nodiscard]] std::array<T, symmetric_row_group> symmetric_rows_tail(const T * rows,
42 std::size_t n,
43 const T * v,
44 T * out,
45 std::size_t i,
46 std::size_t done,
47 std::array<T, symmetric_row_group> sums)
48 {
49 for (std::size_t r {0}; r < symmetric_row_group; ++r) {
50 const T* row {rows + (r * n)};
51 const T factor {v[i + r]};
52 for (std::size_t k {done}; k < i + r; ++k) {
53 sums[r] += row[k] * v[k];
54 out[k] += row[k] * factor;
55 }
56 sums[r] += row[i + r] * factor;
57 }
58 return sums;
59 }
60
63 template <typename T>
64 concept split_complex = std::is_same_v<T, std::complex<float>> || std::is_same_v<T, std::complex<double>>;
65
66#if defined(__GNUC__)
68# if defined(__AVX__)
69 inline constexpr std::size_t native_vector_bytes {32};
70# else
71 inline constexpr std::size_t native_vector_bytes {16};
72# endif
73
75 template <typename T>
76 struct native_vector;
77
79 template <>
80 struct native_vector<float>
81 {
82 using type = float __attribute__((vector_size(native_vector_bytes)));
83 };
84
86 template <>
87 struct native_vector<double>
88 {
89 using type = double __attribute__((vector_size(native_vector_bytes)));
90 };
91
93 template <typename T>
94 concept native_scalar = std::is_same_v<T, float> || std::is_same_v<T, double>;
95
101 template <native_scalar T>
102 [[nodiscard]] T dot_product(const T* x,
103 const T* y,
104 std::size_t count)
105 {
106 using lane = typename native_vector<T>::type;
107 constexpr std::size_t width {native_vector_bytes / sizeof(T)};
108 constexpr std::size_t lanes {4};
109 std::array<lane, lanes> sum {};
110 const std::size_t done {count - (count % (width * lanes))};
111 for (std::size_t i {0}; i < done; i += width * lanes) {
112# pragma GCC unroll 4 // the sums in registers: GCC keeps a rolled array of vectors in memory
113 for (std::size_t l {0}; l < lanes; ++l) {
114 lane xv;
115 lane yv;
116 std::memcpy(&xv, x + i + (l * width), sizeof(lane));
117 std::memcpy(&yv, y + i + (l * width), sizeof(lane));
118 sum[l] += xv * yv;
119 }
120 }
121 T total {};
122 for (std::size_t l {0}; l < lanes; ++l) {
123 for (std::size_t k {0}; k < width; ++k) {
124 total += sum[l][k];
125 }
126 }
127 for (std::size_t i {done}; i < count; ++i) {
128 total += x[i] * y[i];
129 }
130 return total;
131 }
132
142 template <split_complex C>
143 [[nodiscard]] std::array<typename C::value_type, 4> complex_dot_sums(const C* x,
144 const C* y,
145 std::size_t count)
146 {
147 using R = typename C::value_type;
148 using lane = typename native_vector<R>::type;
149 constexpr std::size_t width {native_vector_bytes / sizeof(R)}; // parts in a vector: an even count
150 constexpr std::size_t lanes {4};
151 const auto* const a {reinterpret_cast<const R*>(x)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): a complex number is an array of its two parts
152 const auto* const b {reinterpret_cast<const R*>(y)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): as above
153 const std::size_t parts {2 * count};
154 std::array<lane, lanes> same {};
155 std::array<lane, lanes> cross {};
156 const std::size_t done {parts - (parts % (width * lanes))};
157 for (std::size_t i {0}; i < done; i += width * lanes) {
158# pragma GCC unroll 4 // the sums in registers
159 for (std::size_t l {0}; l < lanes; ++l) {
160 lane av;
161 lane bv;
162 std::memcpy(&av, a + i + (l * width), sizeof(lane));
163 std::memcpy(&bv, b + i + (l * width), sizeof(lane));
164 lane swapped;
165 for (std::size_t k {0}; k < width; ++k) {
166 swapped[k] = bv[k ^ 1U];
167 }
168 same[l] += av * bv;
169 cross[l] += av * swapped;
170 }
171 }
172 R rr {};
173 R ii {};
174 R ri {};
175 R ir {};
176 for (std::size_t l {0}; l < lanes; ++l) {
177 for (std::size_t k {0}; k < width; k += 2) {
178 rr += same[l][k];
179 ii += same[l][k + 1];
180 ri += cross[l][k];
181 ir += cross[l][k + 1];
182 }
183 }
184 for (std::size_t i {done}; i < parts; i += 2) {
185 rr += a[i] * b[i];
186 ii += a[i + 1] * b[i + 1];
187 ri += a[i] * b[i + 1];
188 ir += a[i + 1] * b[i];
189 }
190 return {rr, ii, ri, ir};
191 }
192
194 template <split_complex C>
195 [[nodiscard]] C dot_product(const C* x,
196 const C* y,
197 std::size_t count)
198 {
199 const auto [rr, ii, ri, ir] {complex_dot_sums(x, y, count)};
200 return C {rr - ii, ri + ir};
201 }
202
205 template <split_complex C>
206 [[nodiscard]] C conjugate_dot_product(const C* x,
207 const C* y,
208 std::size_t count)
209 {
210 const auto [rr, ii, ri, ir] {complex_dot_sums(x, y, count)};
211 return C {rr + ii, ri - ir};
212 }
213
221 template <native_scalar T>
222 void row_dots(const T* m,
223 std::size_t ld,
224 std::size_t rows,
225 const T* x,
226 std::size_t count,
227 T* out)
228 {
229 using lane = typename native_vector<T>::type;
230 constexpr std::size_t width {native_vector_bytes / sizeof(T)};
231 constexpr std::size_t lanes {4};
232 constexpr std::size_t group {4};
233 const std::size_t done {count - (count % (width * lanes))};
234 std::size_t r0 {0};
235 for (; r0 + group <= rows; r0 += group) {
236 std::array<std::array<lane, lanes>, group> sum {};
237 for (std::size_t i {0}; i < done; i += width * lanes) {
238# pragma GCC unroll 4 // the sums in registers: GCC keeps a rolled array of vectors in memory
239 for (std::size_t l {0}; l < lanes; ++l) {
240 lane xv;
241 std::memcpy(&xv, x + i + (l * width), sizeof(lane));
242# pragma GCC unroll 4
243 for (std::size_t r {0}; r < group; ++r) {
244 lane mv;
245 std::memcpy(&mv, m + ((r0 + r) * ld) + i + (l * width), sizeof(lane));
246 sum[r][l] += mv * xv;
247 }
248 }
249 }
250 for (std::size_t r {0}; r < group; ++r) {
251 T total {};
252 for (std::size_t l {0}; l < lanes; ++l) {
253 for (std::size_t k {0}; k < width; ++k) {
254 total += sum[r][l][k];
255 }
256 }
257 const T* const row {m + ((r0 + r) * ld)};
258 for (std::size_t i {done}; i < count; ++i) {
259 total += row[i] * x[i];
260 }
261 out[r0 + r] = total;
262 }
263 }
264 for (; r0 < rows; ++r0) {
265 out[r0] = dot_product(m + (r0 * ld), x, count);
266 }
267 }
268
274 template <native_scalar T>
275 [[nodiscard]] T symmetric_row(const T* row,
276 const T* v,
277 const T& f,
278 T* out,
279 std::size_t count)
280 {
281 using lane = typename native_vector<T>::type;
282 constexpr std::size_t width {native_vector_bytes / sizeof(T)};
283 constexpr std::size_t lanes {4};
284 const T factor {f}; // read once: it may be an element of v
285 const lane splat = factor - lane {}; // factor in every lane
286 std::array<lane, lanes> sum {};
287 const std::size_t done {count - (count % (width * lanes))};
288 for (std::size_t i {0}; i < done; i += width * lanes) {
289# pragma GCC unroll 4 // the sums in registers: GCC keeps a rolled array of vectors in memory
290 for (std::size_t l {0}; l < lanes; ++l) {
291 lane rv;
292 lane vv;
293 lane ov;
294 std::memcpy(&rv, row + i + (l * width), sizeof(lane));
295 std::memcpy(&vv, v + i + (l * width), sizeof(lane));
296 std::memcpy(&ov, out + i + (l * width), sizeof(lane));
297 sum[l] += rv * vv;
298 ov += rv * splat;
299 std::memcpy(out + i + (l * width), &ov, sizeof(lane));
300 }
301 }
302 T total {};
303 for (std::size_t l {0}; l < lanes; ++l) {
304 for (std::size_t k {0}; k < width; ++k) {
305 total += sum[l][k];
306 }
307 }
308 for (std::size_t i {done}; i < count; ++i) {
309 total += row[i] * v[i];
310 out[i] += row[i] * factor;
311 }
312 return total;
313 }
314
321 template <native_scalar T>
322 [[nodiscard]] std::array<T, symmetric_row_group> symmetric_rows(const T* rows,
323 std::size_t n,
324 const T* v,
325 T* out,
326 std::size_t i)
327 {
328 using lane = typename native_vector<T>::type;
329 constexpr std::size_t width {native_vector_bytes / sizeof(T)};
330 std::array<lane, symmetric_row_group> sum {};
331 std::array<lane, symmetric_row_group> splat {};
332 for (std::size_t r {0}; r < symmetric_row_group; ++r) {
333 splat[r] = v[i + r] - lane {}; // v[i + r] in every lane
334 }
335 const std::size_t done {i - (i % width)};
336 for (std::size_t j {0}; j < done; j += width) {
337 lane vv;
338 lane ov;
339 std::memcpy(&vv, v + j, sizeof(lane));
340 std::memcpy(&ov, out + j, sizeof(lane));
341# pragma GCC unroll 4 // the sums in registers: GCC keeps a rolled array of vectors in memory
342 for (std::size_t r {0}; r < symmetric_row_group; ++r) {
343 lane rv;
344 std::memcpy(&rv, rows + (r * n) + j, sizeof(lane));
345 sum[r] += rv * vv;
346 ov += rv * splat[r];
347 }
348 std::memcpy(out + j, &ov, sizeof(lane));
349 }
350 std::array<T, symmetric_row_group> sums {};
351 for (std::size_t r {0}; r < symmetric_row_group; ++r) {
352 for (std::size_t k {0}; k < width; ++k) {
353 sums[r] += sum[r][k];
354 }
355 }
356 return symmetric_rows_tail(rows, n, v, out, i, done, sums);
357 }
358
373 template <native_scalar T>
374 std::size_t largest_magnitude_vectors(const T* x,
375 std::size_t count,
376 T& best,
377 std::size_t& at)
378 {
379 using lane = typename native_vector<T>::type;
380 using mask = decltype(lane {} > lane {});
381 using index = std::remove_cvref_t<decltype(mask {}[0])>;
382 constexpr std::size_t width {native_vector_bytes / sizeof(T)};
383 constexpr std::size_t ways {4}; // independent chains: one waits a compare and a select per vector
384 if (count < 1 + (ways * width)) {
385 return 1;
386 }
387 const mask magnitude {~std::bit_cast<mask>(T {-0.0} - lane {})}; // every bit but the sign's
388 std::array<lane, ways> top;
389 std::array<mask, ways> where {}; // index 0 in every lane
390 std::array<mask, ways> next; // each chain's next indices
391 for (std::size_t w {0}; w < ways; ++w) {
392 top[w] = best - lane {}; // best in every lane
393 for (std::size_t k {0}; k < width; ++k) {
394 next[w][k] = static_cast<index>(1 + (w * width) + k);
395 }
396 }
397 std::size_t i {1};
398 for (; i + (ways * width) <= count; i += ways * width) { // selects by masks: each lane takes its own side
399# pragma GCC unroll 4 // the chains in registers: GCC keeps a rolled array of vectors in memory
400 for (std::size_t w {0}; w < ways; ++w) {
401 lane v;
402 std::memcpy(&v, x + i + (w * width), sizeof(lane));
403 const lane v_abs {std::bit_cast<lane>(std::bit_cast<mask>(v) & magnitude)};
404 const mask larger {v_abs > top[w]};
405 top[w] = std::bit_cast<lane>((std::bit_cast<mask>(v_abs) & larger) | (std::bit_cast<mask>(top[w]) & ~larger));
406 where[w] = (next[w] & larger) | (where[w] & ~larger);
407 next[w] += static_cast<index>(ways * width);
408 }
409 }
410 for (std::size_t w {0}; w < ways; ++w) { // each lane holds its subsequence's first largest
411 for (std::size_t k {0}; k < width; ++k) {
412 const auto lane_at {static_cast<std::size_t>(where[w][k])};
413 if (top[w][k] > best || (top[w][k] == best && lane_at < at)) {
414 best = top[w][k];
415 at = lane_at;
416 }
417 }
418 }
419 return i;
420 }
421
422#endif
423
425 inline constexpr std::size_t product_rows {8};
426
427#if defined(__GNUC__) && !defined(__clang__)
428
430 template <native_scalar T>
431 std::size_t subtract_scaled_vectors(T* y,
432 const T* x,
433 const T& f,
434 std::size_t count)
435 {
436 using lane = typename native_vector<T>::type;
437 constexpr std::size_t width {native_vector_bytes / sizeof(T)};
438 const std::size_t done {count - (count % width)};
439 for (std::size_t i {0}; i < done; i += width) {
440 lane yv;
441 lane xv;
442 std::memcpy(&yv, y + i, sizeof(lane));
443 std::memcpy(&xv, x + i, sizeof(lane));
444 yv -= f * xv;
445 std::memcpy(y + i, &yv, sizeof(lane));
446 }
447 return done;
448 }
449
458 template <native_scalar T>
459 std::size_t subtract_combined_vectors(T* y,
460 const T* f,
461 const T* x,
462 std::size_t ldx,
463 std::size_t terms,
464 std::size_t count)
465 {
466 using lane = typename native_vector<T>::type;
467 constexpr std::size_t lanes {4};
468 constexpr std::size_t width {native_vector_bytes / sizeof(T)};
469 const std::size_t done {count - (count % (lanes * width))};
470 for (std::size_t i {0}; i < done; i += lanes * width) {
471 std::array<lane, lanes> held {};
472# pragma GCC unroll 4
473 for (std::size_t u {0}; u < lanes; ++u) {
474 std::memcpy(&held[u], y + i + (u * width), sizeof(lane));
475 }
476 for (std::size_t t {0}; t < terms; ++t) {
477 const T factor {f[t]};
478# pragma GCC unroll 4
479 for (std::size_t u {0}; u < lanes; ++u) {
480 lane row;
481 std::memcpy(&row, x + (t * ldx) + i + (u * width), sizeof(lane));
482 held[u] -= factor * row;
483 }
484 }
485# pragma GCC unroll 4
486 for (std::size_t u {0}; u < lanes; ++u) {
487 std::memcpy(y + i + (u * width), &held[u], sizeof(lane));
488 }
489 }
490 return done;
491 }
492
498 template <native_scalar T>
499 std::size_t subtract_rows_vectors(T* y,
500 std::size_t ldy,
501 const T* f,
502 std::size_t ldf,
503 const T* x,
504 std::size_t ldx,
505 std::size_t terms,
506 std::size_t count)
507 {
508 using lane = typename native_vector<T>::type;
509 constexpr std::size_t width {native_vector_bytes / sizeof(T)};
510 const std::size_t done {count - (count % width)};
511 for (std::size_t i {0}; i < done; i += width) {
512 std::array<lane, product_rows> held {};
513# pragma GCC unroll 8
514 for (std::size_t r {0}; r < product_rows; ++r) {
515 std::memcpy(&held[r], y + (r * ldy) + i, sizeof(lane));
516 }
517 for (std::size_t t {0}; t < terms; ++t) {
518 lane row;
519 std::memcpy(&row, x + (t * ldx) + i, sizeof(lane));
520# pragma GCC unroll 8
521 for (std::size_t r {0}; r < product_rows; ++r) {
522 const T factor {f[(r * ldf) + t]};
523 held[r] -= factor * row;
524 }
525 }
526# pragma GCC unroll 8
527 for (std::size_t r {0}; r < product_rows; ++r) {
528 std::memcpy(y + (r * ldy) + i, &held[r], sizeof(lane));
529 }
530 }
531 return done;
532 }
533
535 template <native_scalar T>
536 std::size_t add_scaled_vectors(T* y,
537 const T* x,
538 const T& f,
539 std::size_t count)
540 {
541 using lane = typename native_vector<T>::type;
542 constexpr std::size_t width {native_vector_bytes / sizeof(T)};
543 const std::size_t done {count - (count % width)};
544 for (std::size_t i {0}; i < done; i += width) {
545 lane yv;
546 lane xv;
547 std::memcpy(&yv, y + i, sizeof(lane));
548 std::memcpy(&xv, x + i, sizeof(lane));
549 yv += f * xv;
550 std::memcpy(y + i, &yv, sizeof(lane));
551 }
552 return done;
553 }
554
561 template <split_complex C>
562 std::size_t add_scaled_vectors(C* y,
563 const C* x,
564 const C& f,
565 std::size_t count)
566 {
567 using R = C::value_type;
568 using lane = typename native_vector<R>::type;
569 constexpr std::size_t width {native_vector_bytes / sizeof(R)}; // parts in a vector: an even count
570 auto* const out {reinterpret_cast<R*>(y)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): a complex number is an array of its two parts
571 const auto* const in {reinterpret_cast<const R*>(x)}; // NOLINT(cppcoreguidelines-pro-type-reinterpret-cast): as above
572 lane real_part;
573 lane imaginary_part;
574 lane sign;
575 for (std::size_t k {0}; k < width; ++k) {
576 real_part[k] = f.real();
577 imaginary_part[k] = f.imag();
578 sign[k] = (k % 2 == 0) ? R {-1} : R {1};
579 }
580 const std::size_t parts {2 * count};
581 const std::size_t done {parts - (parts % width)};
582 for (std::size_t i {0}; i < done; i += width) {
583 lane yv;
584 lane xv;
585 std::memcpy(&yv, out + i, sizeof(lane));
586 std::memcpy(&xv, in + i, sizeof(lane));
587 lane swapped;
588 for (std::size_t k {0}; k < width; ++k) {
589 swapped[k] = xv[k ^ 1U];
590 }
591 const lane straight {real_part * xv}; // re·re, re·im
592 const lane crossed {imaginary_part * swapped}; // im·im, im·re
593 const lane signed_crossed {sign * crossed}; // −im·im, +im·re
594 const lane product {straight + signed_crossed}; // re·re − im·im, re·im + im·re
595 yv = yv + product;
596 std::memcpy(out + i, &yv, sizeof(lane));
597 }
598 return done / 2;
599 }
600
601#endif
602
609 template <std::size_t Rows,
610 std::size_t Cols,
611 typename T,
612 typename Full>
613 void through_edge(T* c,
614 std::size_t ldc,
615 std::size_t rows,
616 std::size_t cols,
617 Full full) noexcept
618 {
619 if (rows == Rows && cols == Cols) {
620 full(c, ldc);
621 return;
622 }
623 std::array<T, Rows * Cols> edge {};
624 for (std::size_t i {0}; i < rows; ++i) {
625 for (std::size_t j {0}; j < cols; ++j) {
626 edge[(i * Cols) + j] = c[(i * ldc) + j];
627 }
628 }
629 full(edge.data(), Cols);
630 for (std::size_t i {0}; i < rows; ++i) {
631 for (std::size_t j {0}; j < cols; ++j) {
632 c[(i * ldc) + j] = edge[(i * Cols) + j];
633 }
634 }
635 }
636
638 template <typename T>
639 std::size_t subtract_scaled_vectors(T* /*y*/,
640 const T* /*x*/,
641 const T& /*f*/,
642 std::size_t /*count*/)
643 {
644 return 0;
645 }
646
648 template <typename T,
649 typename Magnitude>
650 std::size_t largest_magnitude_vectors(const T* /*x*/,
651 std::size_t /*count*/,
652 Magnitude& /*best*/,
653 std::size_t& /*at*/)
654 {
655 return 1;
656 }
657
660 template <typename T>
661 std::size_t subtract_combined_vectors(T* /*y*/,
662 const T* /*f*/,
663 const T* /*x*/,
664 std::size_t /*ldx*/,
665 std::size_t /*terms*/,
666 std::size_t /*count*/)
667 {
668 return 0;
669 }
670
673 template <typename T>
674 std::size_t subtract_rows_vectors(T* /*y*/,
675 std::size_t /*ldy*/,
676 const T* /*f*/,
677 std::size_t /*ldf*/,
678 const T* /*x*/,
679 std::size_t /*ldx*/,
680 std::size_t /*terms*/,
681 std::size_t /*count*/)
682 {
683 return 0;
684 }
685
687 template <typename T>
688 std::size_t add_scaled_vectors(T* /*y*/,
689 const T* /*x*/,
690 const T& /*f*/,
691 std::size_t /*count*/)
692 {
693 return 0;
694 }
695
697 template <typename T>
698 [[nodiscard]] T symmetric_row(const T* row,
699 const T* v,
700 const T& f,
701 T* out,
702 std::size_t count)
703 {
704 const T factor {f}; // read once: it may be an element of v
705 T total {};
706 for (std::size_t i {0}; i < count; ++i) {
707 total += row[i] * v[i];
708 out[i] += row[i] * factor;
709 }
710 return total;
711 }
712
714 template <typename T>
715 [[nodiscard]] std::array<T, symmetric_row_group> symmetric_rows(const T* rows,
716 std::size_t n,
717 const T* v,
718 T* out,
719 std::size_t i)
720 {
721 return symmetric_rows_tail(rows, n, v, out, i, 0, {});
722 }
723
725 template <typename T>
726 [[nodiscard]] T dot_product(const T* x,
727 const T* y,
728 std::size_t count)
729 {
730 T total {};
731 for (std::size_t i {0}; i < count; ++i) {
732 total += x[i] * y[i];
733 }
734 return total;
735 }
736
738 template <typename T>
739 [[nodiscard]] T conjugate_dot_product(const T* x,
740 const T* y,
741 std::size_t count)
742 {
743 T total {};
744 for (std::size_t i {0}; i < count; ++i) {
745 total += std::conj(x[i]) * y[i];
746 }
747 return total;
748 }
749
751 template <typename T>
752 void row_dots(const T* m,
753 std::size_t ld,
754 std::size_t rows,
755 const T* x,
756 std::size_t count,
757 T* out)
758 {
759 for (std::size_t r {0}; r < rows; ++r) {
760 out[r] = dot_product(m + (r * ld), x, count);
761 }
762 }
763
771 template <typename T>
772 [[nodiscard]] NUMERICAL_ALWAYS_INLINE T operator_times(const T& a,
773 const T& b)
774 {
775 return a * b;
776 }
777
778#if defined(_LIBCPP_VERSION)
780 template <typename R>
781 [[nodiscard]] NUMERICAL_ALWAYS_INLINE std::complex<R> operator_times(const std::complex<R>& a,
782 const std::complex<R>& b)
783 {
784 const R ac {a.real() * b.real()};
785 const R bd {a.imag() * b.imag()};
786 const R ad {a.real() * b.imag()};
787 const R bc {a.imag() * b.real()};
788 const R x {ac - bd};
789 const R y {ad + bc};
790 if (std::isnan(x) && std::isnan(y)) [[unlikely]] {
791 return a * b;
792 }
793 return {x, y};
794 }
795
796#endif
797
808 template <typename T>
809 void subtract_product(T* y,
810 const T& f,
811 const T* x,
812 std::size_t count)
813 {
814 const T factor {f}; // read once, as a scalar loop would: it may be an element of y's array
815 for (std::size_t i {subtract_scaled_vectors(y, x, factor, count)}; i < count; ++i) {
816 y[i] -= operator_times(factor, x[i]);
817 }
818 }
819
825 template <typename T>
826 void subtract_products(T* y,
827 const T* f,
828 const T* x,
829 std::size_t ldx,
830 std::size_t terms,
831 std::size_t count)
832 {
833 for (std::size_t t {0}; t < terms; ++t) {
834 subtract_product(y, f[t], x + (t * ldx), count);
835 }
836 }
837
846 template <typename T>
847 requires std::is_floating_point_v<T>
848 void subtract_products(T* y,
849 const T* f,
850 const T* x,
851 std::size_t ldx,
852 std::size_t terms,
853 std::size_t count)
854 {
855 std::size_t i {subtract_combined_vectors(y, f, x, ldx, terms, count)};
856 // Blocks of 32, then of 8, 4, 2 and 1, each element a chain of its own through the terms: four multiply-adds a
857 // cycle wait on their last one, and 8 elements held four chains (arm64, an inverse of order 128: 151 → 137 µs);
858 // a narrower rest still takes the terms across its elements at once, where element by element each waited on
859 // the last term (four columns of an LU solve of order 128 on arm64: 112 → 92 µs, the factorisation 75).
860 const auto blocks = [&]<std::size_t Width>(std::integral_constant<std::size_t, Width>) {
861 for (; i + Width <= count; i += Width) {
862 std::array<T, Width> held {};
863 std::copy_n(y + i, Width, held.begin());
864 for (std::size_t t {0}; t < terms; ++t) {
865 const T factor {f[t]};
866 const T* const row {x + (t * ldx) + i};
867 for (std::size_t u {0}; u < Width; ++u) {
868 held[u] -= factor * row[u];
869 }
870 }
871 std::copy_n(held.begin(), Width, y + i);
872 }
873 };
874 blocks(std::integral_constant<std::size_t, 32> {});
875 blocks(std::integral_constant<std::size_t, 8> {});
876 blocks(std::integral_constant<std::size_t, 4> {});
877 blocks(std::integral_constant<std::size_t, 2> {});
878 blocks(std::integral_constant<std::size_t, 1> {});
879 }
880
888 template <typename T>
889 requires std::is_floating_point_v<T>
890 void subtract_rows(T* y,
891 std::size_t ldy,
892 const T* f,
893 std::size_t ldf,
894 const T* x,
895 std::size_t ldx,
896 std::size_t terms,
897 std::size_t count)
898 {
899 constexpr std::size_t width {4}; // four columns of eight rows: 32 values, which arm64's registers hold
900 std::size_t i {subtract_rows_vectors(y, ldy, f, ldf, x, ldx, terms, count)};
901 const auto columns = [&]<std::size_t Width>(std::integral_constant<std::size_t, Width>) {
902 for (; i + Width <= count; i += Width) {
903 std::array<std::array<T, Width>, product_rows> held {};
904 for (std::size_t r {0}; r < product_rows; ++r) {
905 std::copy_n(y + (r * ldy) + i, Width, held[r].begin());
906 }
907 for (std::size_t t {0}; t < terms; ++t) {
908 const T* const row {x + (t * ldx) + i};
909 for (std::size_t r {0}; r < product_rows; ++r) {
910 const T factor {f[(r * ldf) + t]};
911 for (std::size_t u {0}; u < Width; ++u) {
912 held[r][u] -= factor * row[u];
913 }
914 }
915 }
916 for (std::size_t r {0}; r < product_rows; ++r) {
917 std::copy_n(held[r].begin(), Width, y + (r * ldy) + i);
918 }
919 }
920 };
921 columns(std::integral_constant<std::size_t, width> {});
922 columns(std::integral_constant<std::size_t, 1> {});
923 }
924
926 template <typename T>
927 void subtract_product(T* y,
928 const T* x,
929 const T& f,
930 std::size_t count)
931 {
932 const T factor {f}; // read once, as a scalar loop would: it may be an element of y's array
933 for (std::size_t i {subtract_scaled_vectors(y, x, factor, count)}; i < count; ++i) {
934 y[i] -= operator_times(x[i], factor);
935 }
936 }
937
940 template <split_complex C>
941 void complex_scaled_add(C& y,
942 const C& f,
943 const C& x)
944 {
945 using R = C::value_type;
946 const R ac {f.real() * x.real()};
947 const R bd {f.imag() * x.imag()};
948 const R ad {f.real() * x.imag()};
949 const R bc {f.imag() * x.real()};
950 const R re {ac - bd};
951 const R im {ad + bc};
952 y = C {y.real() + re, y.imag() + im};
953 }
954
957 template <split_complex C>
958 void add_product(C* y,
959 const C& f,
960 const C* x,
961 std::size_t count)
962 {
963 const C factor {f}; // read once, as a scalar loop would: it may be an element of y's array
964 for (std::size_t i {add_scaled_vectors(y, x, factor, count)}; i < count; ++i) {
965 complex_scaled_add(y[i], factor, x[i]);
966 }
967 }
968
970 template <typename T>
971 void add_one_product(T& y,
972 const T& f,
973 const T& x)
974 {
975 y += operator_times(f, x);
976 }
977
979 template <split_complex C>
980 void add_one_product(C& y,
981 const C& f,
982 const C& x)
983 {
984 complex_scaled_add(y, f, x);
985 }
986
988 template <typename T>
989 void add_product(T* y,
990 const T& f,
991 const T* x,
992 std::size_t count)
993 {
994 const T factor {f}; // read once, as a scalar loop would: it may be an element of y's array
995 for (std::size_t i {add_scaled_vectors(y, x, factor, count)}; i < count; ++i) {
996 y[i] += operator_times(factor, x[i]);
997 }
998 }
999} // namespace numerical::detail
1000
1001#endif // NUMERICAL_SIMD_HPP
numerical::detail::split_complex
A complex number of float or double, which gemm_split multiplies through the real kernel and the kern...
Definition simd.hpp:64
numerical::detail::symmetric_row
T symmetric_row(const T *row, const T *v, const T &f, T *out, std::size_t count)
Σ row[i] · v[i] and out[i] += row[i] · f, in order, for a type no native vector holds.
Definition simd.hpp:698
numerical::detail::subtract_rows
void subtract_rows(double *y, std::size_t ldy, const double *f, std::size_t ldf, const double *x, std::size_t ldx, std::size_t terms, std::size_t count)
subtract_rows for double, the units chosen as subtract_products chooses them: each element the value ...
Definition gemm.hpp:217
numerical::detail::conjugate_dot_product
T conjugate_dot_product(const T *x, const T *y, std::size_t count)
Σ conj(x[i]) · y[i], in order, where no native vector holds the parts.
Definition simd.hpp:739
numerical::detail::transpose_tile
constexpr std::size_t transpose_tile
The side of the tiles a transposing copy goes by (an einsum permutation of two axes,...
Definition simd.hpp:31
numerical::detail::add_one_product
void add_one_product(T &y, const T &f, const T &x)
y += f · x of one element, as add_product computes each of its elements.
Definition simd.hpp:971
numerical::detail::subtract_rows_vectors
std::size_t subtract_rows_vectors(T *, std::size_t, const T *, std::size_t, const T *, std::size_t, std::size_t, std::size_t)
No leading columns in vectors, for a type no native vector holds or a compiler that vectorises the sc...
Definition simd.hpp:674
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::operator_times
T operator_times(const T &a, const T &b)
a · b, as std::operator* computes it, in line. libc++'s complex product (the four products and their ...
Definition simd.hpp:772
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::detail::add_product
void add_product(double *y, const double &f, const double *x, std::size_t count)
y[i] += f · x[i] for double, as subtract_product.
Definition gemm.hpp:308
numerical::detail::subtract_scaled_vectors
std::size_t subtract_scaled_vectors(T *, const T *, const T &, std::size_t)
No leading elements in vectors, for a type no native vector holds.
Definition simd.hpp:639
numerical::detail::row_dots
void row_dots(const std::complex< double > *m, std::size_t ld, std::size_t rows, const std::complex< double > *x, std::size_t count, std::complex< double > *out)
row_dots for std::complex<double>, through complex_row_dots.
Definition gemm.hpp:432
numerical::detail::subtract_products
void subtract_products(double *y, const double *f, const double *x, std::size_t ldx, std::size_t terms, std::size_t count)
y[i] -= Σₜ f[t] · x[(t · ldx) + i] for double, the units chosen as subtract_product chooses them,...
Definition gemm.hpp:194
numerical::detail::subtract_product
void subtract_product(double *y, const double &f, const double *x, std::size_t count)
y[i] −= f · x[i] for double on the widest x86-64 units the matrix product's kernel may use (gemm_path...
Definition gemm.hpp:163
numerical::detail::symmetric_rows
std::array< T, symmetric_row_group > symmetric_rows(const T *rows, std::size_t n, const T *v, T *out, std::size_t i)
symmetric_rows one column at a time, for a type no native vector holds.
Definition simd.hpp:715
numerical::detail::real_part
double real_part(const T &x)
The real part of a real (or general) scalar as a double.
Definition decomposition.hpp:280
numerical::detail::subtract_combined_vectors
std::size_t subtract_combined_vectors(T *, const T *, const T *, std::size_t, std::size_t, std::size_t)
No leading blocks in vectors, for a type no native vector holds or a compiler that vectorises the sca...
Definition simd.hpp:661
numerical::detail::dot_product
double dot_product(const double *x, const double *y, std::size_t count)
Σ x[i] · y[i] for double, as subtract_product chooses its units (in another order of summation there,...
Definition gemm.hpp:369
numerical::detail::complex_scaled_add
void complex_scaled_add(C &y, const C &f, const C &x)
y += f · x of one complex number, its part products rounded apart as libc++'s operator computes them ...
Definition simd.hpp:941
numerical::detail::product_rows
constexpr std::size_t product_rows
The rows of y subtract_rows takes at once.
Definition simd.hpp:425
numerical::detail::add_scaled_vectors
std::size_t add_scaled_vectors(T *, const T *, const T &, std::size_t)
No leading elements in vectors, for a type no native vector holds.
Definition simd.hpp:688
numerical::detail::through_edge
void through_edge(T *c, std::size_t ldc, std::size_t rows, std::size_t cols, Full full) noexcept
One tile of a matrix product kernel, rows × cols of C (at most Rows × Cols, leading dimension ldc),...
Definition simd.hpp:613
numerical::detail::largest_magnitude_vectors
std::size_t largest_magnitude_vectors(const T *, std::size_t, Magnitude &, std::size_t &)
Only x[0], the one best and at come in with, for a type no native vector holds.
Definition simd.hpp:650
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:726
numerical::full
tensor< T > full(std::vector< std::size_t > shape, const T &value)
A tensor of shape with every element value.
Definition tensor_ops.hpp:816
NUMERICAL_ALWAYS_INLINE
#define NUMERICAL_ALWAYS_INLINE
Definition simd.hpp:22
Generated by doxygen 1.9.8