Imported Upstream version 3.0
[debian/gnuradio] / gnuradio-core / src / lib / filter / ccomplex_dotprod_sse64.S
1 #
2 # Copyright 2002,2005 Free Software Foundation, Inc.
3
4 # This file is part of GNU Radio
5
6 # GNU Radio is free software; you can redistribute it and/or modify
7 # it under the terms of the GNU General Public License as published by
8 # the Free Software Foundation; either version 2, or (at your option)
9 # any later version.
10
11 # GNU Radio is distributed in the hope that it will be useful,
12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 # GNU General Public License for more details.
15
16 # You should have received a copy of the GNU General Public License
17 # along with GNU Radio; see the file COPYING.  If not, write to
18 # the Free Software Foundation, Inc., 51 Franklin Street,
19 # Boston, MA 02110-1301, USA.
20
21
22
23 # input and taps are guarenteed to be 16 byte aligned.
24 # n_2_ccomplex_blocks is != 0
25 #       
26 #
27 #  ccomplex_dotprod_generic (const float *input,
28 #                         const float *taps, unsigned n_2_ccomplex_blocks, float *result)
29 #  {
30 #    float sum0 = 0;
31 #    float sum1 = 0;
32 #    float sum2 = 0;
33 #    float sum3 = 0;
34 #  
35 #    do {
36 #
37 #      sum0 += input[0] * taps[0] - input[1] * taps[1];
38 #      sum1 += input[0] * taps[1] + input[1] * taps[0];
39 #      sum2 += input[2] * taps[2] - input[3] * taps[3];
40 #      sum3 += input[2] * taps[3] + input[3] * taps[2];
41 #  
42 #      input += 4;
43 #      taps += 4;  
44 #  
45 #    } while (--n_2_ccomplex_blocks != 0);
46 #  
47 #  
48 #    result[0] = sum0 + sum2;
49 #    result[1] = sum1 + sum3;
50 #  }
51 #
52
53 # TODO: prefetch and better scheduling
54
55 #include "assembly.h"
56
57
58         .file   "ccomplex_dotprod_sse64.S"
59         .version        "01.01"
60 .text
61         .p2align 4
62 .globl GLOB_SYMB(ccomplex_dotprod_sse)
63         DEF_FUNC_HEAD(ccomplex_dotprod_sse)
64 GLOB_SYMB(ccomplex_dotprod_sse):
65
66         # intput: rdi, taps: rsi, n_2_ccomplex_blocks: rdx, result: rcx
67
68         mov     %rdx, %rax
69
70         xorps   %xmm6, %xmm6            # zero accumulators
71         
72         movaps  0(%rdi), %xmm0
73
74         xorps   %xmm7, %xmm7            # zero accumulators
75
76         movaps  0(%rsi), %xmm2
77
78         shr     $1, %rax                # rax = n_2_ccomplex_blocks / 2
79
80         jmp     .L1_test
81
82         #
83         # 4 taps / loop
84         # something like ?? cycles / loop
85         #
86         
87         .p2align 4
88 .loop1: 
89
90 # complex prod: C += A * B,  w/ temp Z & Y (or B), xmmPN=$0x8000000080000000
91 #
92 #       movaps  (%rdi), %xmmA
93 #       movaps  (%rsi), %xmmB
94 #
95 #       movaps  %xmmA, %xmmZ
96 #       shufps  $0xb1, %xmmZ, %xmmZ     # swap internals
97 #
98 #       mulps   %xmmB, %xmmA
99 #       mulps   %xmmZ, %xmmB
100 #
101 #       # SSE replacement for: pfpnacc %xmmB, %xmmA
102 #       xorps   %xmmPN, %xmmA
103 #       movaps  %xmmA, %xmmZ
104 #       unpcklps %xmmB, %xmmA
105 #       unpckhps %xmmB, %xmmZ
106 #       movaps  %xmmZ, %xmmY
107 #       shufps  $0x44, %xmmA, %xmmZ     # b01000100
108 #       shufps  $0xee, %xmmY, %xmmA     # b11101110
109 #       addps   %xmmZ, %xmmA
110 #
111 #       addps   %xmmA, %xmmC
112
113 # A=xmm0, B=xmm2, Z=xmm4
114 # A'=xmm1, B'=xmm3, Z'=xmm5
115
116         movaps  16(%rdi), %xmm1
117
118         movaps  %xmm0, %xmm4
119         mulps   %xmm2, %xmm0
120
121         shufps  $0xb1, %xmm4, %xmm4     # swap internals
122         movaps  16(%rsi), %xmm3
123         movaps  %xmm1, %xmm5
124         addps   %xmm0, %xmm6
125         mulps   %xmm3, %xmm1
126         shufps  $0xb1, %xmm5, %xmm5     # swap internals
127         addps   %xmm1, %xmm6
128         mulps   %xmm4, %xmm2
129         movaps  32(%rdi), %xmm0
130         addps   %xmm2, %xmm7
131         mulps   %xmm5, %xmm3
132
133         add     $32, %rdi
134
135         movaps  32(%rsi), %xmm2
136         addps   %xmm3, %xmm7
137
138         add     $32, %rsi
139
140
141
142 .L1_test:
143         dec     %rax
144         jge     .loop1
145
146         # We've handled the bulk of multiplies up to here.
147         # Let's sse if original n_2_ccomplex_blocks was odd.
148         # If so, we've got 2 more taps to do.
149         
150         and     $1, %rdx
151         je      .Leven
152         
153         # The count was odd, do 2 more taps.
154         # Note that we've already got mm0/mm2 preloaded
155         # from the main loop.
156
157         movaps  %xmm0, %xmm4
158         mulps   %xmm2, %xmm0
159         shufps  $0xb1, %xmm4, %xmm4     # swap internals
160         addps   %xmm0, %xmm6
161         mulps   %xmm4, %xmm2
162         addps   %xmm2, %xmm7
163
164
165 .Leven:
166         # neg inversor
167         xorps   %xmm1, %xmm1
168         movl    $0x80000000, -8(%rsp)
169         movss   -8(%rsp), %xmm1
170         shufps  $0x11, %xmm1, %xmm1     # b00010001 # 0 -0 0 -0
171
172         # pfpnacc
173         xorps   %xmm1, %xmm6
174
175         movaps  %xmm6, %xmm2
176         unpcklps %xmm7, %xmm6
177         unpckhps %xmm7, %xmm2
178         movaps  %xmm2, %xmm3
179         shufps  $0x44, %xmm6, %xmm2     # b01000100
180         shufps  $0xee, %xmm3, %xmm6     # b11101110
181         addps   %xmm2, %xmm6
182
183                                         # xmm6 = r1 i2 r3 i4
184         movhlps %xmm6, %xmm4            # xmm4 = r3 i4 ?? ??
185         addps   %xmm4, %xmm6            # xmm6 = r1+r3 i2+i4 ?? ??
186         movlps  %xmm6, (%rcx)           # store low 2x32 bits (complex) to memory
187
188         retq
189
190 FUNC_TAIL(ccomplex_dotprod_sse)
191         .ident  "Hand coded x86_64 SSE assembly"