Imported Upstream version 3.0
[debian/gnuradio] / gnuradio-core / src / lib / filter / fcomplex_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_complex_blocks is != 0
25 #       
26 #
27 #  fcomplex_dotprod_generic (const float *input,
28 #                         const float *taps, unsigned n_2_complex_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];
38 #      sum1 += input[0] * taps[1];
39 #      sum2 += input[1] * taps[2];
40 #      sum3 += input[1] * taps[3];
41 #  
42 #      input += 2;
43 #      taps += 4;
44 #  
45 #    } while (--n_2_complex_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   "fcomplex_dotprod_sse64.S"
59         .version        "01.01"
60 .text
61         .p2align 4
62 .globl GLOB_SYMB(fcomplex_dotprod_sse)
63         DEF_FUNC_HEAD(fcomplex_dotprod_sse)
64 GLOB_SYMB(fcomplex_dotprod_sse):
65
66         # intput: rdi, taps: rsi, n_2_ccomplex_blocks: rdx, result: rcx
67
68         mov     %rdx, %rax
69         
70         # xmm0 xmm1 xmm2 xmm3 are used to hold taps and the result of mults
71         # xmm4 xmm5 xmm6 xmm7 are used to hold the accumulated results
72
73         xorps   %xmm4, %xmm4            # zero two accumulators
74         xorps   %xmm5, %xmm5            # xmm5 holds zero for use below
75
76         # first handle any non-zero remainder of (n_2_complex_blocks % 4)
77
78         and     $0x3, %rax
79         jmp     .L1_test
80
81         .p2align 4
82 .loop1: 
83
84         movlps  0(%rdi), %xmm0
85         shufps  $0x50, %xmm0, %xmm0     # b01010000
86
87         mulps   (%rsi), %xmm0
88         add     $0x10, %rsi
89         add     $8, %rdi
90         addps   %xmm0, %xmm4
91 .L1_test:       
92         dec     %rax
93         jge     .loop1
94
95         
96         # set up for primary loop which is unrolled 4 times
97         
98         movaps  %xmm5, %xmm6            # zero remaining accumulators
99         movaps  %xmm5, %xmm7 
100
101         shr     $2, %rdx                # n_2_complex_blocks / 4
102         je      .cleanup                # if zero, take short path
103
104         # finish setup and loop priming
105
106         movlps  0(%rdi), %xmm0
107
108         movaps  %xmm5, %xmm2
109         movaps  %xmm5, %xmm3
110
111         movlps  8(%rdi), %xmm1
112         shufps  $0x50, %xmm0, %xmm0
113
114         shufps  $0x50, %xmm1, %xmm1
115
116         # we know rdx is not zero, we checked above,
117         # hence enter loop at top
118
119         .p2align 4
120 .loop2:
121         addps   %xmm2, %xmm6
122         movlps  0x10(%rdi), %xmm2
123
124         addps   %xmm3, %xmm7
125
126         mulps   (%rsi), %xmm0
127
128         movlps  0x18(%rdi), %xmm3
129         shufps  $0x50, %xmm2, %xmm2
130
131         mulps   0x10(%rsi), %xmm1
132
133         shufps  $0x50, %xmm3, %xmm3
134
135         addps   %xmm0, %xmm4
136         movlps  0x20(%rdi), %xmm0
137
138         addps   %xmm1, %xmm5
139
140         mulps   0x20(%rsi), %xmm2
141         
142         movlps  0x28(%rdi), %xmm1
143         shufps  $0x50, %xmm0, %xmm0
144
145         mulps   0x30(%rsi), %xmm3
146
147         shufps  $0x50, %xmm1, %xmm1
148
149         add     $0x40, %rsi
150         add     $0x20, %rdi
151         dec     %rdx
152         jne     .loop2
153
154         # OK, now we've done with all the multiplies, but
155         # we still need to handle the unaccumulated
156         # products in xmm2 and xmm3
157
158         addps   %xmm2, %xmm6
159         addps   %xmm3, %xmm7
160
161         # now we want to add all accumulators into xmm4
162
163         addps   %xmm5, %xmm4
164         addps   %xmm6, %xmm7
165         addps   %xmm7, %xmm4
166
167         
168         # At this point, xmm4 contains 2x2 partial sums.  We need
169         # to compute a "horizontal complex add" across xmm4.  
170         
171 .cleanup:                               # xmm4 = r1 i2 r3 i4
172         movhlps %xmm4, %xmm0            # xmm0 = ?? ?? r1 r2
173         addps   %xmm4, %xmm0            # xmm0 = ?? ?? r1+r3 i2+i4
174         movlps  %xmm0, (%rcx)           # store low 2x32 bits (complex) to memory
175
176         retq
177
178 FUNC_TAIL(fcomplex_dotprod_sse)
179         .ident  "Hand coded x86_64 SSE assembly"