/* References: saffah_codex_6a_agg3 https://duck.ac/submission/121243,
reused the matrix engine and inherited references below. No independent
license declared.
Idea: Address output rows in pairs: store the second row via base+stride,
then advance the base by two strides. This replaces five serial pointer adds
with two updates in each six-row microkernel, preserving all store order.
Purpose: Reduce output-address arithmetic after the FMA loop.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/121181,
reused the matrix engine and all inherited source references below.
No independent license declared.
Idea: Specialize recursion for the known 4096, 2048, 1024 and 512 sizes. Each
level can discard unreachable size branches and use constant transform strides
and extents, preserving the same product order and reconstruction schedule.
Purpose: Remove dynamic size arithmetic and dispatch above the matrix leaves.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/121167,
reused the matrix engine and all inherited references below. No independent
license declared.
Idea: Every recursive multiplication in this fixed 4096-square task terminates
at a 512-square leaf. Compile that leaf extent as a constant to simplify packing
bounds, panel strides, tail selection and streaming-store conditions.
Purpose: Remove general-size arithmetic from the remaining leaf hot loops.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/121142,
reused the matrix engine and all inherited references below. No independent
license declared.
Idea: Specialize leaf reconstruction independently for each of the three
levels active at entry. Deferred tasks can retire but cannot start within a
leaf; omit checks for the levels that were absent when this leaf began.
Purpose: Extend the measured benefit of inactive-task specialization to leaves
with only one or two active reconstruction levels.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/120888,
reused the current pipelined matrix engine and all inherited source credits below.
No independent license declared.
Idea: Select a leaf implementation once based on whether deferred reconstruction
is active. No reconstruction task can begin inside a leaf, so a leaf entering
with none can omit every per-microkernel task check and row-frontier update.
The other implementation retains all existing task progression and ordering.
Purpose: Remove repeated inactive-task dispatch from most microkernel calls.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/120615,
reused the current matrix engine and plain-kernel software B-load pipeline.
All earlier source attributions remain below; no independent license declared.
Idea: Apply the one-step B-load pipeline to the A-packing microkernel too.
Its spare register is only needed for packing after the last inner FMA, so
alternating that register with the current B vector preserves all values.
Purpose: Hide B load latency in the remaining carried-packing kernel calls.
*/
/* References: pdoom https://duck.ac/submission/118506: reused FP64
recursive matrix engine. saffah_cc_v41_agg1 https://duck.ac/submission/119315:
reused plain-kernel B-load software pipeline and deferred prefetch scheduling.
saffah_codex_6a_agg3 https://duck.ac/submission/120547: B-panel reuse and merge64.
All inherited source credits retained below; no independent licenses declared.
Idea: Apply the one-step-ahead B-load scheduling to the larger matrix kernel,
using the fixed-size output path and preserving each FMA accumulation order.
Purpose: Hide B load-to-use latency across more than100000 microkernel calls.
*/
/* References:
[1] saffah_codex_6a_agg3, https://duck.ac/submission/120521.
Reused the complete FP64 matrix engine; all inherited source citations remain
below. No separate license is declared on the public source.
Idea: Reconstruct64 positions per bottom merge step instead of32. Both divide
the512-wide rows, preserving frontiers and operation grouping while amortizing
step checks and cursor updates across twice as much contiguous memory traffic.
Purpose: Measure lower reconstruction bookkeeping in the fused compute schedule.
*/
/* References:
[1] pdoom, https://duck.ac/submission/118506.
Reused the4k FP64 engine and its deferred root reconstruction, preserving its
complete inherited attribution list below. No separate license is declared.
[2] saffah_cc_v41_agg1, https://duck.ac/submission/119315.
Copied the1k fused packed-B helpers and two-reuse leaf ordering.
[3] saffah_codex_6a_agg3, https://duck.ac/submission/120514.
Reused our validated2k integration of those helpers; no separate license.
Idea: Apply two packed-B reuse positions to every bottom group of the4k
recursion, removing repeated source traversal when a pure panel is still packed.
Purpose: Measure the accumulated traffic reduction at the largest FP64 size.
*/
/*
Credits for inherited implementation and ideas:
- saffah_cc_v41_agg1, https://duck.ac/submission/105696: immediate engine
base; safely extend the upper deferred-combine window across siblings.
- saffah_codex_6s_agg2, https://duck.ac/submission/104821: preceding engine,
compiler model and upper workspace layout, including its inherited work.
- saffah_cc_v41_agg1, https://duck.ac/submission/101193: two-level deferred
reconstruction, product order and panel-layout improvements.
- saffah_cc_v41_agg1, https://duck.ac/submission/103010: produced-row
frontier and early reconstruction; saffah_codex_6s_agg2,
https://duck.ac/submission/103224: its4k integration and leaf ordering.
- saffah_codex_6s_agg2, https://duck.ac/submission/103194: eight-K plain
microkernel; https://duck.ac/submission/103275: power-of-two addressing;
https://duck.ac/submission/104809: leaf scratch-plane layout.
- saffah_cc_v41_agg1, https://duck.ac/submission/100776: dedicated two-row
tail; https://duck.ac/submission/100511 and https://duck.ac/submission/100475:
carry next-A packing in the FMA loop, double-buffer A and grouped B packing.
- saffah_cc_v41_260924, https://duck.ac/submission/99846: fused operand sums
and packing, one-pass sums and reconstruction.
- saffah_cc_v41_agg1, https://duck.ac/submission/99875 and
https://duck.ac/submission/99950: non-temporal intermediate stores.
- saffah_cc_v41_260924, https://duck.ac/submission/99141, and
saffah_cc_v41_agg1, https://duck.ac/submission/99167:6x8 AVX2 blocking,
explicit vector memory instructions, prefetch and packing-layout baseline.
- Earlier pdoom work: hybrid Winograd/Strassen recurrence, output-backed
scratch placement, delayed reconstruction, merged32-row tail and2x16
kernel, and page-spacing experiments.
New changes in this version:
1. Specialize the full6x8 kernels and hot upper reconstruction geometry for
the fixed4096 problem. Two-row output has its separate2x16 kernel.
2. The root's last product M3 only contributes to C12. Once the preceding
M1 is ready, reconstruct C11,C21,C22 and the C12 base (M1+M6)+M5 during
M3's remaining five1024 subproblems. Retired M5 storage holds that base.
3. Add the base while reconstructing M3's four quadrants. The root task
processes32 elements per step only when bottom reconstruction is inactive.
This moves320MiB of logical root traffic into the last multiplication
subtree; total nominal traffic and arithmetic count do not decrease.
All arithmetic remains binary64, with the original expression grouping.
No input quantization. Specialized to the problem's n=4096 contract.
*/
#pragma GCC optimize("O3,unroll-loops,no-strict-aliasing")
#include <sys/mman.h>
#pragma GCC target("avx2,fma,arch=skylake")
#include <immintrin.h>
#ifndef B2M_NO_BARA
#define B2M_NO_BARA 1
#endif
#define MD9_T2BASE 1
#include <stddef.h>
#ifndef BASE
#define BASE 512
#endif
#ifndef MC
#define MC 30
#endif
#ifndef NTMIN
#define NTMIN 512
#endif
#ifndef UPPER_PAGE_PAD
#define UPPER_PAGE_PAD 15872
#endif
#ifndef LEAF_PAD
#define LEAF_PAD 1152
#endif
#ifndef SUMNT
#define SUMNT 1024
#endif
static constexpr int MR=6,NR=8;
alignas(4096) static double workspace[48*1024*1024];
alignas(64) static double ap[(MC+MR)*(BASE+24)],ap2[(MC+MR)*(BASE+24)],bp[BASE*(BASE+24)];
/* m2s/md5 fused-pack support: a zero row (for pad rows and for pure copies) and the
two sign masks (+0.0 / -0.0) used to apply Strassen's operand sign in-register. */
alignas(64) static double pk_zero[BASE+8];
alignas(64) static double pk_mp[4]={0.0,0.0,0.0,0.0}, pk_mn[4]={-0.0,-0.0,-0.0,-0.0};
struct Mat { const double *x; int ld; const double *y=nullptr; int ly=0; int sign=0; };
static inline __m256d ld(Mat a,int i,int j) {
__m256d x=_mm256_load_pd(a.x+(size_t)i*a.ld+j);
if(a.sign) {
__m256d y=_mm256_load_pd(a.y+(size_t)i*a.ly+j);
x=a.sign>0?_mm256_add_pd(x,y):_mm256_sub_pd(x,y);
}
return x;
}
static inline double scalar(Mat a,int i,int k) {
double x=a.x[(size_t)i*a.ld+k];
if(a.sign)x+=a.sign>0?a.y[(size_t)i*a.ly+k]:-a.y[(size_t)i*a.ly+k];
return x;
}
template<int SG> static __attribute__((always_inline)) inline void kernel(int n,const double *a,const double *b,double *c,int ldc,int rows, const double *pks,const double *pky,double *pkd) {
double *out=c;
ptrdiff_t stride=(ptrdiff_t)ldc*8;
int k=n;
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
"vxorpd %%ymm4, %%ymm4, %%ymm4\n\t"
"vxorpd %%ymm5, %%ymm5, %%ymm5\n\t"
"vxorpd %%ymm6, %%ymm6, %%ymm6\n\t"
"vxorpd %%ymm7, %%ymm7, %%ymm7\n\t"
"vxorpd %%ymm8, %%ymm8, %%ymm8\n\t"
"vxorpd %%ymm9, %%ymm9, %%ymm9\n\t"
"vxorpd %%ymm10, %%ymm10, %%ymm10\n\t"
"vxorpd %%ymm11, %%ymm11, %%ymm11\n\t"
"vmovapd 0(%[b]), %%ymm12\n\t"
".p2align 6\n\t"
"1:\n\t"
"vmovapd 64(%[b]), %%ymm15\n\t"
"vmovapd 32(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 128(%[b]), %%ymm12\n\t"
"vmovapd 96(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 704(%[b])\n\t"
"vmovapd 192(%[b]), %%ymm15\n\t"
"vmovapd 160(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 256(%[b]), %%ymm12\n\t"
"vmovapd 224(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 768(%[b])\n\t"
"prefetcht0 832(%[b])\n\t"
"prefetcht0 896(%[b])\n\t"
"prefetcht0 960(%[b])\n\t"
"vmovapd 320(%[b]), %%ymm15\n\t"
"vmovapd 288(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 384(%[b]), %%ymm12\n\t"
"vmovapd 352(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 448(%[b]), %%ymm15\n\t"
"vmovapd 416(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 512(%[b]), %%ymm12\n\t"
"vmovapd 480(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovupd 0(%[pks]), %%ymm15\n\t"
".if %c[sg] == 1\n\t"
"vaddpd 0(%[pky]), %%ymm15, %%ymm15\n\t"
".elseif %c[sg] == -1\n\t"
"vsubpd 0(%[pky]), %%ymm15, %%ymm15\n\t"
".endif\n\t"
"vmovupd %%ymm15, 0(%[pkd])\n\t"
"add $32, %[pks]\n\t"
"add $32, %[pky]\n\t"
"add $32, %[pkd]\n\t"
"add $64, %[a]\n\t"
"add $512, %[b]\n\t"
"sub $8, %[k]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"vmovntpd %%ymm2, (%[out],%[stride])\n\t"
"vmovntpd %%ymm3, 32(%[out],%[stride])\n\t"
"lea (%[out],%[stride],2), %[out]\n\t"
"vmovntpd %%ymm4, (%[out])\n\t"
"vmovntpd %%ymm5, 32(%[out])\n\t"
"vmovntpd %%ymm6, (%[out],%[stride])\n\t"
"vmovntpd %%ymm7, 32(%[out],%[stride])\n\t"
"lea (%[out],%[stride],2), %[out]\n\t"
"vmovntpd %%ymm8, (%[out])\n\t"
"vmovntpd %%ymm9, 32(%[out])\n\t"
"vmovntpd %%ymm10, (%[out],%[stride])\n\t"
"vmovntpd %%ymm11, 32(%[out],%[stride])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k), [pks] "+&r"(pks), [pky] "+&r"(pky), [pkd] "+&r"(pkd)
: [sg] "i"(SG), [stride] "r"(stride), [o0] "i"(0*(BASE+24)*8), [o1] "i"(1*(BASE+24)*8), [o2] "i"(2*(BASE+24)*8), [o3] "i"(3*(BASE+24)*8), [o4] "i"(4*(BASE+24)*8), [o5] "i"(5*(BASE+24)*8)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
static __attribute__((always_inline)) inline void kernel_plain(int n,const double *a,const double *b,double *c,int ldc,int rows) {
double *out=c;
ptrdiff_t stride=(ptrdiff_t)ldc*8;
int k=n;
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
"vxorpd %%ymm4, %%ymm4, %%ymm4\n\t"
"vxorpd %%ymm5, %%ymm5, %%ymm5\n\t"
"vxorpd %%ymm6, %%ymm6, %%ymm6\n\t"
"vxorpd %%ymm7, %%ymm7, %%ymm7\n\t"
"vxorpd %%ymm8, %%ymm8, %%ymm8\n\t"
"vxorpd %%ymm9, %%ymm9, %%ymm9\n\t"
"vxorpd %%ymm10, %%ymm10, %%ymm10\n\t"
"vxorpd %%ymm11, %%ymm11, %%ymm11\n\t"
"vmovapd 0(%[b]), %%ymm12\n\t"
".p2align 6\n\t"
"1:\n\t"
"vmovapd 64(%[b]), %%ymm15\n\t"
"vmovapd 32(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 128(%[b]), %%ymm12\n\t"
"vmovapd 96(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 704(%[b])\n\t"
"vmovapd 192(%[b]), %%ymm15\n\t"
"vmovapd 160(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 256(%[b]), %%ymm12\n\t"
"vmovapd 224(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"prefetcht0 768(%[b])\n\t"
"prefetcht0 832(%[b])\n\t"
"prefetcht0 896(%[b])\n\t"
"prefetcht0 960(%[b])\n\t"
"vmovapd 320(%[b]), %%ymm15\n\t"
"vmovapd 288(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+32(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 384(%[b]), %%ymm12\n\t"
"vmovapd 352(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+40(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 448(%[b]), %%ymm15\n\t"
"vmovapd 416(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+48(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"vmovapd 512(%[b]), %%ymm12\n\t"
"vmovapd 480(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vbroadcastsd %c[o2]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm4\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm5\n\t"
"vbroadcastsd %c[o3]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm6\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm7\n\t"
"vbroadcastsd %c[o4]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm8\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm9\n\t"
"vbroadcastsd %c[o5]+56(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm15, %%ymm14, %%ymm10\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm11\n\t"
"add $64, %[a]\n\t"
"add $512, %[b]\n\t"
"sub $8, %[k]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"vmovntpd %%ymm2, (%[out],%[stride])\n\t"
"vmovntpd %%ymm3, 32(%[out],%[stride])\n\t"
"lea (%[out],%[stride],2), %[out]\n\t"
"vmovntpd %%ymm4, (%[out])\n\t"
"vmovntpd %%ymm5, 32(%[out])\n\t"
"vmovntpd %%ymm6, (%[out],%[stride])\n\t"
"vmovntpd %%ymm7, 32(%[out],%[stride])\n\t"
"lea (%[out],%[stride],2), %[out]\n\t"
"vmovntpd %%ymm8, (%[out])\n\t"
"vmovntpd %%ymm9, 32(%[out])\n\t"
"vmovntpd %%ymm10, (%[out],%[stride])\n\t"
"vmovntpd %%ymm11, 32(%[out],%[stride])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k)
: [stride] "r"(stride), [o0] "i"(0*(BASE+24)*8), [o1] "i"(1*(BASE+24)*8), [o2] "i"(2*(BASE+24)*8), [o3] "i"(3*(BASE+24)*8), [o4] "i"(4*(BASE+24)*8), [o5] "i"(5*(BASE+24)*8)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
/* ===== mdA_dp2: deferred h=1024 combine (Winograd form, exact same arithmetic) ===== */
struct Pend2 { double *c; int ldc; const double *m; int h; int pz; int pos; bool active; const double *add=nullptr; };
static Pend2 pend2={0,0,0,0,0,0,false};
static constexpr int STEP2=16;
template<int H=1024> static inline void pend2_prefetch(void){
if(!pend2.active)return;
int i=pend2.pos >> __builtin_ctz((unsigned)H), j=pend2.pos & (H-1);
const double *c=pend2.c+(size_t)i*pend2.ldc+j;
const double *m=pend2.m+pend2.pos;
for(int t=0;t<STEP2;t+=8) {
_mm_prefetch((const char*)(c+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+H+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)H*pend2.ldc+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)H*pend2.ldc+H+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+(H*H+UPPER_PAGE_PAD)+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+2*(H*H+UPPER_PAGE_PAD)+t),_MM_HINT_T1);
}
}
template<int H=1024,bool ADD=false> static inline void pend2_step(void){
int i=pend2.pos >> __builtin_ctz((unsigned)H), j=pend2.pos & (H-1);
for(int t=0;t<STEP2;t+=4) {
double *c11=pend2.c+(size_t)i*pend2.ldc+j+t,*c12=c11+H,
*c21=c11+(size_t)H*pend2.ldc,*c22=c21+H;
const double *m=pend2.m+pend2.pos+t;
__m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),
m3=_mm256_load_pd(c12),m4=_mm256_load_pd(c22),
m5=_mm256_load_pd(m),m6=_mm256_load_pd(m+(H*H+UPPER_PAGE_PAD)),m7=_mm256_load_pd(m+2*(H*H+UPPER_PAGE_PAD));
__m256d u2=_mm256_add_pd(m1,m6),u3=_mm256_add_pd(u2,m7),u4=_mm256_add_pd(u2,m5);
__m256d r11=_mm256_add_pd(m1,m2),r12=_mm256_add_pd(u4,m3),
r21=_mm256_sub_pd(u3,m4),r22=_mm256_add_pd(u3,m5);
if(ADD) {
const double *p=pend2.add+(size_t)i*2048+j+t;
r11=_mm256_add_pd(_mm256_load_pd(p),r11);
r12=_mm256_add_pd(_mm256_load_pd(p+H),r12);
r21=_mm256_add_pd(_mm256_load_pd(p+(size_t)H*2048),r21);
r22=_mm256_add_pd(_mm256_load_pd(p+(size_t)H*2048+H),r22);
}
_mm256_store_pd(c11,r11);_mm256_store_pd(c12,r12);
_mm256_store_pd(c21,r21);_mm256_store_pd(c22,r22);
}
pend2.pos+=STEP2;
if(pend2.pos>=H*H)pend2.active=false;
}
static void pend2_finish(void){
if(pend2.add) {while(pend2.active)pend2_step<1024,true>();}
else if(pend2.h==1024) {while(pend2.active)pend2_step<1024>();}
else {while(pend2.active)pend2_step<2048>();}
}
struct Pending { double *c; int ldc; const double *m; int pos=0; bool active=false; int rowlim=0; bool prog=false; };
static Pending pending;
alignas(4096) static double pending_work[3*(BASE*BASE+LEAF_PAD)+LEAF_PAD];
static unsigned pending_tick=0;
static constexpr int STEP=64, PERIOD=1;
static inline void pending_prefetch() {
if(!pending.active)return;
if(pending.prog && pending.pos/BASE>=pending.rowlim)return;
int i=pending.pos/BASE,j=pending.pos%BASE;
const double *c=pending.c+(size_t)i*pending.ldc+j;
const double *m=pending.m+pending.pos;
for(int t=0;t<STEP;t+=8) {
_mm_prefetch((const char*)(c+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+BASE+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)BASE*pending.ldc+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)BASE*pending.ldc+BASE+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+(BASE*BASE+LEAF_PAD)+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+2*(BASE*BASE+LEAF_PAD)+t),_MM_HINT_T1);
}
}
static inline void pending_step() {
if(pending.prog && pending.pos/BASE>=pending.rowlim)return;
int i=pending.pos/BASE,j=pending.pos%BASE;
for(int t=0;t<STEP;t+=4) {
double *c11=pending.c+(size_t)i*pending.ldc+j+t,
*c12=c11+BASE,*c21=c11+(size_t)BASE*pending.ldc,*c22=c21+BASE;
const double *m=pending.m+pending.pos+t;
__m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),
m3=_mm256_load_pd(c12),m4=_mm256_load_pd(c22),
m5=_mm256_load_pd(m),m6=_mm256_load_pd(m+(BASE*BASE+LEAF_PAD)),m7=_mm256_load_pd(m+2*(BASE*BASE+LEAF_PAD));
__m256d r11=_mm256_add_pd(_mm256_sub_pd(_mm256_add_pd(m1,m4),m5),m7),
r12=_mm256_add_pd(m3,m5),r21=_mm256_add_pd(m2,m4),
r22=_mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(m1,m2),m3),m6);
_mm256_store_pd(c11,r11);_mm256_store_pd(c12,r12);
_mm256_store_pd(c21,r21);_mm256_store_pd(c22,r22);
}
pending.pos+=STEP;
if(pending.pos==BASE*BASE)pending.active=false;
}
static void pending_finish() {
pending.prog=false;
while(pending.active)pending_step();
}
static __attribute__((always_inline)) inline void kernel2(int n,const double *a,const double *b,double *c,int ldc) {
int k=n;
double *out=c;
ptrdiff_t stride=(ptrdiff_t)ldc*8;
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
".p2align 6\n\t" "1:\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 704(%[b])\n\t"
"vmovapd 0(%[b]), %%ymm12\n\t"
"vmovapd 32(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vmovapd 64(%[b]), %%ymm12\n\t"
"vmovapd 96(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vmovapd 128(%[b]), %%ymm12\n\t"
"vmovapd 160(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"vmovapd 192(%[b]), %%ymm12\n\t"
"vmovapd 224(%[b]), %%ymm13\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm0\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm1\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm14\n\t"
"vfmadd231pd %%ymm12, %%ymm14, %%ymm2\n\t"
"vfmadd231pd %%ymm13, %%ymm14, %%ymm3\n\t"
"add $32, %[a]\n\t"
"add $256, %[b]\n\t"
"sub $4, %[k]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, (%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm2, (%[out])\n\t"
"vmovntpd %%ymm3, 32(%[out])\n\t"
: [a] "+&r"(a), [b] "+&r"(b), [out] "+&r"(out), [k] "+&r"(k)
: [stride] "r"(stride), [o0] "i"(0*(BASE+24)*8), [o1] "i"(1*(BASE+24)*8)
: "memory", "cc", "ymm0","ymm1","ymm2","ymm3","ymm12","ymm13","ymm14");
}
static __attribute__((noinline)) void kernel_tail16(int n,const double *a,const double *b,double *out,int ldc) {
const double *b1=b+8*(n+1);
ptrdiff_t stride=(ptrdiff_t)ldc*8;
int kk=n;
asm volatile (
"vxorpd %%ymm0, %%ymm0, %%ymm0\n\t"
"vxorpd %%ymm1, %%ymm1, %%ymm1\n\t"
"vxorpd %%ymm2, %%ymm2, %%ymm2\n\t"
"vxorpd %%ymm3, %%ymm3, %%ymm3\n\t"
"vxorpd %%ymm4, %%ymm4, %%ymm4\n\t"
"vxorpd %%ymm5, %%ymm5, %%ymm5\n\t"
"vxorpd %%ymm6, %%ymm6, %%ymm6\n\t"
"vxorpd %%ymm7, %%ymm7, %%ymm7\n\t"
".p2align 6\n\t"
"1:\n\t"
"prefetcht0 512(%[b])\n\t"
"prefetcht0 512(%[b1])\n\t"
"vbroadcastsd %c[o0]+0(%[a]), %%ymm12\n\t"
"vbroadcastsd %c[o1]+0(%[a]), %%ymm13\n\t"
"vmovapd 0(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm0\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm4\n\t"
"vmovapd 32(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm5\n\t"
"vmovapd 0(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm6\n\t"
"vmovapd 32(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"prefetcht0 576(%[b])\n\t"
"prefetcht0 576(%[b1])\n\t"
"vbroadcastsd %c[o0]+8(%[a]), %%ymm12\n\t"
"vbroadcastsd %c[o1]+8(%[a]), %%ymm13\n\t"
"vmovapd 64(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm0\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm4\n\t"
"vmovapd 96(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm5\n\t"
"vmovapd 64(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm6\n\t"
"vmovapd 96(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"prefetcht0 640(%[b])\n\t"
"prefetcht0 640(%[b1])\n\t"
"vbroadcastsd %c[o0]+16(%[a]), %%ymm12\n\t"
"vbroadcastsd %c[o1]+16(%[a]), %%ymm13\n\t"
"vmovapd 128(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm0\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm4\n\t"
"vmovapd 160(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm5\n\t"
"vmovapd 128(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm6\n\t"
"vmovapd 160(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"prefetcht0 704(%[b])\n\t"
"prefetcht0 704(%[b1])\n\t"
"vbroadcastsd %c[o0]+24(%[a]), %%ymm12\n\t"
"vbroadcastsd %c[o1]+24(%[a]), %%ymm13\n\t"
"vmovapd 192(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm0\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm4\n\t"
"vmovapd 224(%[b]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm1\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm5\n\t"
"vmovapd 192(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm2\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm6\n\t"
"vmovapd 224(%[b1]), %%ymm14\n\t"
"vfmadd231pd %%ymm14, %%ymm12, %%ymm3\n\t"
"vfmadd231pd %%ymm14, %%ymm13, %%ymm7\n\t"
"add $32, %[a]\n\t"
"add $256, %[b]\n\t"
"add $256, %[b1]\n\t"
"sub $4, %[kk]\n\t"
"jnz 1b\n\t"
"vmovntpd %%ymm0, 0(%[out])\n\t"
"vmovntpd %%ymm1, 32(%[out])\n\t"
"vmovntpd %%ymm2, 64(%[out])\n\t"
"vmovntpd %%ymm3, 96(%[out])\n\t"
"add %[stride], %[out]\n\t"
"vmovntpd %%ymm4, 0(%[out])\n\t"
"vmovntpd %%ymm5, 32(%[out])\n\t"
"vmovntpd %%ymm6, 64(%[out])\n\t"
"vmovntpd %%ymm7, 96(%[out])\n\t"
: [a] "+&r"(a),[b] "+&r"(b),[b1] "+&r"(b1),[out] "+&r"(out),[kk] "+&r"(kk)
: [stride] "r"(stride),[o0] "i"(0),[o1] "i"((BASE+24)*8)
: "memory","cc","ymm0","ymm1","ymm2","ymm3","ymm4","ymm5","ymm6","ymm7","ymm8","ymm9","ymm10","ymm11","ymm12","ymm13","ymm14","ymm15");
}
// Root M3 is last and contributes only to C12. The other three outputs and
// the C12 base (M1+M6)+M5 are independent of M3.
struct RootPending { double*c;int ldc;double*m;int pos;bool active; };
static RootPending root_pending={nullptr,0,nullptr,0,false};
static constexpr int ROOT_STEP=32;
static inline void root_prefetch() {
if(!root_pending.active)return;
int i=root_pending.pos>>11,j=root_pending.pos&2047;
const double*c=root_pending.c+(size_t)i*root_pending.ldc+j;
const double*m=root_pending.m+root_pending.pos;
constexpr size_t PZ=2048*2048+UPPER_PAGE_PAD;
for(int t=0;t<ROOT_STEP;t+=8){
_mm_prefetch((const char*)(c+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)2048*root_pending.ldc+t),_MM_HINT_T1);
_mm_prefetch((const char*)(c+(size_t)2048*root_pending.ldc+2048+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+PZ+t),_MM_HINT_T1);
_mm_prefetch((const char*)(m+2*PZ+t),_MM_HINT_T1);
}
}
static inline void root_step() {
int i=root_pending.pos>>11,j=root_pending.pos&2047;
constexpr size_t PZ=2048*2048+UPPER_PAGE_PAD;
for(int t=0;t<ROOT_STEP;t+=4){
double*c11=root_pending.c+(size_t)i*root_pending.ldc+j+t;
double*c21=c11+(size_t)2048*root_pending.ldc,*c22=c21+2048;
double*m=root_pending.m+root_pending.pos+t;
__m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),m4=_mm256_load_pd(c22);
__m256d m5=_mm256_load_pd(m),m6=_mm256_load_pd(m+PZ),m7=_mm256_load_pd(m+2*PZ);
__m256d u2=_mm256_add_pd(m1,m6),u3=_mm256_add_pd(u2,m7),u4=_mm256_add_pd(u2,m5);
_mm256_store_pd(c11,_mm256_add_pd(m1,m2));
_mm256_store_pd(c21,_mm256_sub_pd(u3,m4));
_mm256_store_pd(c22,_mm256_add_pd(u3,m5));
_mm256_store_pd(m,u4);
}
root_pending.pos+=ROOT_STEP;
if(root_pending.pos==2048*2048)root_pending.active=false;
}
static void root_finish(){while(root_pending.active)root_step();}
static const int lperm[7]={4,0,1,3,2,5,6};
static int g_fuse_sd[7], g_fuse_sgn[7];
static int g_fuse_sd_cur=0, g_fuse_sgn_cur=0;
static const int d1a2_bi[7]={0,0,1,2,3,0,2}, d1a2_bj[7]={3,-1,3,0,-1,1,3}, d1a2_bs[7]={1,0,-1,-1,0,1,1};
static void d1a2_fuse_tab(void){
for(int r=0;r<7;r++){g_fuse_sd[r]=0;g_fuse_sgn[r]=0;}
for(int r=1;r<7;r++){
int q=lperm[r], pq=lperm[r-1];
if(d1a2_bs[pq]!=0) continue;
if(d1a2_bj[q]==d1a2_bi[pq]){ g_fuse_sd[r]=1; g_fuse_sgn[r]=d1a2_bs[q]; }
else if(d1a2_bi[q]==d1a2_bi[pq]){ g_fuse_sd[r]=2; g_fuse_sgn[r]=d1a2_bs[q]; }
}
}
__attribute__((noinline)) static void pk_fuse(int n,Mat b,double *BP,int sgn,int side){
const double *sx=b.x; const size_t sl=b.ld;
for(int k=0;k<n;k+=2)for(int j=0;j<n;j+=8) {
double *p=BP+(size_t)j*(n+1)+k*8;
__m256d q0=_mm256_load_pd(p+0),q1=_mm256_load_pd(p+4);
__m256d q2=_mm256_load_pd(p+8),q3=_mm256_load_pd(p+12);
__m256d s0=_mm256_load_pd(sx+(size_t)(k+0)*sl+j), s1=_mm256_load_pd(sx+(size_t)(k+0)*sl+j+4);
__m256d s2=_mm256_load_pd(sx+(size_t)(k+1)*sl+j), s3=_mm256_load_pd(sx+(size_t)(k+1)*sl+j+4);
if(side==1){
if(sgn>0){q0=_mm256_add_pd(s0,q0);q1=_mm256_add_pd(s1,q1);q2=_mm256_add_pd(s2,q2);q3=_mm256_add_pd(s3,q3);}
else {q0=_mm256_sub_pd(s0,q0);q1=_mm256_sub_pd(s1,q1);q2=_mm256_sub_pd(s2,q2);q3=_mm256_sub_pd(s3,q3);}
} else {
if(sgn>0){q0=_mm256_add_pd(q0,s0);q1=_mm256_add_pd(q1,s1);q2=_mm256_add_pd(q2,s2);q3=_mm256_add_pd(q3,s3);}
else {q0=_mm256_sub_pd(q0,s0);q1=_mm256_sub_pd(q1,s1);q2=_mm256_sub_pd(q2,s2);q3=_mm256_sub_pd(q3,s3);}
}
_mm256_store_pd(p+0,q0);_mm256_store_pd(p+4,q1);
_mm256_store_pd(p+8,q2);_mm256_store_pd(p+12,q3);
}
}
static void d1a2_pk(int n,Mat b,double *BP){
if(g_fuse_sd_cur) pk_fuse(n,b,BP,g_fuse_sgn_cur,g_fuse_sd_cur);
else for(int k=0;k<n;k+=2)for(int j=0;j<n;j+=8) {
double *p=BP+(size_t)j*(n+1)+k*8;
_mm256_store_pd(p+0,ld(b,k+0,j));_mm256_store_pd(p+4,ld(b,k+0,j+4));
_mm256_store_pd(p+8,ld(b,k+1,j));_mm256_store_pd(p+12,ld(b,k+1,j+4));
}
}
template<int TASKS> static void leaf_impl(int leaf_n,Mat a,Mat b,double *c,int ldc) {
constexpr int n=BASE;
d1a2_pk(n,b,bp);
if constexpr(TASKS&1)pending_prefetch();
double *pkb[2]={ap,ap2};
int cur=0;
{ /* pack the first ic-block the ordinary way */
int mc0=n<MC?n:MC;
for(int r=0;r<mc0;r+=MR) {
int rows=mc0-r<MR?mc0-r:MR;
for(int i=0;i<rows;i++)for(int k=0;k<n;k+=4)
_mm256_store_pd(pkb[0]+(size_t)(r+i)*(BASE+24)+k,ld(a,r+i,k));
for(int i=rows;i<MR;i++)for(int k=0;k<n;k+=4)
_mm256_store_pd(pkb[0]+(size_t)(r+i)*(BASE+24)+k,_mm256_setzero_pd());
}
}
for(int ic=0;ic<n;) {
int mc=n-ic<=MC+MR?n-ic:MC;
if(mc==MC+2) {
for(int j=0;j<n;j+=16) {
for(int jj=j;jj<j+16;jj+=NR)for(int r=0;r<MC;r+=MR) {
kernel_plain(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)jj*(n+1),c+(size_t)(ic+r)*ldc+jj,ldc,MR);
if((TASKS&1) && pending.active && (++pending_tick%PERIOD)==0) {pending_step();if(pending.prog)pending_step();pending_prefetch();}
if((TASKS&2) && pend2.active) {pend2_step();pend2_prefetch();}
if((TASKS&4) && root_pending.active && !pending.active) {root_step();root_prefetch();}
}
kernel_tail16(n,pkb[cur]+(size_t)MC*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+MC)*ldc+j,ldc);
}
ic+=mc;
if((TASKS&1) && pending.prog)pending.rowlim=ic;
continue;
}
/* ---- the NEXT ic-block's A panel is carried by this block's kernel calls ---- */
const int ic2=ic+mc;
const int nxt=(ic2<n)?(n-ic2<=MC+MR?n-ic2:MC):0;
const int rows2=((nxt+MR-1)/MR)*MR; /* rows incl. the zero pad rows */
const int nch=rows2*2; /* 2 chunks per row */
const int qd=n/2; /* doubles per chunk */
int ck=0; int skip=(n/NR)*((mc+MR-1)/MR)-nch;
for(int j=0;j<n;j+=NR)for(int r=0;r<mc;r+=MR) {
if(mc-r==2) { /* md9_t2: 2-row tail body (was 6-row body + tmp copy) */
kernel2(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc);
continue;
}
const double *pks=pk_zero,*pky=pk_zero,*pkm=pk_mp;
double *pkd=pk_zero;
int pkc=0;
if(ck<nch && --skip<0) {
const int row=ck/2, off=(ck%2)*qd;
if(row<nxt) {
pks=a.x+(size_t)(ic2+row)*a.ld+(size_t)off;
pky=a.sign?(a.y+(size_t)(ic2+row)*a.ly+(size_t)off):pk_zero;
pkm=a.sign>0?pk_mp:pk_mn;
}
pkd=pkb[cur^1]+(size_t)row*(BASE+24)+(size_t)off;
pkc=8;
ck++;
}
if(!pkc) kernel_plain(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR);
else if(a.sign==0) kernel<2>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
else if(a.sign>0) kernel<1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
else kernel<-1>(n,pkb[cur]+(size_t)r*(BASE+24),bp+(size_t)j*(n+1),c+(size_t)(ic+r)*ldc+j,ldc,mc-r<MR?mc-r:MR,pks,pky,pkd);
if((TASKS&1) && pending.active && (++pending_tick%PERIOD)==0) {
pending_step(); if(pending.prog)pending_step(); pending_prefetch();
}
if((TASKS&2) && pend2.active) { pend2_step(); pend2_prefetch(); }
if((TASKS&4) && root_pending.active && !pending.active) {root_step();root_prefetch();}
}
/* ---- finish whatever the kernel loop did not carry ---- */
for(int cc=ck;cc<nch;cc++) {
const int row=cc/2, off=(cc%2)*qd;
if(row<nxt) {
for(int k=off;k<off+qd;k+=4)
_mm256_store_pd(pkb[cur^1]+(size_t)row*(BASE+24)+k, ld(a,ic2+row,k));
} else {
for(int k=off;k<off+qd;k+=4)
_mm256_store_pd(pkb[cur^1]+(size_t)row*(BASE+24)+k, _mm256_setzero_pd());
}
}
cur^=1;
ic+=mc;
if((TASKS&1) && pending.prog)pending.rowlim=ic;
}
if(n>=NTMIN)_mm_sfence();
}
static void leaf(int n,Mat a,Mat b,double*c,int ldc){int mode=pending.active+2*pend2.active+4*root_pending.active;switch(mode){case 0:leaf_impl<0>(n,a,b,c,ldc);break;case 1:leaf_impl<1>(n,a,b,c,ldc);break;case 2:leaf_impl<2>(n,a,b,c,ldc);break;case 3:leaf_impl<3>(n,a,b,c,ldc);break;case 4:leaf_impl<4>(n,a,b,c,ldc);break;case 5:leaf_impl<5>(n,a,b,c,ldc);break;case 6:leaf_impl<6>(n,a,b,c,ldc);break;case 7:leaf_impl<7>(n,a,b,c,ldc);break;}}
static void sum(int h,double *dst,Mat a,Mat b,int sign) {
for(int i=0;i<h;i++)for(int j=0;j<h;j+=8) {
__m256d v=sign>0?_mm256_add_pd(ld(a,i,j),ld(b,i,j)):_mm256_sub_pd(ld(a,i,j),ld(b,i,j));
__m256d w=sign>0?_mm256_add_pd(ld(a,i,j+4),ld(b,i,j+4)):_mm256_sub_pd(ld(a,i,j+4),ld(b,i,j+4));
asm volatile("" : "+x"(v), "+x"(w) : : "memory");
if(h>=SUMNT) {
_mm256_stream_pd(dst+(size_t)i*h+j,v);_mm256_stream_pd(dst+(size_t)i*h+j+4,w);
} else {
_mm256_store_pd(dst+(size_t)i*h+j,v);_mm256_store_pd(dst+(size_t)i*h+j+4,w);
}
}
if(h>=SUMNT)_mm_sfence();
}
static void combine(int h,double *c,int ldc,const double *m,bool wino=false) {
size_t z=(size_t)h*h+(wino?UPPER_PAGE_PAD:0);
for(int i=0;i<h;i++)for(int j=0;j<h;j+=4) {
size_t t=(size_t)i*h+j;
double *c11=c+(size_t)i*ldc+j,*c12=c11+h,*c21=c+(size_t)(i+h)*ldc+j,*c22=c21+h;
__m256d m1=_mm256_load_pd(c11),m2=_mm256_load_pd(c21),m3=_mm256_load_pd(c12),m4=_mm256_load_pd(c22),
m5=_mm256_load_pd(m+t),m6=_mm256_load_pd(m+z+t),m7=_mm256_load_pd(m+2*z+t);
__m256d r11,r12,r21,r22;
if(wino) {
__m256d u2=_mm256_add_pd(m1,m6),u3=_mm256_add_pd(u2,m7),u4=_mm256_add_pd(u2,m5);
r11=_mm256_add_pd(m1,m2);r12=_mm256_add_pd(u4,m3);r21=_mm256_sub_pd(u3,m4);r22=_mm256_add_pd(u3,m5);
} else {
r11=_mm256_add_pd(_mm256_sub_pd(_mm256_add_pd(m1,m4),m5),m7);r12=_mm256_add_pd(m3,m5);r21=_mm256_add_pd(m2,m4);
r22=_mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(m1,m2),m3),m6);
}
_mm256_store_pd(c11,r11);_mm256_store_pd(c12,r12);_mm256_store_pd(c21,r21);_mm256_store_pd(c22,r22);
}
}
static inline __m256d raw_pre_ld(Mat a,int i,int j){return _mm256_load_pd(a.x+(size_t)i*a.ld+j);}
static void pre(int h,Mat a,double *s1,double *s2,double *s4,double *s3,int ld3,bool right) {
for(int i=0;i<h;i++)for(int j=0;j<h;j+=8) {
__m256d a0=raw_pre_ld(a,i,j),a1=raw_pre_ld(a,i,j+4),b0=raw_pre_ld(a,i,j+h),b1=raw_pre_ld(a,i,j+h+4),
c0=raw_pre_ld(a,i+h,j),c1=raw_pre_ld(a,i+h,j+4),d0=raw_pre_ld(a,i+h,j+h),d1=raw_pre_ld(a,i+h,j+h+4);
__m256d p0,p1,q0,q1,r0,r1,t0,t1;
if(!right) {
p0=_mm256_add_pd(c0,d0);p1=_mm256_add_pd(c1,d1);
q0=_mm256_sub_pd(p0,a0);q1=_mm256_sub_pd(p1,a1);
r0=_mm256_sub_pd(b0,q0);r1=_mm256_sub_pd(b1,q1);
t0=_mm256_sub_pd(a0,c0);t1=_mm256_sub_pd(a1,c1);
}else {
p0=_mm256_sub_pd(b0,a0);p1=_mm256_sub_pd(b1,a1);
q0=_mm256_sub_pd(d0,p0);q1=_mm256_sub_pd(d1,p1);
r0=_mm256_sub_pd(q0,c0);r1=_mm256_sub_pd(q1,c1);
t0=_mm256_sub_pd(d0,b0);t1=_mm256_sub_pd(d1,b1);
}
size_t off=(size_t)i*h+j,off3=(size_t)i*ld3+j;
_mm256_stream_pd(s1+off,p0);_mm256_stream_pd(s1+off+4,p1);
_mm256_stream_pd(s2+off,q0);_mm256_stream_pd(s2+off+4,q1);
_mm256_stream_pd(s4+off,r0);_mm256_stream_pd(s4+off+4,r1);
_mm256_stream_pd(s3+off3,t0);_mm256_stream_pd(s3+off3+4,t1);
}
_mm_sfence();
}
template<int N> static void rec(int rec_n,Mat a,Mat b,double *c,int ldc,double *work,bool root_last=false) {
constexpr int n=N;
if constexpr(n<=BASE){leaf(n,a,b,c,ldc);return;} else {
constexpr int h=n/2;size_t z=(size_t)h*h;
Mat aa[4]={{a.x,a.ld},{a.x+h,a.ld},{a.x+(size_t)h*a.ld,a.ld},{a.x+(size_t)h*a.ld+h,a.ld}};
Mat bb[4]={{b.x,b.ld},{b.x+h,b.ld},{b.x+(size_t)h*b.ld,b.ld},{b.x+(size_t)h*b.ld+h,b.ld}};
double *m=n==2*BASE?work+LEAF_PAD:work;
if(n==2*BASE && pending.active && pending.m==m)m=pending_work+LEAF_PAD;
size_t pz=z+(n==2*BASE?LEAF_PAD:UPPER_PAGE_PAD);
double *p=work+3*z,*q=p+z,*sub=q+z;
double *out[7]={c,c+(size_t)h*ldc,c+h,c+(size_t)h*ldc+h,m,m+pz,m+2*pz};
if(n>BASE*2) {
pending_finish();
double *pa0=m+3*pz,*pa1=pa0+pz,*pa2=pa1+pz,*pb0=pa2+pz,*pb1=pb0+pz,*pb2=pb1+pz;
sub=pb2+pz;
pre(h,a,pa0,pa1,pa2,c,ldc,false);
pre(h,b,pb0,pb1,pb2,c+h,ldc,true);
/* mdA_dp2: child order [1,3,4,5,6,0,2] -- constraint "6 before 0 and 2" is the only
hard one; this order puts the three m-panel writers (4,5,6) as late as possible,
which IS the deferral window for the previous sibling's deferred h=1024 combine. */
rec<h>(h,aa[1],bb[2],out[1],ldc,sub);
rec<h>(h,aa[3],{pb2,h},out[3],ldc,sub);
if(root_last) {
pend2_finish(); // M1 of the root is now fully reconstructed.
root_pending.active=true;
root_prefetch();
}
#if !B2M_NO_BARA
pend2_finish(); /* barrier: next call writes m, the deferred products */
#endif
rec<h>(h,{pa0,h},{pb0,h},out[4],h,sub);
rec<h>(h,{pa1,h},{pb1,h},out[5],h,sub);
rec<h>(h,{c,ldc},{c+h,ldc},out[6],h,sub);
rec<h>(h,aa[0],bb[0],out[0],ldc,sub);
if(n==4096) {
root_pending={c,ldc,m,0,false};
rec<h>(h,{pa2,h},bb[3],out[2],ldc,sub,true);
pending_finish();root_finish();pend2_finish();
return;
}
rec<h>(h,{pa2,h},bb[3],out[2],ldc,sub);
pending_finish();
pend2_finish(); /* finish an older deferred one */
if(root_last)root_finish();
pend2={c,ldc,m,h,(int)pz,0,true}; /* defer THIS combine (wino form);
m-plane stride = the engine's pz */
if(root_last)pend2.add=root_pending.m;
return;
}
const int ai[7]={0,2,0,3,0,2,1},aj[7]={3,3,-1,-1,1,0,3},as[7]={1,1,0,0,1,-1,-1};
const int bi[7]={0,0,1,2,3,0,2},bj[7]={3,-1,3,0,-1,1,3},bs[7]={1,0,-1,-1,0,1,1};
for(int r=0;r<7;r++) {
Mat x=aa[ai[r]],y=bb[bi[r]];
if(h<=BASE) {
/* Independent product order; M7 is last so reconstruction may
advance behind the produced-row frontier of that leaf. */
const int q=lperm[r];
g_fuse_sd_cur=g_fuse_sd[r];g_fuse_sgn_cur=g_fuse_sgn[r];
Mat x=aa[ai[q]],y=bb[bi[q]];
if(as[q]){x.y=aa[aj[q]].x;x.ly=a.ld;x.sign=as[q];}
if(bs[q]){y.y=bb[bj[q]].x;y.ly=b.ld;y.sign=bs[q];}
if(r==6) {
pending_finish();
pending={c,ldc,m,0,true,0,true};
leaf(h,x,y,out[q],q<4?ldc:h);
pending.prog=false;
} else leaf(h,x,y,out[q],q<4?ldc:h);
} else {
if(as[r]){sum(h,p,x,aa[aj[r]],as[r]);x={p,h};}
if(bs[r]){sum(h,q,y,bb[bj[r]],bs[r]);y={q,h};}
rec<h>(h,x,y,out[r],r<4?ldc:h,sub);
}
}
if(h==BASE) {
/* Keep any final portion deferred into the next sibling. */
} else {
pending_finish();
combine(h,c,ldc,m);
}
}
}
void matrix_multiply(int n,const double *a,const double *b,double *c) {
d1a2_fuse_tab();
madvise((void*)workspace,sizeof(workspace),MADV_HUGEPAGE);
madvise((void*)pending_work,sizeof(pending_work),MADV_HUGEPAGE);
rec<4096>(n,{a,n},{b,n},c,n,workspace);
pending_finish();
pend2_finish();
_mm_sfence();
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 1.878 s | 502 MB + 336 KB | Accepted | Score: 100 | 显示更多 |