poulpy/poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/spqlios/cplx/cplx_fft_avx512.c
poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/spqlios/cplx/cplx_fft_avx512.cRaw
1#include <immintrin.h>
2
3#include "cplx_fft_internal.h"
4#include "cplx_fft_private.h"
5
6typedef double D2MEM[2];
7typedef double D4MEM[4];
8typedef double D8MEM[8];
9
10EXPORT void cplx_fftvec_addmul_avx512(const CPLX_FFTVEC_ADDMUL_PRECOMP* precomp, void* r, const void* a,
11 const void* b) {
12 const uint32_t m = precomp->m;
13 const D8MEM* aa = (D8MEM*)a;
14 const D8MEM* bb = (D8MEM*)b;
15 D8MEM* rr = (D8MEM*)r;
16 const D8MEM* const aend = aa + (m >> 2);
17 do {
18 /*
19BEGIN_TEMPLATE
20const __m512d ari% = _mm512_loadu_pd(aa[%]);
21const __m512d bri% = _mm512_loadu_pd(bb[%]);
22const __m512d rri% = _mm512_loadu_pd(rr[%]);
23const __m512d bir% = _mm512_shuffle_pd(bri%,bri%, 0b01010101);
24const __m512d aii% = _mm512_shuffle_pd(ari%,ari%, 0b11111111);
25const __m512d pro% = _mm512_fmaddsub_pd(aii%,bir%,rri%);
26const __m512d arr% = _mm512_shuffle_pd(ari%,ari%, 0b00000000);
27const __m512d res% = _mm512_fmaddsub_pd(arr%,bri%,pro%);
28_mm512_storeu_pd(rr[%],res%);
29rr += @; // ONCE
30aa += @; // ONCE
31bb += @; // ONCE
32END_TEMPLATE
33 */
34 // BEGIN_INTERLEAVE 2
35 const __m512d ari0 = _mm512_loadu_pd(aa[0]);
36 const __m512d ari1 = _mm512_loadu_pd(aa[1]);
37 const __m512d bri0 = _mm512_loadu_pd(bb[0]);
38 const __m512d bri1 = _mm512_loadu_pd(bb[1]);
39 const __m512d rri0 = _mm512_loadu_pd(rr[0]);
40 const __m512d rri1 = _mm512_loadu_pd(rr[1]);
41 const __m512d bir0 = _mm512_shuffle_pd(bri0, bri0, 0b01010101);
42 const __m512d bir1 = _mm512_shuffle_pd(bri1, bri1, 0b01010101);
43 const __m512d aii0 = _mm512_shuffle_pd(ari0, ari0, 0b11111111);
44 const __m512d aii1 = _mm512_shuffle_pd(ari1, ari1, 0b11111111);
45 const __m512d pro0 = _mm512_fmaddsub_pd(aii0, bir0, rri0);
46 const __m512d pro1 = _mm512_fmaddsub_pd(aii1, bir1, rri1);
47 const __m512d arr0 = _mm512_shuffle_pd(ari0, ari0, 0b00000000);
48 const __m512d arr1 = _mm512_shuffle_pd(ari1, ari1, 0b00000000);
49 const __m512d res0 = _mm512_fmaddsub_pd(arr0, bri0, pro0);
50 const __m512d res1 = _mm512_fmaddsub_pd(arr1, bri1, pro1);
51 _mm512_storeu_pd(rr[0], res0);
52 _mm512_storeu_pd(rr[1], res1);
53 rr += 2; // ONCE
54 aa += 2; // ONCE
55 bb += 2; // ONCE
56 // END_INTERLEAVE
57 } while (aa < aend);
58}
59
60#if 0
61EXPORT void cplx_fftvec_mul_fma(uint32_t m, void* r, const void* a, const void* b) {
62 const double(*aa)[4] = (double(*)[4])a;
63 const double(*bb)[4] = (double(*)[4])b;
64 double(*rr)[4] = (double(*)[4])r;
65 const double(*const aend)[4] = aa + (m >> 1);
66 do {
67 /*
68BEGIN_TEMPLATE
69const __m256d ari% = _mm256_loadu_pd(aa[%]);
70const __m256d bri% = _mm256_loadu_pd(bb[%]);
71const __m256d bir% = _mm256_shuffle_pd(bri%,bri%, 5); // conj of b
72const __m256d aii% = _mm256_shuffle_pd(ari%,ari%, 15); // im of a
73const __m256d pro% = _mm256_mul_pd(aii%,bir%);
74const __m256d arr% = _mm256_shuffle_pd(ari%,ari%, 0); // rr of a
75const __m256d res% = _mm256_fmaddsub_pd(arr%,bri%,pro%);
76_mm256_storeu_pd(rr[%],res%);
77END_TEMPLATE
78 */
79 // BEGIN_INTERLEAVE 4
80const __m256d ari0 = _mm256_loadu_pd(aa[0]);
81const __m256d ari1 = _mm256_loadu_pd(aa[1]);
82const __m256d ari2 = _mm256_loadu_pd(aa[2]);
83const __m256d ari3 = _mm256_loadu_pd(aa[3]);
84const __m256d bri0 = _mm256_loadu_pd(bb[0]);
85const __m256d bri1 = _mm256_loadu_pd(bb[1]);
86const __m256d bri2 = _mm256_loadu_pd(bb[2]);
87const __m256d bri3 = _mm256_loadu_pd(bb[3]);
88const __m256d bir0 = _mm256_shuffle_pd(bri0,bri0, 5); // conj of b
89const __m256d bir1 = _mm256_shuffle_pd(bri1,bri1, 5); // conj of b
90const __m256d bir2 = _mm256_shuffle_pd(bri2,bri2, 5); // conj of b
91const __m256d bir3 = _mm256_shuffle_pd(bri3,bri3, 5); // conj of b
92const __m256d aii0 = _mm256_shuffle_pd(ari0,ari0, 15); // im of a
93const __m256d aii1 = _mm256_shuffle_pd(ari1,ari1, 15); // im of a
94const __m256d aii2 = _mm256_shuffle_pd(ari2,ari2, 15); // im of a
95const __m256d aii3 = _mm256_shuffle_pd(ari3,ari3, 15); // im of a
96const __m256d pro0 = _mm256_mul_pd(aii0,bir0);
97const __m256d pro1 = _mm256_mul_pd(aii1,bir1);
98const __m256d pro2 = _mm256_mul_pd(aii2,bir2);
99const __m256d pro3 = _mm256_mul_pd(aii3,bir3);
100const __m256d arr0 = _mm256_shuffle_pd(ari0,ari0, 0); // rr of a
101const __m256d arr1 = _mm256_shuffle_pd(ari1,ari1, 0); // rr of a
102const __m256d arr2 = _mm256_shuffle_pd(ari2,ari2, 0); // rr of a
103const __m256d arr3 = _mm256_shuffle_pd(ari3,ari3, 0); // rr of a
104const __m256d res0 = _mm256_fmaddsub_pd(arr0,bri0,pro0);
105const __m256d res1 = _mm256_fmaddsub_pd(arr1,bri1,pro1);
106const __m256d res2 = _mm256_fmaddsub_pd(arr2,bri2,pro2);
107const __m256d res3 = _mm256_fmaddsub_pd(arr3,bri3,pro3);
108_mm256_storeu_pd(rr[0],res0);
109_mm256_storeu_pd(rr[1],res1);
110_mm256_storeu_pd(rr[2],res2);
111_mm256_storeu_pd(rr[3],res3);
112 // END_INTERLEAVE
113 rr += 4;
114 aa += 4;
115 bb += 4;
116 } while (aa < aend);
117}
118
119EXPORT void cplx_fftvec_add_fma(uint32_t m, void* r, const void* a, const void* b) {
120 const double(*aa)[4] = (double(*)[4])a;
121 const double(*bb)[4] = (double(*)[4])b;
122 double(*rr)[4] = (double(*)[4])r;
123 const double(*const aend)[4] = aa + (m >> 1);
124 do {
125 /*
126BEGIN_TEMPLATE
127const __m256d ari% = _mm256_loadu_pd(aa[%]);
128const __m256d bri% = _mm256_loadu_pd(bb[%]);
129const __m256d res% = _mm256_add_pd(ari%,bri%);
130_mm256_storeu_pd(rr[%],res%);
131END_TEMPLATE
132 */
133 // BEGIN_INTERLEAVE 4
134const __m256d ari0 = _mm256_loadu_pd(aa[0]);
135const __m256d ari1 = _mm256_loadu_pd(aa[1]);
136const __m256d ari2 = _mm256_loadu_pd(aa[2]);
137const __m256d ari3 = _mm256_loadu_pd(aa[3]);
138const __m256d bri0 = _mm256_loadu_pd(bb[0]);
139const __m256d bri1 = _mm256_loadu_pd(bb[1]);
140const __m256d bri2 = _mm256_loadu_pd(bb[2]);
141const __m256d bri3 = _mm256_loadu_pd(bb[3]);
142const __m256d res0 = _mm256_add_pd(ari0,bri0);
143const __m256d res1 = _mm256_add_pd(ari1,bri1);
144const __m256d res2 = _mm256_add_pd(ari2,bri2);
145const __m256d res3 = _mm256_add_pd(ari3,bri3);
146_mm256_storeu_pd(rr[0],res0);
147_mm256_storeu_pd(rr[1],res1);
148_mm256_storeu_pd(rr[2],res2);
149_mm256_storeu_pd(rr[3],res3);
150 // END_INTERLEAVE
151 rr += 4;
152 aa += 4;
153 bb += 4;
154 } while (aa < aend);
155}
156
157EXPORT void cplx_fftvec_sub2_to_fma(uint32_t m, void* r, const void* a, const void* b) {
158 const double(*aa)[4] = (double(*)[4])a;
159 const double(*bb)[4] = (double(*)[4])b;
160 double(*rr)[4] = (double(*)[4])r;
161 const double(*const aend)[4] = aa + (m >> 1);
162 do {
163 /*
164BEGIN_TEMPLATE
165const __m256d ari% = _mm256_loadu_pd(aa[%]);
166const __m256d bri% = _mm256_loadu_pd(bb[%]);
167const __m256d sum% = _mm256_add_pd(ari%,bri%);
168const __m256d rri% = _mm256_loadu_pd(rr[%]);
169const __m256d res% = _mm256_sub_pd(rri%,sum%);
170_mm256_storeu_pd(rr[%],res%);
171END_TEMPLATE
172 */
173 // BEGIN_INTERLEAVE 4
174const __m256d ari0 = _mm256_loadu_pd(aa[0]);
175const __m256d ari1 = _mm256_loadu_pd(aa[1]);
176const __m256d ari2 = _mm256_loadu_pd(aa[2]);
177const __m256d ari3 = _mm256_loadu_pd(aa[3]);
178const __m256d bri0 = _mm256_loadu_pd(bb[0]);
179const __m256d bri1 = _mm256_loadu_pd(bb[1]);
180const __m256d bri2 = _mm256_loadu_pd(bb[2]);
181const __m256d bri3 = _mm256_loadu_pd(bb[3]);
182const __m256d sum0 = _mm256_add_pd(ari0,bri0);
183const __m256d sum1 = _mm256_add_pd(ari1,bri1);
184const __m256d sum2 = _mm256_add_pd(ari2,bri2);
185const __m256d sum3 = _mm256_add_pd(ari3,bri3);
186const __m256d rri0 = _mm256_loadu_pd(rr[0]);
187const __m256d rri1 = _mm256_loadu_pd(rr[1]);
188const __m256d rri2 = _mm256_loadu_pd(rr[2]);
189const __m256d rri3 = _mm256_loadu_pd(rr[3]);
190const __m256d res0 = _mm256_sub_pd(rri0,sum0);
191const __m256d res1 = _mm256_sub_pd(rri1,sum1);
192const __m256d res2 = _mm256_sub_pd(rri2,sum2);
193const __m256d res3 = _mm256_sub_pd(rri3,sum3);
194_mm256_storeu_pd(rr[0],res0);
195_mm256_storeu_pd(rr[1],res1);
196_mm256_storeu_pd(rr[2],res2);
197_mm256_storeu_pd(rr[3],res3);
198 // END_INTERLEAVE
199 rr += 4;
200 aa += 4;
201 bb += 4;
202 } while (aa < aend);
203}
204
205EXPORT void cplx_fftvec_copy_fma(uint32_t m, void* r, const void* a) {
206 const double(*aa)[4] = (double(*)[4])a;
207 double(*rr)[4] = (double(*)[4])r;
208 const double(*const aend)[4] = aa + (m >> 1);
209 do {
210 /*
211BEGIN_TEMPLATE
212const __m256d ari% = _mm256_loadu_pd(aa[%]);
213_mm256_storeu_pd(rr[%],ari%);
214END_TEMPLATE
215 */
216 // BEGIN_INTERLEAVE 4
217const __m256d ari0 = _mm256_loadu_pd(aa[0]);
218const __m256d ari1 = _mm256_loadu_pd(aa[1]);
219const __m256d ari2 = _mm256_loadu_pd(aa[2]);
220const __m256d ari3 = _mm256_loadu_pd(aa[3]);
221_mm256_storeu_pd(rr[0],ari0);
222_mm256_storeu_pd(rr[1],ari1);
223_mm256_storeu_pd(rr[2],ari2);
224_mm256_storeu_pd(rr[3],ari3);
225 // END_INTERLEAVE
226 rr += 4;
227 aa += 4;
228 } while (aa < aend);
229}
230
231EXPORT void cplx_fftvec_twiddle_fma(uint32_t m, void* a, void* b, const void* omg) {
232 double(*aa)[4] = (double(*)[4])a;
233 double(*bb)[4] = (double(*)[4])b;
234 const double(*const aend)[4] = aa + (m >> 1);
235 const __m256d om = _mm256_loadu_pd(omg);
236 const __m256d omrr = _mm256_shuffle_pd(om, om, 0);
237 const __m256d omii = _mm256_shuffle_pd(om, om, 15);
238 do {
239 /*
240BEGIN_TEMPLATE
241const __m256d bri% = _mm256_loadu_pd(bb[%]);
242const __m256d bir% = _mm256_shuffle_pd(bri%,bri%,5);
243__m256d p% = _mm256_mul_pd(bir%,omii);
244p% = _mm256_fmaddsub_pd(bri%,omrr,p%);
245const __m256d ari% = _mm256_loadu_pd(aa[%]);
246_mm256_storeu_pd(aa[%],_mm256_add_pd(ari%,p%));
247_mm256_storeu_pd(bb[%],_mm256_sub_pd(ari%,p%));
248END_TEMPLATE
249 */
250 // BEGIN_INTERLEAVE 4
251const __m256d bri0 = _mm256_loadu_pd(bb[0]);
252const __m256d bri1 = _mm256_loadu_pd(bb[1]);
253const __m256d bri2 = _mm256_loadu_pd(bb[2]);
254const __m256d bri3 = _mm256_loadu_pd(bb[3]);
255const __m256d bir0 = _mm256_shuffle_pd(bri0,bri0,5);
256const __m256d bir1 = _mm256_shuffle_pd(bri1,bri1,5);
257const __m256d bir2 = _mm256_shuffle_pd(bri2,bri2,5);
258const __m256d bir3 = _mm256_shuffle_pd(bri3,bri3,5);
259__m256d p0 = _mm256_mul_pd(bir0,omii);
260__m256d p1 = _mm256_mul_pd(bir1,omii);
261__m256d p2 = _mm256_mul_pd(bir2,omii);
262__m256d p3 = _mm256_mul_pd(bir3,omii);
263p0 = _mm256_fmaddsub_pd(bri0,omrr,p0);
264p1 = _mm256_fmaddsub_pd(bri1,omrr,p1);
265p2 = _mm256_fmaddsub_pd(bri2,omrr,p2);
266p3 = _mm256_fmaddsub_pd(bri3,omrr,p3);
267const __m256d ari0 = _mm256_loadu_pd(aa[0]);
268const __m256d ari1 = _mm256_loadu_pd(aa[1]);
269const __m256d ari2 = _mm256_loadu_pd(aa[2]);
270const __m256d ari3 = _mm256_loadu_pd(aa[3]);
271_mm256_storeu_pd(aa[0],_mm256_add_pd(ari0,p0));
272_mm256_storeu_pd(aa[1],_mm256_add_pd(ari1,p1));
273_mm256_storeu_pd(aa[2],_mm256_add_pd(ari2,p2));
274_mm256_storeu_pd(aa[3],_mm256_add_pd(ari3,p3));
275_mm256_storeu_pd(bb[0],_mm256_sub_pd(ari0,p0));
276_mm256_storeu_pd(bb[1],_mm256_sub_pd(ari1,p1));
277_mm256_storeu_pd(bb[2],_mm256_sub_pd(ari2,p2));
278_mm256_storeu_pd(bb[3],_mm256_sub_pd(ari3,p3));
279 // END_INTERLEAVE
280 bb += 4;
281 aa += 4;
282 } while (aa < aend);
283}
284#endif
285
286EXPORT void cplx_fftvec_twiddle_avx512(const CPLX_FFTVEC_TWIDDLE_PRECOMP* precomp, void* a, void* b, const void* omg) {
287 const uint32_t m = precomp->m;
288 D8MEM* aa = (D8MEM*)a;
289 D8MEM* bb = (D8MEM*)b;
290 D8MEM* const aend = aa + (m >> 2);
291 const __m512d om = _mm512_broadcast_f64x4(_mm256_loadu_pd(omg));
292 const __m512d omrr = _mm512_shuffle_pd(om, om, 0b00000000);
293 const __m512d omii = _mm512_shuffle_pd(om, om, 0b11111111);
294 do {
295 /*
296 BEGIN_TEMPLATE
297 const __m512d bri% = _mm512_loadu_pd(bb[%]);
298 const __m512d bir% = _mm512_shuffle_pd(bri%,bri%,0b10011001);
299 __m512d p% = _mm512_mul_pd(bir%,omii);
300 p% = _mm512_fmaddsub_pd(bri%,omrr,p%);
301 const __m512d ari% = _mm512_loadu_pd(aa[%]);
302 _mm512_storeu_pd(aa[%],_mm512_add_pd(ari%,p%));
303 _mm512_storeu_pd(bb[%],_mm512_sub_pd(ari%,p%));
304 bb += @; // ONCE
305 aa += @; // ONCE
306 END_TEMPLATE
307 */
308 // BEGIN_INTERLEAVE 4
309 const __m512d bri0 = _mm512_loadu_pd(bb[0]);
310 const __m512d bri1 = _mm512_loadu_pd(bb[1]);
311 const __m512d bri2 = _mm512_loadu_pd(bb[2]);
312 const __m512d bri3 = _mm512_loadu_pd(bb[3]);
313 const __m512d bir0 = _mm512_shuffle_pd(bri0, bri0, 0b10011001);
314 const __m512d bir1 = _mm512_shuffle_pd(bri1, bri1, 0b10011001);
315 const __m512d bir2 = _mm512_shuffle_pd(bri2, bri2, 0b10011001);
316 const __m512d bir3 = _mm512_shuffle_pd(bri3, bri3, 0b10011001);
317 __m512d p0 = _mm512_mul_pd(bir0, omii);
318 __m512d p1 = _mm512_mul_pd(bir1, omii);
319 __m512d p2 = _mm512_mul_pd(bir2, omii);
320 __m512d p3 = _mm512_mul_pd(bir3, omii);
321 p0 = _mm512_fmaddsub_pd(bri0, omrr, p0);
322 p1 = _mm512_fmaddsub_pd(bri1, omrr, p1);
323 p2 = _mm512_fmaddsub_pd(bri2, omrr, p2);
324 p3 = _mm512_fmaddsub_pd(bri3, omrr, p3);
325 const __m512d ari0 = _mm512_loadu_pd(aa[0]);
326 const __m512d ari1 = _mm512_loadu_pd(aa[1]);
327 const __m512d ari2 = _mm512_loadu_pd(aa[2]);
328 const __m512d ari3 = _mm512_loadu_pd(aa[3]);
329 _mm512_storeu_pd(aa[0], _mm512_add_pd(ari0, p0));
330 _mm512_storeu_pd(aa[1], _mm512_add_pd(ari1, p1));
331 _mm512_storeu_pd(aa[2], _mm512_add_pd(ari2, p2));
332 _mm512_storeu_pd(aa[3], _mm512_add_pd(ari3, p3));
333 _mm512_storeu_pd(bb[0], _mm512_sub_pd(ari0, p0));
334 _mm512_storeu_pd(bb[1], _mm512_sub_pd(ari1, p1));
335 _mm512_storeu_pd(bb[2], _mm512_sub_pd(ari2, p2));
336 _mm512_storeu_pd(bb[3], _mm512_sub_pd(ari3, p3));
337 bb += 4; // ONCE
338 aa += 4; // ONCE
339 // END_INTERLEAVE
340 } while (aa < aend);
341}
342
343EXPORT void cplx_fftvec_bitwiddle_avx512(const CPLX_FFTVEC_BITWIDDLE_PRECOMP* precomp, void* a, uint64_t slicea,
344 const void* omg) {
345 const uint32_t m = precomp->m;
346 const uint64_t OFFSET = slicea / sizeof(D8MEM);
347 D8MEM* aa = (D8MEM*)a;
348 const D8MEM* aend = aa + (m >> 2);
349 const __m512d om = _mm512_broadcast_f64x4(_mm256_loadu_pd(omg));
350 const __m512d om1rr = _mm512_shuffle_pd(om, om, 0);
351 const __m512d om1ii = _mm512_shuffle_pd(om, om, 15);
352 const __m512d om2rr = _mm512_shuffle_pd(om, om, 0);
353 const __m512d om2ii = _mm512_shuffle_pd(om, om, 0);
354 const __m512d om3rr = _mm512_shuffle_pd(om, om, 15);
355 const __m512d om3ii = _mm512_shuffle_pd(om, om, 15);
356 do {
357 /*
358 BEGIN_TEMPLATE
359 __m512d ari% = _mm512_loadu_pd(aa[%]);
360 __m512d bri% = _mm512_loadu_pd((aa+OFFSET)[%]);
361 __m512d cri% = _mm512_loadu_pd((aa+2*OFFSET)[%]);
362 __m512d dri% = _mm512_loadu_pd((aa+3*OFFSET)[%]);
363 __m512d pa% = _mm512_shuffle_pd(cri%,cri%,5);
364 __m512d pb% = _mm512_shuffle_pd(dri%,dri%,5);
365 pa% = _mm512_mul_pd(pa%,om1ii);
366 pb% = _mm512_mul_pd(pb%,om1ii);
367 pa% = _mm512_fmaddsub_pd(cri%,om1rr,pa%);
368 pb% = _mm512_fmaddsub_pd(dri%,om1rr,pb%);
369 cri% = _mm512_sub_pd(ari%,pa%);
370 dri% = _mm512_sub_pd(bri%,pb%);
371 ari% = _mm512_add_pd(ari%,pa%);
372 bri% = _mm512_add_pd(bri%,pb%);
373 pa% = _mm512_shuffle_pd(bri%,bri%,5);
374 pb% = _mm512_shuffle_pd(dri%,dri%,5);
375 pa% = _mm512_mul_pd(pa%,om2ii);
376 pb% = _mm512_mul_pd(pb%,om3ii);
377 pa% = _mm512_fmaddsub_pd(bri%,om2rr,pa%);
378 pb% = _mm512_fmaddsub_pd(dri%,om3rr,pb%);
379 bri% = _mm512_sub_pd(ari%,pa%);
380 dri% = _mm512_sub_pd(cri%,pb%);
381 ari% = _mm512_add_pd(ari%,pa%);
382 cri% = _mm512_add_pd(cri%,pb%);
383 _mm512_storeu_pd(aa[%], ari%);
384 _mm512_storeu_pd((aa+OFFSET)[%],bri%);
385 _mm512_storeu_pd((aa+2*OFFSET)[%],cri%);
386 _mm512_storeu_pd((aa+3*OFFSET)[%],dri%);
387 aa += @; // ONCE
388 END_TEMPLATE
389 */
390 // BEGIN_INTERLEAVE 2
391 __m512d ari0 = _mm512_loadu_pd(aa[0]);
392 __m512d ari1 = _mm512_loadu_pd(aa[1]);
393 __m512d bri0 = _mm512_loadu_pd((aa + OFFSET)[0]);
394 __m512d bri1 = _mm512_loadu_pd((aa + OFFSET)[1]);
395 __m512d cri0 = _mm512_loadu_pd((aa + 2 * OFFSET)[0]);
396 __m512d cri1 = _mm512_loadu_pd((aa + 2 * OFFSET)[1]);
397 __m512d dri0 = _mm512_loadu_pd((aa + 3 * OFFSET)[0]);
398 __m512d dri1 = _mm512_loadu_pd((aa + 3 * OFFSET)[1]);
399 __m512d pa0 = _mm512_shuffle_pd(cri0, cri0, 5);
400 __m512d pa1 = _mm512_shuffle_pd(cri1, cri1, 5);
401 __m512d pb0 = _mm512_shuffle_pd(dri0, dri0, 5);
402 __m512d pb1 = _mm512_shuffle_pd(dri1, dri1, 5);
403 pa0 = _mm512_mul_pd(pa0, om1ii);
404 pa1 = _mm512_mul_pd(pa1, om1ii);
405 pb0 = _mm512_mul_pd(pb0, om1ii);
406 pb1 = _mm512_mul_pd(pb1, om1ii);
407 pa0 = _mm512_fmaddsub_pd(cri0, om1rr, pa0);
408 pa1 = _mm512_fmaddsub_pd(cri1, om1rr, pa1);
409 pb0 = _mm512_fmaddsub_pd(dri0, om1rr, pb0);
410 pb1 = _mm512_fmaddsub_pd(dri1, om1rr, pb1);
411 cri0 = _mm512_sub_pd(ari0, pa0);
412 cri1 = _mm512_sub_pd(ari1, pa1);
413 dri0 = _mm512_sub_pd(bri0, pb0);
414 dri1 = _mm512_sub_pd(bri1, pb1);
415 ari0 = _mm512_add_pd(ari0, pa0);
416 ari1 = _mm512_add_pd(ari1, pa1);
417 bri0 = _mm512_add_pd(bri0, pb0);
418 bri1 = _mm512_add_pd(bri1, pb1);
419 pa0 = _mm512_shuffle_pd(bri0, bri0, 5);
420 pa1 = _mm512_shuffle_pd(bri1, bri1, 5);
421 pb0 = _mm512_shuffle_pd(dri0, dri0, 5);
422 pb1 = _mm512_shuffle_pd(dri1, dri1, 5);
423 pa0 = _mm512_mul_pd(pa0, om2ii);
424 pa1 = _mm512_mul_pd(pa1, om2ii);
425 pb0 = _mm512_mul_pd(pb0, om3ii);
426 pb1 = _mm512_mul_pd(pb1, om3ii);
427 pa0 = _mm512_fmaddsub_pd(bri0, om2rr, pa0);
428 pa1 = _mm512_fmaddsub_pd(bri1, om2rr, pa1);
429 pb0 = _mm512_fmaddsub_pd(dri0, om3rr, pb0);
430 pb1 = _mm512_fmaddsub_pd(dri1, om3rr, pb1);
431 bri0 = _mm512_sub_pd(ari0, pa0);
432 bri1 = _mm512_sub_pd(ari1, pa1);
433 dri0 = _mm512_sub_pd(cri0, pb0);
434 dri1 = _mm512_sub_pd(cri1, pb1);
435 ari0 = _mm512_add_pd(ari0, pa0);
436 ari1 = _mm512_add_pd(ari1, pa1);
437 cri0 = _mm512_add_pd(cri0, pb0);
438 cri1 = _mm512_add_pd(cri1, pb1);
439 _mm512_storeu_pd(aa[0], ari0);
440 _mm512_storeu_pd(aa[1], ari1);
441 _mm512_storeu_pd((aa + OFFSET)[0], bri0);
442 _mm512_storeu_pd((aa + OFFSET)[1], bri1);
443 _mm512_storeu_pd((aa + 2 * OFFSET)[0], cri0);
444 _mm512_storeu_pd((aa + 2 * OFFSET)[1], cri1);
445 _mm512_storeu_pd((aa + 3 * OFFSET)[0], dri0);
446 _mm512_storeu_pd((aa + 3 * OFFSET)[1], dri1);
447 aa += 2; // ONCE
448 // END_INTERLEAVE
449 } while (aa < aend);
450 _mm256_zeroupper();
451}