GNU Radio Manual and C++ API Reference  3.7.6.1
The Free & Open Software Radio Ecosystem
volk_32fc_x2_dot_prod_32fc.h
Go to the documentation of this file.
1 /* -*- c++ -*- */
2 /*
3  * Copyright 2014 Free Software Foundation, Inc.
4  *
5  * This file is part of GNU Radio
6  *
7  * GNU Radio is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 3, or (at your option)
10  * any later version.
11  *
12  * GNU Radio is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with GNU Radio; see the file COPYING. If not, write to
19  * the Free Software Foundation, Inc., 51 Franklin Street,
20  * Boston, MA 02110-1301, USA.
21  */
22 
23 #ifndef INCLUDED_volk_32fc_x2_dot_prod_32fc_u_H
24 #define INCLUDED_volk_32fc_x2_dot_prod_32fc_u_H
25 
26 #include <volk/volk_common.h>
27 #include <volk/volk_complex.h>
28 #include <stdio.h>
29 #include <string.h>
30 
31 
32 #ifdef LV_HAVE_GENERIC
33 
34 
35 static inline void volk_32fc_x2_dot_prod_32fc_generic(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
36 
37  float * res = (float*) result;
38  float * in = (float*) input;
39  float * tp = (float*) taps;
40  unsigned int n_2_ccomplex_blocks = num_points/2;
41  unsigned int isodd = num_points & 1;
42 
43  float sum0[2] = {0,0};
44  float sum1[2] = {0,0};
45  unsigned int i = 0;
46 
47  for(i = 0; i < n_2_ccomplex_blocks; ++i) {
48  sum0[0] += in[0] * tp[0] - in[1] * tp[1];
49  sum0[1] += in[0] * tp[1] + in[1] * tp[0];
50  sum1[0] += in[2] * tp[2] - in[3] * tp[3];
51  sum1[1] += in[2] * tp[3] + in[3] * tp[2];
52 
53  in += 4;
54  tp += 4;
55  }
56 
57  res[0] = sum0[0] + sum1[0];
58  res[1] = sum0[1] + sum1[1];
59 
60  // Cleanup if we had an odd number of points
61  for(i = 0; i < isodd; ++i) {
62  *result += input[num_points - 1] * taps[num_points - 1];
63  }
64 }
65 
66 #endif /*LV_HAVE_GENERIC*/
67 
68 
69 
70 #if LV_HAVE_SSE && LV_HAVE_64
71 
72 static inline void volk_32fc_x2_dot_prod_32fc_u_sse_64(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
73 
74  const unsigned int num_bytes = num_points*8;
75  unsigned int isodd = num_points & 1;
76 
77  asm
78  (
79  "# ccomplex_dotprod_generic (float* result, const float *input,\n\t"
80  "# const float *taps, unsigned num_bytes)\n\t"
81  "# float sum0 = 0;\n\t"
82  "# float sum1 = 0;\n\t"
83  "# float sum2 = 0;\n\t"
84  "# float sum3 = 0;\n\t"
85  "# do {\n\t"
86  "# sum0 += input[0] * taps[0] - input[1] * taps[1];\n\t"
87  "# sum1 += input[0] * taps[1] + input[1] * taps[0];\n\t"
88  "# sum2 += input[2] * taps[2] - input[3] * taps[3];\n\t"
89  "# sum3 += input[2] * taps[3] + input[3] * taps[2];\n\t"
90  "# input += 4;\n\t"
91  "# taps += 4; \n\t"
92  "# } while (--n_2_ccomplex_blocks != 0);\n\t"
93  "# result[0] = sum0 + sum2;\n\t"
94  "# result[1] = sum1 + sum3;\n\t"
95  "# TODO: prefetch and better scheduling\n\t"
96  " xor %%r9, %%r9\n\t"
97  " xor %%r10, %%r10\n\t"
98  " movq %%rcx, %%rax\n\t"
99  " movq %%rcx, %%r8\n\t"
100  " movq %[rsi], %%r9\n\t"
101  " movq %[rdx], %%r10\n\t"
102  " xorps %%xmm6, %%xmm6 # zero accumulators\n\t"
103  " movups 0(%%r9), %%xmm0\n\t"
104  " xorps %%xmm7, %%xmm7 # zero accumulators\n\t"
105  " movups 0(%%r10), %%xmm2\n\t"
106  " shr $5, %%rax # rax = n_2_ccomplex_blocks / 2\n\t"
107  " shr $4, %%r8\n\t"
108  " jmp .%=L1_test\n\t"
109  " # 4 taps / loop\n\t"
110  " # something like ?? cycles / loop\n\t"
111  ".%=Loop1: \n\t"
112  "# complex prod: C += A * B, w/ temp Z & Y (or B), xmmPN=$0x8000000080000000\n\t"
113  "# movups (%%r9), %%xmmA\n\t"
114  "# movups (%%r10), %%xmmB\n\t"
115  "# movups %%xmmA, %%xmmZ\n\t"
116  "# shufps $0xb1, %%xmmZ, %%xmmZ # swap internals\n\t"
117  "# mulps %%xmmB, %%xmmA\n\t"
118  "# mulps %%xmmZ, %%xmmB\n\t"
119  "# # SSE replacement for: pfpnacc %%xmmB, %%xmmA\n\t"
120  "# xorps %%xmmPN, %%xmmA\n\t"
121  "# movups %%xmmA, %%xmmZ\n\t"
122  "# unpcklps %%xmmB, %%xmmA\n\t"
123  "# unpckhps %%xmmB, %%xmmZ\n\t"
124  "# movups %%xmmZ, %%xmmY\n\t"
125  "# shufps $0x44, %%xmmA, %%xmmZ # b01000100\n\t"
126  "# shufps $0xee, %%xmmY, %%xmmA # b11101110\n\t"
127  "# addps %%xmmZ, %%xmmA\n\t"
128  "# addps %%xmmA, %%xmmC\n\t"
129  "# A=xmm0, B=xmm2, Z=xmm4\n\t"
130  "# A'=xmm1, B'=xmm3, Z'=xmm5\n\t"
131  " movups 16(%%r9), %%xmm1\n\t"
132  " movups %%xmm0, %%xmm4\n\t"
133  " mulps %%xmm2, %%xmm0\n\t"
134  " shufps $0xb1, %%xmm4, %%xmm4 # swap internals\n\t"
135  " movups 16(%%r10), %%xmm3\n\t"
136  " movups %%xmm1, %%xmm5\n\t"
137  " addps %%xmm0, %%xmm6\n\t"
138  " mulps %%xmm3, %%xmm1\n\t"
139  " shufps $0xb1, %%xmm5, %%xmm5 # swap internals\n\t"
140  " addps %%xmm1, %%xmm6\n\t"
141  " mulps %%xmm4, %%xmm2\n\t"
142  " movups 32(%%r9), %%xmm0\n\t"
143  " addps %%xmm2, %%xmm7\n\t"
144  " mulps %%xmm5, %%xmm3\n\t"
145  " add $32, %%r9\n\t"
146  " movups 32(%%r10), %%xmm2\n\t"
147  " addps %%xmm3, %%xmm7\n\t"
148  " add $32, %%r10\n\t"
149  ".%=L1_test:\n\t"
150  " dec %%rax\n\t"
151  " jge .%=Loop1\n\t"
152  " # We've handled the bulk of multiplies up to here.\n\t"
153  " # Let's sse if original n_2_ccomplex_blocks was odd.\n\t"
154  " # If so, we've got 2 more taps to do.\n\t"
155  " and $1, %%r8\n\t"
156  " je .%=Leven\n\t"
157  " # The count was odd, do 2 more taps.\n\t"
158  " # Note that we've already got mm0/mm2 preloaded\n\t"
159  " # from the main loop.\n\t"
160  " movups %%xmm0, %%xmm4\n\t"
161  " mulps %%xmm2, %%xmm0\n\t"
162  " shufps $0xb1, %%xmm4, %%xmm4 # swap internals\n\t"
163  " addps %%xmm0, %%xmm6\n\t"
164  " mulps %%xmm4, %%xmm2\n\t"
165  " addps %%xmm2, %%xmm7\n\t"
166  ".%=Leven:\n\t"
167  " # neg inversor\n\t"
168  " xorps %%xmm1, %%xmm1\n\t"
169  " mov $0x80000000, %%r9\n\t"
170  " movd %%r9, %%xmm1\n\t"
171  " shufps $0x11, %%xmm1, %%xmm1 # b00010001 # 0 -0 0 -0\n\t"
172  " # pfpnacc\n\t"
173  " xorps %%xmm1, %%xmm6\n\t"
174  " movups %%xmm6, %%xmm2\n\t"
175  " unpcklps %%xmm7, %%xmm6\n\t"
176  " unpckhps %%xmm7, %%xmm2\n\t"
177  " movups %%xmm2, %%xmm3\n\t"
178  " shufps $0x44, %%xmm6, %%xmm2 # b01000100\n\t"
179  " shufps $0xee, %%xmm3, %%xmm6 # b11101110\n\t"
180  " addps %%xmm2, %%xmm6\n\t"
181  " # xmm6 = r1 i2 r3 i4\n\t"
182  " movhlps %%xmm6, %%xmm4 # xmm4 = r3 i4 ?? ??\n\t"
183  " addps %%xmm4, %%xmm6 # xmm6 = r1+r3 i2+i4 ?? ??\n\t"
184  " movlps %%xmm6, (%[rdi]) # store low 2x32 bits (complex) to memory\n\t"
185  :
186  :[rsi] "r" (input), [rdx] "r" (taps), "c" (num_bytes), [rdi] "r" (result)
187  :"rax", "r8", "r9", "r10"
188  );
189 
190 
191  if(isodd) {
192  *result += input[num_points - 1] * taps[num_points - 1];
193  }
194 
195  return;
196 
197 }
198 
199 #endif /* LV_HAVE_SSE && LV_HAVE_64 */
200 
201 
202 
203 
204 #ifdef LV_HAVE_SSE3
205 
206 #include <pmmintrin.h>
207 
208 static inline void volk_32fc_x2_dot_prod_32fc_u_sse3(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
209 
210  lv_32fc_t dotProduct;
211  memset(&dotProduct, 0x0, 2*sizeof(float));
212 
213  unsigned int number = 0;
214  const unsigned int halfPoints = num_points/2;
215  unsigned int isodd = num_points & 1;
216 
217  __m128 x, y, yl, yh, z, tmp1, tmp2, dotProdVal;
218 
219  const lv_32fc_t* a = input;
220  const lv_32fc_t* b = taps;
221 
222  dotProdVal = _mm_setzero_ps();
223 
224  for(;number < halfPoints; number++){
225 
226  x = _mm_loadu_ps((float*)a); // Load the ar + ai, br + bi as ar,ai,br,bi
227  y = _mm_loadu_ps((float*)b); // Load the cr + ci, dr + di as cr,ci,dr,di
228 
229  yl = _mm_moveldup_ps(y); // Load yl with cr,cr,dr,dr
230  yh = _mm_movehdup_ps(y); // Load yh with ci,ci,di,di
231 
232  tmp1 = _mm_mul_ps(x,yl); // tmp1 = ar*cr,ai*cr,br*dr,bi*dr
233 
234  x = _mm_shuffle_ps(x,x,0xB1); // Re-arrange x to be ai,ar,bi,br
235 
236  tmp2 = _mm_mul_ps(x,yh); // tmp2 = ai*ci,ar*ci,bi*di,br*di
237 
238  z = _mm_addsub_ps(tmp1,tmp2); // ar*cr-ai*ci, ai*cr+ar*ci, br*dr-bi*di, bi*dr+br*di
239 
240  dotProdVal = _mm_add_ps(dotProdVal, z); // Add the complex multiplication results together
241 
242  a += 2;
243  b += 2;
244  }
245 
246  __VOLK_ATTR_ALIGNED(16) lv_32fc_t dotProductVector[2];
247 
248  _mm_storeu_ps((float*)dotProductVector,dotProdVal); // Store the results back into the dot product vector
249 
250  dotProduct += ( dotProductVector[0] + dotProductVector[1] );
251 
252  if(isodd) {
253  dotProduct += input[num_points - 1] * taps[num_points - 1];
254  }
255 
256  *result = dotProduct;
257 }
258 
259 #endif /*LV_HAVE_SSE3*/
260 
261 #ifdef LV_HAVE_SSE4_1
262 
263 #include <smmintrin.h>
264 
265 static inline void volk_32fc_x2_dot_prod_32fc_u_sse4_1(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
266 
267  unsigned int i = 0;
268  const unsigned int qtr_points = num_points/4;
269  const unsigned int isodd = num_points & 3;
270 
271  __m128 xmm0, xmm1, xmm2, xmm3, xmm4, xmm5, xmm6, xmm7, real0, real1, im0, im1;
272  float *p_input, *p_taps;
273  __m64 *p_result;
274 
275  p_result = (__m64*)result;
276  p_input = (float*)input;
277  p_taps = (float*)taps;
278 
279  static const __m128i neg = {0x000000000000000080000000};
280 
281  real0 = _mm_setzero_ps();
282  real1 = _mm_setzero_ps();
283  im0 = _mm_setzero_ps();
284  im1 = _mm_setzero_ps();
285 
286  for(; i < qtr_points; ++i) {
287  xmm0 = _mm_loadu_ps(p_input);
288  xmm1 = _mm_loadu_ps(p_taps);
289 
290  p_input += 4;
291  p_taps += 4;
292 
293  xmm2 = _mm_loadu_ps(p_input);
294  xmm3 = _mm_loadu_ps(p_taps);
295 
296  p_input += 4;
297  p_taps += 4;
298 
299  xmm4 = _mm_unpackhi_ps(xmm0, xmm2);
300  xmm5 = _mm_unpackhi_ps(xmm1, xmm3);
301  xmm0 = _mm_unpacklo_ps(xmm0, xmm2);
302  xmm2 = _mm_unpacklo_ps(xmm1, xmm3);
303 
304  //imaginary vector from input
305  xmm1 = _mm_unpackhi_ps(xmm0, xmm4);
306  //real vector from input
307  xmm3 = _mm_unpacklo_ps(xmm0, xmm4);
308  //imaginary vector from taps
309  xmm0 = _mm_unpackhi_ps(xmm2, xmm5);
310  //real vector from taps
311  xmm2 = _mm_unpacklo_ps(xmm2, xmm5);
312 
313  xmm4 = _mm_dp_ps(xmm3, xmm2, 0xf1);
314  xmm5 = _mm_dp_ps(xmm1, xmm0, 0xf1);
315 
316  xmm6 = _mm_dp_ps(xmm3, xmm0, 0xf2);
317  xmm7 = _mm_dp_ps(xmm1, xmm2, 0xf2);
318 
319  real0 = _mm_add_ps(xmm4, real0);
320  real1 = _mm_add_ps(xmm5, real1);
321  im0 = _mm_add_ps(xmm6, im0);
322  im1 = _mm_add_ps(xmm7, im1);
323  }
324 
325  real1 = _mm_xor_ps(real1, bit128_p(&neg)->float_vec);
326 
327  im0 = _mm_add_ps(im0, im1);
328  real0 = _mm_add_ps(real0, real1);
329 
330  im0 = _mm_add_ps(im0, real0);
331 
332  _mm_storel_pi(p_result, im0);
333 
334  for(i = num_points-isodd; i < num_points; i++) {
335  *result += input[i] * taps[i];
336  }
337 }
338 
339 #endif /*LV_HAVE_SSE4_1*/
340 
341 #ifdef LV_HAVE_AVX
342 
343 #include <immintrin.h>
344 
345 static inline void volk_32fc_x2_dot_prod_32fc_u_avx(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
346 
347  unsigned int isodd = num_points & 3;
348  unsigned int i = 0;
349  lv_32fc_t dotProduct;
350  memset(&dotProduct, 0x0, 2*sizeof(float));
351 
352  unsigned int number = 0;
353  const unsigned int quarterPoints = num_points / 4;
354 
355  __m256 x, y, yl, yh, z, tmp1, tmp2, dotProdVal;
356 
357  const lv_32fc_t* a = input;
358  const lv_32fc_t* b = taps;
359 
360  dotProdVal = _mm256_setzero_ps();
361 
362  for(;number < quarterPoints; number++){
363  x = _mm256_loadu_ps((float*)a); // Load a,b,e,f as ar,ai,br,bi,er,ei,fr,fi
364  y = _mm256_loadu_ps((float*)b); // Load c,d,g,h as cr,ci,dr,di,gr,gi,hr,hi
365 
366  yl = _mm256_moveldup_ps(y); // Load yl with cr,cr,dr,dr,gr,gr,hr,hr
367  yh = _mm256_movehdup_ps(y); // Load yh with ci,ci,di,di,gi,gi,hi,hi
368 
369  tmp1 = _mm256_mul_ps(x,yl); // tmp1 = ar*cr,ai*cr,br*dr,bi*dr ...
370 
371  x = _mm256_shuffle_ps(x,x,0xB1); // Re-arrange x to be ai,ar,bi,br,ei,er,fi,fr
372 
373  tmp2 = _mm256_mul_ps(x,yh); // tmp2 = ai*ci,ar*ci,bi*di,br*di ...
374 
375  z = _mm256_addsub_ps(tmp1,tmp2); // ar*cr-ai*ci, ai*cr+ar*ci, br*dr-bi*di, bi*dr+br*di
376 
377  dotProdVal = _mm256_add_ps(dotProdVal, z); // Add the complex multiplication results together
378 
379  a += 4;
380  b += 4;
381  }
382 
383  __VOLK_ATTR_ALIGNED(32) lv_32fc_t dotProductVector[4];
384 
385  _mm256_storeu_ps((float*)dotProductVector,dotProdVal); // Store the results back into the dot product vector
386 
387  dotProduct += ( dotProductVector[0] + dotProductVector[1] + dotProductVector[2] + dotProductVector[3]);
388 
389  for(i = num_points-isodd; i < num_points; i++) {
390  dotProduct += input[i] * taps[i];
391  }
392 
393  *result = dotProduct;
394 }
395 
396 #endif /*LV_HAVE_AVX*/
397 
398 
399 #endif /*INCLUDED_volk_32fc_x2_dot_prod_32fc_u_H*/
400 
401 #ifndef INCLUDED_volk_32fc_x2_dot_prod_32fc_a_H
402 #define INCLUDED_volk_32fc_x2_dot_prod_32fc_a_H
403 
404 #include <volk/volk_common.h>
405 #include <volk/volk_complex.h>
406 #include <stdio.h>
407 #include <string.h>
408 
409 
410 #ifdef LV_HAVE_GENERIC
411 
412 
413 static inline void volk_32fc_x2_dot_prod_32fc_a_generic(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
414 
415  const unsigned int num_bytes = num_points*8;
416 
417  float * res = (float*) result;
418  float * in = (float*) input;
419  float * tp = (float*) taps;
420  unsigned int n_2_ccomplex_blocks = num_bytes >> 4;
421  unsigned int isodd = num_points & 1;
422 
423  float sum0[2] = {0,0};
424  float sum1[2] = {0,0};
425  unsigned int i = 0;
426 
427  for(i = 0; i < n_2_ccomplex_blocks; ++i) {
428  sum0[0] += in[0] * tp[0] - in[1] * tp[1];
429  sum0[1] += in[0] * tp[1] + in[1] * tp[0];
430  sum1[0] += in[2] * tp[2] - in[3] * tp[3];
431  sum1[1] += in[2] * tp[3] + in[3] * tp[2];
432 
433  in += 4;
434  tp += 4;
435  }
436 
437  res[0] = sum0[0] + sum1[0];
438  res[1] = sum0[1] + sum1[1];
439 
440  for(i = 0; i < isodd; ++i) {
441  *result += input[num_points - 1] * taps[num_points - 1];
442  }
443 }
444 
445 #endif /*LV_HAVE_GENERIC*/
446 
447 
448 #if LV_HAVE_SSE && LV_HAVE_64
449 
450 
451 static inline void volk_32fc_x2_dot_prod_32fc_a_sse_64(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
452 
453  const unsigned int num_bytes = num_points*8;
454  unsigned int isodd = num_points & 1;
455 
456  asm
457  (
458  "# ccomplex_dotprod_generic (float* result, const float *input,\n\t"
459  "# const float *taps, unsigned num_bytes)\n\t"
460  "# float sum0 = 0;\n\t"
461  "# float sum1 = 0;\n\t"
462  "# float sum2 = 0;\n\t"
463  "# float sum3 = 0;\n\t"
464  "# do {\n\t"
465  "# sum0 += input[0] * taps[0] - input[1] * taps[1];\n\t"
466  "# sum1 += input[0] * taps[1] + input[1] * taps[0];\n\t"
467  "# sum2 += input[2] * taps[2] - input[3] * taps[3];\n\t"
468  "# sum3 += input[2] * taps[3] + input[3] * taps[2];\n\t"
469  "# input += 4;\n\t"
470  "# taps += 4; \n\t"
471  "# } while (--n_2_ccomplex_blocks != 0);\n\t"
472  "# result[0] = sum0 + sum2;\n\t"
473  "# result[1] = sum1 + sum3;\n\t"
474  "# TODO: prefetch and better scheduling\n\t"
475  " xor %%r9, %%r9\n\t"
476  " xor %%r10, %%r10\n\t"
477  " movq %%rcx, %%rax\n\t"
478  " movq %%rcx, %%r8\n\t"
479  " movq %[rsi], %%r9\n\t"
480  " movq %[rdx], %%r10\n\t"
481  " xorps %%xmm6, %%xmm6 # zero accumulators\n\t"
482  " movaps 0(%%r9), %%xmm0\n\t"
483  " xorps %%xmm7, %%xmm7 # zero accumulators\n\t"
484  " movaps 0(%%r10), %%xmm2\n\t"
485  " shr $5, %%rax # rax = n_2_ccomplex_blocks / 2\n\t"
486  " shr $4, %%r8\n\t"
487  " jmp .%=L1_test\n\t"
488  " # 4 taps / loop\n\t"
489  " # something like ?? cycles / loop\n\t"
490  ".%=Loop1: \n\t"
491  "# complex prod: C += A * B, w/ temp Z & Y (or B), xmmPN=$0x8000000080000000\n\t"
492  "# movaps (%%r9), %%xmmA\n\t"
493  "# movaps (%%r10), %%xmmB\n\t"
494  "# movaps %%xmmA, %%xmmZ\n\t"
495  "# shufps $0xb1, %%xmmZ, %%xmmZ # swap internals\n\t"
496  "# mulps %%xmmB, %%xmmA\n\t"
497  "# mulps %%xmmZ, %%xmmB\n\t"
498  "# # SSE replacement for: pfpnacc %%xmmB, %%xmmA\n\t"
499  "# xorps %%xmmPN, %%xmmA\n\t"
500  "# movaps %%xmmA, %%xmmZ\n\t"
501  "# unpcklps %%xmmB, %%xmmA\n\t"
502  "# unpckhps %%xmmB, %%xmmZ\n\t"
503  "# movaps %%xmmZ, %%xmmY\n\t"
504  "# shufps $0x44, %%xmmA, %%xmmZ # b01000100\n\t"
505  "# shufps $0xee, %%xmmY, %%xmmA # b11101110\n\t"
506  "# addps %%xmmZ, %%xmmA\n\t"
507  "# addps %%xmmA, %%xmmC\n\t"
508  "# A=xmm0, B=xmm2, Z=xmm4\n\t"
509  "# A'=xmm1, B'=xmm3, Z'=xmm5\n\t"
510  " movaps 16(%%r9), %%xmm1\n\t"
511  " movaps %%xmm0, %%xmm4\n\t"
512  " mulps %%xmm2, %%xmm0\n\t"
513  " shufps $0xb1, %%xmm4, %%xmm4 # swap internals\n\t"
514  " movaps 16(%%r10), %%xmm3\n\t"
515  " movaps %%xmm1, %%xmm5\n\t"
516  " addps %%xmm0, %%xmm6\n\t"
517  " mulps %%xmm3, %%xmm1\n\t"
518  " shufps $0xb1, %%xmm5, %%xmm5 # swap internals\n\t"
519  " addps %%xmm1, %%xmm6\n\t"
520  " mulps %%xmm4, %%xmm2\n\t"
521  " movaps 32(%%r9), %%xmm0\n\t"
522  " addps %%xmm2, %%xmm7\n\t"
523  " mulps %%xmm5, %%xmm3\n\t"
524  " add $32, %%r9\n\t"
525  " movaps 32(%%r10), %%xmm2\n\t"
526  " addps %%xmm3, %%xmm7\n\t"
527  " add $32, %%r10\n\t"
528  ".%=L1_test:\n\t"
529  " dec %%rax\n\t"
530  " jge .%=Loop1\n\t"
531  " # We've handled the bulk of multiplies up to here.\n\t"
532  " # Let's sse if original n_2_ccomplex_blocks was odd.\n\t"
533  " # If so, we've got 2 more taps to do.\n\t"
534  " and $1, %%r8\n\t"
535  " je .%=Leven\n\t"
536  " # The count was odd, do 2 more taps.\n\t"
537  " # Note that we've already got mm0/mm2 preloaded\n\t"
538  " # from the main loop.\n\t"
539  " movaps %%xmm0, %%xmm4\n\t"
540  " mulps %%xmm2, %%xmm0\n\t"
541  " shufps $0xb1, %%xmm4, %%xmm4 # swap internals\n\t"
542  " addps %%xmm0, %%xmm6\n\t"
543  " mulps %%xmm4, %%xmm2\n\t"
544  " addps %%xmm2, %%xmm7\n\t"
545  ".%=Leven:\n\t"
546  " # neg inversor\n\t"
547  " xorps %%xmm1, %%xmm1\n\t"
548  " mov $0x80000000, %%r9\n\t"
549  " movd %%r9, %%xmm1\n\t"
550  " shufps $0x11, %%xmm1, %%xmm1 # b00010001 # 0 -0 0 -0\n\t"
551  " # pfpnacc\n\t"
552  " xorps %%xmm1, %%xmm6\n\t"
553  " movaps %%xmm6, %%xmm2\n\t"
554  " unpcklps %%xmm7, %%xmm6\n\t"
555  " unpckhps %%xmm7, %%xmm2\n\t"
556  " movaps %%xmm2, %%xmm3\n\t"
557  " shufps $0x44, %%xmm6, %%xmm2 # b01000100\n\t"
558  " shufps $0xee, %%xmm3, %%xmm6 # b11101110\n\t"
559  " addps %%xmm2, %%xmm6\n\t"
560  " # xmm6 = r1 i2 r3 i4\n\t"
561  " movhlps %%xmm6, %%xmm4 # xmm4 = r3 i4 ?? ??\n\t"
562  " addps %%xmm4, %%xmm6 # xmm6 = r1+r3 i2+i4 ?? ??\n\t"
563  " movlps %%xmm6, (%[rdi]) # store low 2x32 bits (complex) to memory\n\t"
564  :
565  :[rsi] "r" (input), [rdx] "r" (taps), "c" (num_bytes), [rdi] "r" (result)
566  :"rax", "r8", "r9", "r10"
567  );
568 
569 
570  if(isodd) {
571  *result += input[num_points - 1] * taps[num_points - 1];
572  }
573 
574  return;
575 
576 }
577 
578 #endif
579 
580 #if LV_HAVE_SSE && LV_HAVE_32
581 
582 static inline void volk_32fc_x2_dot_prod_32fc_a_sse_32(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
583 
584  volk_32fc_x2_dot_prod_32fc_a_generic(result, input, taps, num_points);
585 
586 #if 0
587  const unsigned int num_bytes = num_points*8;
588  unsigned int isodd = num_points & 1;
589 
590  asm volatile
591  (
592  " #pushl %%ebp\n\t"
593  " #movl %%esp, %%ebp\n\t"
594  " movl 12(%%ebp), %%eax # input\n\t"
595  " movl 16(%%ebp), %%edx # taps\n\t"
596  " movl 20(%%ebp), %%ecx # n_bytes\n\t"
597  " xorps %%xmm6, %%xmm6 # zero accumulators\n\t"
598  " movaps 0(%%eax), %%xmm0\n\t"
599  " xorps %%xmm7, %%xmm7 # zero accumulators\n\t"
600  " movaps 0(%%edx), %%xmm2\n\t"
601  " shrl $5, %%ecx # ecx = n_2_ccomplex_blocks / 2\n\t"
602  " jmp .%=L1_test\n\t"
603  " # 4 taps / loop\n\t"
604  " # something like ?? cycles / loop\n\t"
605  ".%=Loop1: \n\t"
606  "# complex prod: C += A * B, w/ temp Z & Y (or B), xmmPN=$0x8000000080000000\n\t"
607  "# movaps (%%eax), %%xmmA\n\t"
608  "# movaps (%%edx), %%xmmB\n\t"
609  "# movaps %%xmmA, %%xmmZ\n\t"
610  "# shufps $0xb1, %%xmmZ, %%xmmZ # swap internals\n\t"
611  "# mulps %%xmmB, %%xmmA\n\t"
612  "# mulps %%xmmZ, %%xmmB\n\t"
613  "# # SSE replacement for: pfpnacc %%xmmB, %%xmmA\n\t"
614  "# xorps %%xmmPN, %%xmmA\n\t"
615  "# movaps %%xmmA, %%xmmZ\n\t"
616  "# unpcklps %%xmmB, %%xmmA\n\t"
617  "# unpckhps %%xmmB, %%xmmZ\n\t"
618  "# movaps %%xmmZ, %%xmmY\n\t"
619  "# shufps $0x44, %%xmmA, %%xmmZ # b01000100\n\t"
620  "# shufps $0xee, %%xmmY, %%xmmA # b11101110\n\t"
621  "# addps %%xmmZ, %%xmmA\n\t"
622  "# addps %%xmmA, %%xmmC\n\t"
623  "# A=xmm0, B=xmm2, Z=xmm4\n\t"
624  "# A'=xmm1, B'=xmm3, Z'=xmm5\n\t"
625  " movaps 16(%%eax), %%xmm1\n\t"
626  " movaps %%xmm0, %%xmm4\n\t"
627  " mulps %%xmm2, %%xmm0\n\t"
628  " shufps $0xb1, %%xmm4, %%xmm4 # swap internals\n\t"
629  " movaps 16(%%edx), %%xmm3\n\t"
630  " movaps %%xmm1, %%xmm5\n\t"
631  " addps %%xmm0, %%xmm6\n\t"
632  " mulps %%xmm3, %%xmm1\n\t"
633  " shufps $0xb1, %%xmm5, %%xmm5 # swap internals\n\t"
634  " addps %%xmm1, %%xmm6\n\t"
635  " mulps %%xmm4, %%xmm2\n\t"
636  " movaps 32(%%eax), %%xmm0\n\t"
637  " addps %%xmm2, %%xmm7\n\t"
638  " mulps %%xmm5, %%xmm3\n\t"
639  " addl $32, %%eax\n\t"
640  " movaps 32(%%edx), %%xmm2\n\t"
641  " addps %%xmm3, %%xmm7\n\t"
642  " addl $32, %%edx\n\t"
643  ".%=L1_test:\n\t"
644  " decl %%ecx\n\t"
645  " jge .%=Loop1\n\t"
646  " # We've handled the bulk of multiplies up to here.\n\t"
647  " # Let's sse if original n_2_ccomplex_blocks was odd.\n\t"
648  " # If so, we've got 2 more taps to do.\n\t"
649  " movl 20(%%ebp), %%ecx # n_2_ccomplex_blocks\n\t"
650  " shrl $4, %%ecx\n\t"
651  " andl $1, %%ecx\n\t"
652  " je .%=Leven\n\t"
653  " # The count was odd, do 2 more taps.\n\t"
654  " # Note that we've already got mm0/mm2 preloaded\n\t"
655  " # from the main loop.\n\t"
656  " movaps %%xmm0, %%xmm4\n\t"
657  " mulps %%xmm2, %%xmm0\n\t"
658  " shufps $0xb1, %%xmm4, %%xmm4 # swap internals\n\t"
659  " addps %%xmm0, %%xmm6\n\t"
660  " mulps %%xmm4, %%xmm2\n\t"
661  " addps %%xmm2, %%xmm7\n\t"
662  ".%=Leven:\n\t"
663  " # neg inversor\n\t"
664  " movl 8(%%ebp), %%eax \n\t"
665  " xorps %%xmm1, %%xmm1\n\t"
666  " movl $0x80000000, (%%eax)\n\t"
667  " movss (%%eax), %%xmm1\n\t"
668  " shufps $0x11, %%xmm1, %%xmm1 # b00010001 # 0 -0 0 -0\n\t"
669  " # pfpnacc\n\t"
670  " xorps %%xmm1, %%xmm6\n\t"
671  " movaps %%xmm6, %%xmm2\n\t"
672  " unpcklps %%xmm7, %%xmm6\n\t"
673  " unpckhps %%xmm7, %%xmm2\n\t"
674  " movaps %%xmm2, %%xmm3\n\t"
675  " shufps $0x44, %%xmm6, %%xmm2 # b01000100\n\t"
676  " shufps $0xee, %%xmm3, %%xmm6 # b11101110\n\t"
677  " addps %%xmm2, %%xmm6\n\t"
678  " # xmm6 = r1 i2 r3 i4\n\t"
679  " #movl 8(%%ebp), %%eax # @result\n\t"
680  " movhlps %%xmm6, %%xmm4 # xmm4 = r3 i4 ?? ??\n\t"
681  " addps %%xmm4, %%xmm6 # xmm6 = r1+r3 i2+i4 ?? ??\n\t"
682  " movlps %%xmm6, (%%eax) # store low 2x32 bits (complex) to memory\n\t"
683  " #popl %%ebp\n\t"
684  :
685  :
686  : "eax", "ecx", "edx"
687  );
688 
689 
690  int getem = num_bytes % 16;
691 
692  if(isodd) {
693  *result += (input[num_points - 1] * taps[num_points - 1]);
694  }
695 
696  return;
697 #endif
698 }
699 
700 #endif /*LV_HAVE_SSE*/
701 
702 #ifdef LV_HAVE_SSE3
703 
704 #include <pmmintrin.h>
705 
706 static inline void volk_32fc_x2_dot_prod_32fc_a_sse3(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
707 
708  const unsigned int num_bytes = num_points*8;
709  unsigned int isodd = num_points & 1;
710 
711  lv_32fc_t dotProduct;
712  memset(&dotProduct, 0x0, 2*sizeof(float));
713 
714  unsigned int number = 0;
715  const unsigned int halfPoints = num_bytes >> 4;
716 
717  __m128 x, y, yl, yh, z, tmp1, tmp2, dotProdVal;
718 
719  const lv_32fc_t* a = input;
720  const lv_32fc_t* b = taps;
721 
722  dotProdVal = _mm_setzero_ps();
723 
724  for(;number < halfPoints; number++){
725 
726  x = _mm_load_ps((float*)a); // Load the ar + ai, br + bi as ar,ai,br,bi
727  y = _mm_load_ps((float*)b); // Load the cr + ci, dr + di as cr,ci,dr,di
728 
729  yl = _mm_moveldup_ps(y); // Load yl with cr,cr,dr,dr
730  yh = _mm_movehdup_ps(y); // Load yh with ci,ci,di,di
731 
732  tmp1 = _mm_mul_ps(x,yl); // tmp1 = ar*cr,ai*cr,br*dr,bi*dr
733 
734  x = _mm_shuffle_ps(x,x,0xB1); // Re-arrange x to be ai,ar,bi,br
735 
736  tmp2 = _mm_mul_ps(x,yh); // tmp2 = ai*ci,ar*ci,bi*di,br*di
737 
738  z = _mm_addsub_ps(tmp1,tmp2); // ar*cr-ai*ci, ai*cr+ar*ci, br*dr-bi*di, bi*dr+br*di
739 
740  dotProdVal = _mm_add_ps(dotProdVal, z); // Add the complex multiplication results together
741 
742  a += 2;
743  b += 2;
744  }
745 
746  __VOLK_ATTR_ALIGNED(16) lv_32fc_t dotProductVector[2];
747 
748  _mm_store_ps((float*)dotProductVector,dotProdVal); // Store the results back into the dot product vector
749 
750  dotProduct += ( dotProductVector[0] + dotProductVector[1] );
751 
752  if(isodd) {
753  dotProduct += input[num_points - 1] * taps[num_points - 1];
754  }
755 
756  *result = dotProduct;
757 }
758 
759 #endif /*LV_HAVE_SSE3*/
760 
761 
762 #ifdef LV_HAVE_SSE4_1
763 
764 #include <smmintrin.h>
765 
766 static inline void volk_32fc_x2_dot_prod_32fc_a_sse4_1(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
767 
768  unsigned int i = 0;
769  const unsigned int qtr_points = num_points/4;
770  const unsigned int isodd = num_points & 3;
771 
772  __m128 xmm0, xmm1, xmm2, xmm3, xmm4, xmm5, xmm6, xmm7, real0, real1, im0, im1;
773  float *p_input, *p_taps;
774  __m64 *p_result;
775 
776  static const __m128i neg = {0x000000000000000080000000};
777 
778  p_result = (__m64*)result;
779  p_input = (float*)input;
780  p_taps = (float*)taps;
781 
782  real0 = _mm_setzero_ps();
783  real1 = _mm_setzero_ps();
784  im0 = _mm_setzero_ps();
785  im1 = _mm_setzero_ps();
786 
787  for(; i < qtr_points; ++i) {
788  xmm0 = _mm_load_ps(p_input);
789  xmm1 = _mm_load_ps(p_taps);
790 
791  p_input += 4;
792  p_taps += 4;
793 
794  xmm2 = _mm_load_ps(p_input);
795  xmm3 = _mm_load_ps(p_taps);
796 
797  p_input += 4;
798  p_taps += 4;
799 
800  xmm4 = _mm_unpackhi_ps(xmm0, xmm2);
801  xmm5 = _mm_unpackhi_ps(xmm1, xmm3);
802  xmm0 = _mm_unpacklo_ps(xmm0, xmm2);
803  xmm2 = _mm_unpacklo_ps(xmm1, xmm3);
804 
805  //imaginary vector from input
806  xmm1 = _mm_unpackhi_ps(xmm0, xmm4);
807  //real vector from input
808  xmm3 = _mm_unpacklo_ps(xmm0, xmm4);
809  //imaginary vector from taps
810  xmm0 = _mm_unpackhi_ps(xmm2, xmm5);
811  //real vector from taps
812  xmm2 = _mm_unpacklo_ps(xmm2, xmm5);
813 
814  xmm4 = _mm_dp_ps(xmm3, xmm2, 0xf1);
815  xmm5 = _mm_dp_ps(xmm1, xmm0, 0xf1);
816 
817  xmm6 = _mm_dp_ps(xmm3, xmm0, 0xf2);
818  xmm7 = _mm_dp_ps(xmm1, xmm2, 0xf2);
819 
820  real0 = _mm_add_ps(xmm4, real0);
821  real1 = _mm_add_ps(xmm5, real1);
822  im0 = _mm_add_ps(xmm6, im0);
823  im1 = _mm_add_ps(xmm7, im1);
824  }
825 
826  real1 = _mm_xor_ps(real1, bit128_p(&neg)->float_vec);
827 
828  im0 = _mm_add_ps(im0, im1);
829  real0 = _mm_add_ps(real0, real1);
830 
831  im0 = _mm_add_ps(im0, real0);
832 
833  _mm_storel_pi(p_result, im0);
834 
835  for(i = num_points-isodd; i < num_points; i++) {
836  *result += input[i] * taps[i];
837  }
838 }
839 
840 #endif /*LV_HAVE_SSE4_1*/
841 
842 #ifdef LV_HAVE_NEON
843 #include <arm_neon.h>
844 
845 static inline void volk_32fc_x2_dot_prod_32fc_neon(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
846 
847  unsigned int quarter_points = num_points / 4;
848  unsigned int number;
849 
850  lv_32fc_t* a_ptr = (lv_32fc_t*) taps;
851  lv_32fc_t* b_ptr = (lv_32fc_t*) input;
852  // for 2-lane vectors, 1st lane holds the real part,
853  // 2nd lane holds the imaginary part
854  float32x4x2_t a_val, b_val, c_val, accumulator;
855  float32x4x2_t tmp_real, tmp_imag;
856  accumulator.val[0] = vdupq_n_f32(0);
857  accumulator.val[1] = vdupq_n_f32(0);
858 
859  for(number = 0; number < quarter_points; ++number) {
860  a_val = vld2q_f32((float*)a_ptr); // a0r|a1r|a2r|a3r || a0i|a1i|a2i|a3i
861  b_val = vld2q_f32((float*)b_ptr); // b0r|b1r|b2r|b3r || b0i|b1i|b2i|b3i
862  __builtin_prefetch(a_ptr+8);
863  __builtin_prefetch(b_ptr+8);
864 
865  // multiply the real*real and imag*imag to get real result
866  // a0r*b0r|a1r*b1r|a2r*b2r|a3r*b3r
867  tmp_real.val[0] = vmulq_f32(a_val.val[0], b_val.val[0]);
868  // a0i*b0i|a1i*b1i|a2i*b2i|a3i*b3i
869  tmp_real.val[1] = vmulq_f32(a_val.val[1], b_val.val[1]);
870 
871  // Multiply cross terms to get the imaginary result
872  // a0r*b0i|a1r*b1i|a2r*b2i|a3r*b3i
873  tmp_imag.val[0] = vmulq_f32(a_val.val[0], b_val.val[1]);
874  // a0i*b0r|a1i*b1r|a2i*b2r|a3i*b3r
875  tmp_imag.val[1] = vmulq_f32(a_val.val[1], b_val.val[0]);
876 
877  c_val.val[0] = vsubq_f32(tmp_real.val[0], tmp_real.val[1]);
878  c_val.val[1] = vaddq_f32(tmp_imag.val[0], tmp_imag.val[1]);
879 
880  accumulator.val[0] = vaddq_f32(accumulator.val[0], c_val.val[0]);
881  accumulator.val[1] = vaddq_f32(accumulator.val[1], c_val.val[1]);
882 
883  a_ptr += 4;
884  b_ptr += 4;
885  }
886  lv_32fc_t accum_result[4];
887  vst2q_f32((float*)accum_result, accumulator);
888  *result = accum_result[0] + accum_result[1] + accum_result[2] + accum_result[3];
889 
890  // tail case
891  for(number = quarter_points*4; number < num_points; ++number) {
892  *result += (*a_ptr++) * (*b_ptr++);
893  }
894 
895 }
896 #endif /*LV_HAVE_NEON*/
897 
898 #ifdef LV_HAVE_NEON
899 #include <arm_neon.h>
900 static inline void volk_32fc_x2_dot_prod_32fc_neon_opttests(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
901 
902  unsigned int quarter_points = num_points / 4;
903  unsigned int number;
904 
905  lv_32fc_t* a_ptr = (lv_32fc_t*) taps;
906  lv_32fc_t* b_ptr = (lv_32fc_t*) input;
907  // for 2-lane vectors, 1st lane holds the real part,
908  // 2nd lane holds the imaginary part
909  float32x4x2_t a_val, b_val, accumulator;
910  float32x4x2_t tmp_imag;
911  accumulator.val[0] = vdupq_n_f32(0);
912  accumulator.val[1] = vdupq_n_f32(0);
913 
914  for(number = 0; number < quarter_points; ++number) {
915  a_val = vld2q_f32((float*)a_ptr); // a0r|a1r|a2r|a3r || a0i|a1i|a2i|a3i
916  b_val = vld2q_f32((float*)b_ptr); // b0r|b1r|b2r|b3r || b0i|b1i|b2i|b3i
917  __builtin_prefetch(a_ptr+8);
918  __builtin_prefetch(b_ptr+8);
919 
920  // do the first multiply
921  tmp_imag.val[1] = vmulq_f32(a_val.val[1], b_val.val[0]);
922  tmp_imag.val[0] = vmulq_f32(a_val.val[0], b_val.val[0]);
923 
924  // use multiply accumulate/subtract to get result
925  tmp_imag.val[1] = vmlaq_f32(tmp_imag.val[1], a_val.val[0], b_val.val[1]);
926  tmp_imag.val[0] = vmlsq_f32(tmp_imag.val[0], a_val.val[1], b_val.val[1]);
927 
928  accumulator.val[0] = vaddq_f32(accumulator.val[0], tmp_imag.val[0]);
929  accumulator.val[1] = vaddq_f32(accumulator.val[1], tmp_imag.val[1]);
930 
931  // increment pointers
932  a_ptr += 4;
933  b_ptr += 4;
934  }
935  lv_32fc_t accum_result[4];
936  vst2q_f32((float*)accum_result, accumulator);
937  *result = accum_result[0] + accum_result[1] + accum_result[2] + accum_result[3];
938 
939  // tail case
940  for(number = quarter_points*4; number < num_points; ++number) {
941  *result += (*a_ptr++) * (*b_ptr++);
942  }
943 
944 }
945 #endif /*LV_HAVE_NEON*/
946 
947 #ifdef LV_HAVE_NEON
948 static inline void volk_32fc_x2_dot_prod_32fc_neon_optfma(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
949 
950  unsigned int quarter_points = num_points / 4;
951  unsigned int number;
952 
953  lv_32fc_t* a_ptr = (lv_32fc_t*) taps;
954  lv_32fc_t* b_ptr = (lv_32fc_t*) input;
955  // for 2-lane vectors, 1st lane holds the real part,
956  // 2nd lane holds the imaginary part
957  float32x4x2_t a_val, b_val, accumulator1, accumulator2;
958  accumulator1.val[0] = vdupq_n_f32(0);
959  accumulator1.val[1] = vdupq_n_f32(0);
960  accumulator2.val[0] = vdupq_n_f32(0);
961  accumulator2.val[1] = vdupq_n_f32(0);
962 
963  for(number = 0; number < quarter_points; ++number) {
964  a_val = vld2q_f32((float*)a_ptr); // a0r|a1r|a2r|a3r || a0i|a1i|a2i|a3i
965  b_val = vld2q_f32((float*)b_ptr); // b0r|b1r|b2r|b3r || b0i|b1i|b2i|b3i
966  __builtin_prefetch(a_ptr+8);
967  __builtin_prefetch(b_ptr+8);
968 
969  // use 2 accumulators to remove inter-instruction data dependencies
970  accumulator1.val[0] = vmlaq_f32(accumulator1.val[0], a_val.val[0], b_val.val[0]);
971  accumulator1.val[1] = vmlaq_f32(accumulator1.val[1], a_val.val[0], b_val.val[1]);
972  accumulator2.val[0] = vmlsq_f32(accumulator2.val[0], a_val.val[1], b_val.val[1]);
973  accumulator2.val[1] = vmlaq_f32(accumulator2.val[1], a_val.val[1], b_val.val[0]);
974  // increment pointers
975  a_ptr += 4;
976  b_ptr += 4;
977  }
978  accumulator1.val[0] = vaddq_f32(accumulator1.val[0], accumulator2.val[0]);
979  accumulator1.val[1] = vaddq_f32(accumulator1.val[1], accumulator2.val[1]);
980  lv_32fc_t accum_result[4];
981  vst2q_f32((float*)accum_result, accumulator1);
982  *result = accum_result[0] + accum_result[1] + accum_result[2] + accum_result[3];
983 
984  // tail case
985  for(number = quarter_points*4; number < num_points; ++number) {
986  *result += (*a_ptr++) * (*b_ptr++);
987  }
988 
989 }
990 #endif /*LV_HAVE_NEON*/
991 
992 #ifdef LV_HAVE_NEON
993 static inline void volk_32fc_x2_dot_prod_32fc_neon_optfmaunroll(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
994 // NOTE: GCC does a poor job with this kernel, but the euivalent ASM code is very fast
995 
996  unsigned int quarter_points = num_points / 8;
997  unsigned int number;
998 
999  lv_32fc_t* a_ptr = (lv_32fc_t*) taps;
1000  lv_32fc_t* b_ptr = (lv_32fc_t*) input;
1001  // for 2-lane vectors, 1st lane holds the real part,
1002  // 2nd lane holds the imaginary part
1003  float32x4x4_t a_val, b_val, accumulator1, accumulator2;
1004  float32x4x2_t reduced_accumulator;
1005  accumulator1.val[0] = vdupq_n_f32(0);
1006  accumulator1.val[1] = vdupq_n_f32(0);
1007  accumulator1.val[2] = vdupq_n_f32(0);
1008  accumulator1.val[3] = vdupq_n_f32(0);
1009  accumulator2.val[0] = vdupq_n_f32(0);
1010  accumulator2.val[1] = vdupq_n_f32(0);
1011  accumulator2.val[2] = vdupq_n_f32(0);
1012  accumulator2.val[3] = vdupq_n_f32(0);
1013 
1014  // 8 input regs, 8 accumulators -> 16/16 neon regs are used
1015  for(number = 0; number < quarter_points; ++number) {
1016  a_val = vld4q_f32((float*)a_ptr); // a0r|a1r|a2r|a3r || a0i|a1i|a2i|a3i
1017  b_val = vld4q_f32((float*)b_ptr); // b0r|b1r|b2r|b3r || b0i|b1i|b2i|b3i
1018  __builtin_prefetch(a_ptr+8);
1019  __builtin_prefetch(b_ptr+8);
1020 
1021  // use 2 accumulators to remove inter-instruction data dependencies
1022  accumulator1.val[0] = vmlaq_f32(accumulator1.val[0], a_val.val[0], b_val.val[0]);
1023  accumulator1.val[1] = vmlaq_f32(accumulator1.val[1], a_val.val[0], b_val.val[1]);
1024 
1025  accumulator1.val[2] = vmlaq_f32(accumulator1.val[2], a_val.val[2], b_val.val[2]);
1026  accumulator1.val[3] = vmlaq_f32(accumulator1.val[3], a_val.val[2], b_val.val[3]);
1027 
1028  accumulator2.val[0] = vmlsq_f32(accumulator2.val[0], a_val.val[1], b_val.val[1]);
1029  accumulator2.val[1] = vmlaq_f32(accumulator2.val[1], a_val.val[1], b_val.val[0]);
1030 
1031  accumulator2.val[2] = vmlsq_f32(accumulator2.val[2], a_val.val[3], b_val.val[3]);
1032  accumulator2.val[3] = vmlaq_f32(accumulator2.val[3], a_val.val[3], b_val.val[2]);
1033  // increment pointers
1034  a_ptr += 8;
1035  b_ptr += 8;
1036  }
1037  // reduce 8 accumulator lanes down to 2 (1 real and 1 imag)
1038  accumulator1.val[0] = vaddq_f32(accumulator1.val[0], accumulator1.val[2]);
1039  accumulator1.val[1] = vaddq_f32(accumulator1.val[1], accumulator1.val[3]);
1040  accumulator2.val[0] = vaddq_f32(accumulator2.val[0], accumulator2.val[2]);
1041  accumulator2.val[1] = vaddq_f32(accumulator2.val[1], accumulator2.val[3]);
1042  reduced_accumulator.val[0] = vaddq_f32(accumulator1.val[0], accumulator2.val[0]);
1043  reduced_accumulator.val[1] = vaddq_f32(accumulator1.val[1], accumulator2.val[1]);
1044  // now reduce accumulators to scalars
1045  lv_32fc_t accum_result[4];
1046  vst2q_f32((float*)accum_result, reduced_accumulator);
1047  *result = accum_result[0] + accum_result[1] + accum_result[2] + accum_result[3];
1048 
1049  // tail case
1050  for(number = quarter_points*8; number < num_points; ++number) {
1051  *result += (*a_ptr++) * (*b_ptr++);
1052  }
1053 
1054 }
1055 #endif /*LV_HAVE_NEON*/
1056 
1057 
1058 #ifdef LV_HAVE_AVX
1059 
1060 #include <immintrin.h>
1061 
1062 static inline void volk_32fc_x2_dot_prod_32fc_a_avx(lv_32fc_t* result, const lv_32fc_t* input, const lv_32fc_t* taps, unsigned int num_points) {
1063 
1064  unsigned int isodd = num_points & 3;
1065  unsigned int i = 0;
1066  lv_32fc_t dotProduct;
1067  memset(&dotProduct, 0x0, 2*sizeof(float));
1068 
1069  unsigned int number = 0;
1070  const unsigned int quarterPoints = num_points / 4;
1071 
1072  __m256 x, y, yl, yh, z, tmp1, tmp2, dotProdVal;
1073 
1074  const lv_32fc_t* a = input;
1075  const lv_32fc_t* b = taps;
1076 
1077  dotProdVal = _mm256_setzero_ps();
1078 
1079  for(;number < quarterPoints; number++){
1080 
1081  x = _mm256_load_ps((float*)a); // Load a,b,e,f as ar,ai,br,bi,er,ei,fr,fi
1082  y = _mm256_load_ps((float*)b); // Load c,d,g,h as cr,ci,dr,di,gr,gi,hr,hi
1083 
1084  yl = _mm256_moveldup_ps(y); // Load yl with cr,cr,dr,dr,gr,gr,hr,hr
1085  yh = _mm256_movehdup_ps(y); // Load yh with ci,ci,di,di,gi,gi,hi,hi
1086 
1087  tmp1 = _mm256_mul_ps(x,yl); // tmp1 = ar*cr,ai*cr,br*dr,bi*dr ...
1088 
1089  x = _mm256_shuffle_ps(x,x,0xB1); // Re-arrange x to be ai,ar,bi,br,ei,er,fi,fr
1090 
1091  tmp2 = _mm256_mul_ps(x,yh); // tmp2 = ai*ci,ar*ci,bi*di,br*di ...
1092 
1093  z = _mm256_addsub_ps(tmp1,tmp2); // ar*cr-ai*ci, ai*cr+ar*ci, br*dr-bi*di, bi*dr+br*di
1094 
1095  dotProdVal = _mm256_add_ps(dotProdVal, z); // Add the complex multiplication results together
1096 
1097  a += 4;
1098  b += 4;
1099  }
1100 
1101  __VOLK_ATTR_ALIGNED(32) lv_32fc_t dotProductVector[4];
1102 
1103  _mm256_store_ps((float*)dotProductVector,dotProdVal); // Store the results back into the dot product vector
1104 
1105  dotProduct += ( dotProductVector[0] + dotProductVector[1] + dotProductVector[2] + dotProductVector[3]);
1106 
1107  for(i = num_points-isodd; i < num_points; i++) {
1108  dotProduct += input[i] * taps[i];
1109  }
1110 
1111  *result = dotProduct;
1112 }
1113 
1114 #endif /*LV_HAVE_AVX*/
1115 
1116 #endif /*INCLUDED_volk_32fc_x2_dot_prod_32fc_a_H*/
bit128_p
#define bit128_p(x)
Definition: volk_common.h:94
volk_complex.h
__VOLK_ATTR_ALIGNED
#define __VOLK_ATTR_ALIGNED(x)
Definition: volk_common.h:27
taps
static const float taps[NSTEPS+1][NTAPS]
Definition: interpolator_taps.h:9
lv_32fc_t
float complex lv_32fc_t
Definition: volk_complex.h:56
volk_common.h