Chris@42
|
1 /*
|
Chris@42
|
2 * Copyright (c) 2003, 2007-14 Matteo Frigo
|
Chris@42
|
3 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
|
Chris@42
|
4 *
|
Chris@42
|
5 * VSX SIMD implementation added 2015 Erik Lindahl.
|
Chris@42
|
6 * Erik Lindahl places his modifications in the public domain.
|
Chris@42
|
7 *
|
Chris@42
|
8 * This program is free software; you can redistribute it and/or modify
|
Chris@42
|
9 * it under the terms of the GNU General Public License as published by
|
Chris@42
|
10 * the Free Software Foundation; either version 2 of the License, or
|
Chris@42
|
11 * (at your option) any later version.
|
Chris@42
|
12 *
|
Chris@42
|
13 * This program is distributed in the hope that it will be useful,
|
Chris@42
|
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
Chris@42
|
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
Chris@42
|
16 * GNU General Public License for more details.
|
Chris@42
|
17 *
|
Chris@42
|
18 * You should have received a copy of the GNU General Public License
|
Chris@42
|
19 * along with this program; if not, write to the Free Software
|
Chris@42
|
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
|
Chris@42
|
21 *
|
Chris@42
|
22 */
|
Chris@42
|
23
|
Chris@42
|
24 #if defined(FFTW_LDOUBLE) || defined(FFTW_QUAD)
|
Chris@42
|
25 # error "VSX only works in single or double precision"
|
Chris@42
|
26 #endif
|
Chris@42
|
27
|
Chris@42
|
28 #ifdef FFTW_SINGLE
|
Chris@42
|
29 # define DS(d,s) s /* single-precision option */
|
Chris@42
|
30 # define SUFF(name) name ## s
|
Chris@42
|
31 #else
|
Chris@42
|
32 # define DS(d,s) d /* double-precision option */
|
Chris@42
|
33 # define SUFF(name) name ## d
|
Chris@42
|
34 #endif
|
Chris@42
|
35
|
Chris@42
|
36 #define SIMD_SUFFIX _vsx /* for renaming */
|
Chris@42
|
37 #define VL DS(1,2) /* SIMD vector length, in term of complex numbers */
|
Chris@42
|
38 #define SIMD_VSTRIDE_OKA(x) DS(1,((x) == 2))
|
Chris@42
|
39 #define SIMD_STRIDE_OKPAIR SIMD_STRIDE_OK
|
Chris@42
|
40
|
Chris@42
|
41 #include <altivec.h>
|
Chris@42
|
42 #include <stdio.h>
|
Chris@42
|
43
|
Chris@42
|
44 typedef DS(vector double,vector float) V;
|
Chris@42
|
45
|
Chris@42
|
46 #define VADD(a,b) vec_add(a,b)
|
Chris@42
|
47 #define VSUB(a,b) vec_sub(a,b)
|
Chris@42
|
48 #define VMUL(a,b) vec_mul(a,b)
|
Chris@42
|
49 #define VXOR(a,b) vec_xor(a,b)
|
Chris@42
|
50 #define UNPCKL(a,b) vec_mergel(a,b)
|
Chris@42
|
51 #define UNPCKH(a,b) vec_mergeh(a,b)
|
Chris@42
|
52 #ifdef FFTW_SINGLE
|
Chris@42
|
53 # define VDUPL(a) ({ const vector unsigned char perm = {0,1,2,3,0,1,2,3,8,9,10,11,8,9,10,11}; vec_perm(a,a,perm); })
|
Chris@42
|
54 # define VDUPH(a) ({ const vector unsigned char perm = {4,5,6,7,4,5,6,7,12,13,14,15,12,13,14,15}; vec_perm(a,a,perm); })
|
Chris@42
|
55 #else
|
Chris@42
|
56 # define VDUPL(a) ({ const vector unsigned char perm = {0,1,2,3,4,5,6,7,0,1,2,3,4,5,6,7}; vec_perm(a,a,perm); })
|
Chris@42
|
57 # define VDUPH(a) ({ const vector unsigned char perm = {8,9,10,11,12,13,14,15,8,9,10,11,12,13,14,15}; vec_perm(a,a,perm); })
|
Chris@42
|
58 #endif
|
Chris@42
|
59
|
Chris@42
|
60 static inline V LDK(R f) { return vec_splats(f); }
|
Chris@42
|
61
|
Chris@42
|
62 #define DVK(var, val) const R var = K(val)
|
Chris@42
|
63
|
Chris@42
|
64 static inline V VCONJ(V x)
|
Chris@42
|
65 {
|
Chris@42
|
66 const V pmpm = vec_mergel(vec_splats((R)0.0),-(vec_splats((R)0.0)));
|
Chris@42
|
67 return vec_xor(x, pmpm);
|
Chris@42
|
68 }
|
Chris@42
|
69
|
Chris@42
|
70 static inline V LDA(const R *x, INT ivs, const R *aligned_like)
|
Chris@42
|
71 {
|
Chris@42
|
72 #ifdef __ibmxl__
|
Chris@42
|
73 return vec_xl(0,(DS(double,float) *)x);
|
Chris@42
|
74 #else
|
Chris@42
|
75 return (*(const V *)(x));
|
Chris@42
|
76 #endif
|
Chris@42
|
77 }
|
Chris@42
|
78
|
Chris@42
|
79 static inline void STA(R *x, V v, INT ovs, const R *aligned_like)
|
Chris@42
|
80 {
|
Chris@42
|
81 #ifdef __ibmxl__
|
Chris@42
|
82 vec_xst(v,0,x);
|
Chris@42
|
83 #else
|
Chris@42
|
84 *(V *)x = v;
|
Chris@42
|
85 #endif
|
Chris@42
|
86 }
|
Chris@42
|
87
|
Chris@42
|
88 static inline V FLIP_RI(V x)
|
Chris@42
|
89 {
|
Chris@42
|
90 #ifdef FFTW_SINGLE
|
Chris@42
|
91 const vector unsigned char perm = { 4,5,6,7,0,1,2,3,12,13,14,15,8,9,10,11 };
|
Chris@42
|
92 #else
|
Chris@42
|
93 const vector unsigned char perm = { 8,9,10,11,12,13,14,15,0,1,2,3,4,5,6,7 };
|
Chris@42
|
94 #endif
|
Chris@42
|
95 return vec_perm(x,x,perm);
|
Chris@42
|
96 }
|
Chris@42
|
97
|
Chris@42
|
98 #ifdef FFTW_SINGLE
|
Chris@42
|
99
|
Chris@42
|
100 static inline V LD(const R *x, INT ivs, const R *aligned_like)
|
Chris@42
|
101 {
|
Chris@42
|
102 const vector unsigned char perm = {0,1,2,3,4,5,6,7,16,17,18,19,20,21,22,23};
|
Chris@42
|
103
|
Chris@42
|
104 return vec_perm((vector float)vec_splats(*(double *)(x)),
|
Chris@42
|
105 (vector float)vec_splats(*(double *)(x+ivs)),perm);
|
Chris@42
|
106 }
|
Chris@42
|
107
|
Chris@42
|
108 static inline void ST(R *x, V v, INT ovs, const R *aligned_like)
|
Chris@42
|
109 {
|
Chris@42
|
110 *(double *)(x+ovs) = vec_extract( (vector double)v, 1 );
|
Chris@42
|
111 *(double *)x = vec_extract( (vector double)v, 0 );
|
Chris@42
|
112 }
|
Chris@42
|
113 #else
|
Chris@42
|
114 /* DOUBLE */
|
Chris@42
|
115
|
Chris@42
|
116 # define LD LDA
|
Chris@42
|
117 # define ST STA
|
Chris@42
|
118
|
Chris@42
|
119 #endif
|
Chris@42
|
120
|
Chris@42
|
121 #define STM2 DS(STA,ST)
|
Chris@42
|
122 #define STN2(x, v0, v1, ovs) /* nop */
|
Chris@42
|
123
|
Chris@42
|
124 #ifdef FFTW_SINGLE
|
Chris@42
|
125
|
Chris@42
|
126 # define STM4(x, v, ovs, aligned_like) /* no-op */
|
Chris@42
|
127 static inline void STN4(R *x, V v0, V v1, V v2, V v3, int ovs)
|
Chris@42
|
128 {
|
Chris@42
|
129 V xxx0, xxx1, xxx2, xxx3;
|
Chris@42
|
130 xxx0 = vec_mergeh(v0,v1);
|
Chris@42
|
131 xxx1 = vec_mergel(v0,v1);
|
Chris@42
|
132 xxx2 = vec_mergeh(v2,v3);
|
Chris@42
|
133 xxx3 = vec_mergel(v2,v3);
|
Chris@42
|
134 *(double *)x = vec_extract( (vector double)xxx0, 0 );
|
Chris@42
|
135 *(double *)(x+ovs) = vec_extract( (vector double)xxx0, 1 );
|
Chris@42
|
136 *(double *)(x+2*ovs) = vec_extract( (vector double)xxx1, 0 );
|
Chris@42
|
137 *(double *)(x+3*ovs) = vec_extract( (vector double)xxx1, 1 );
|
Chris@42
|
138 *(double *)(x+2) = vec_extract( (vector double)xxx2, 0 );
|
Chris@42
|
139 *(double *)(x+ovs+2) = vec_extract( (vector double)xxx2, 1 );
|
Chris@42
|
140 *(double *)(x+2*ovs+2) = vec_extract( (vector double)xxx3, 0 );
|
Chris@42
|
141 *(double *)(x+3*ovs+2) = vec_extract( (vector double)xxx3, 1 );
|
Chris@42
|
142 }
|
Chris@42
|
143 #else /* !FFTW_SINGLE */
|
Chris@42
|
144
|
Chris@42
|
145 static inline void STM4(R *x, V v, INT ovs, const R *aligned_like)
|
Chris@42
|
146 {
|
Chris@42
|
147 (void)aligned_like; /* UNUSED */
|
Chris@42
|
148 x[0] = vec_extract(v,0);
|
Chris@42
|
149 x[ovs] = vec_extract(v,1);
|
Chris@42
|
150 }
|
Chris@42
|
151 # define STN4(x, v0, v1, v2, v3, ovs) /* nothing */
|
Chris@42
|
152 #endif
|
Chris@42
|
153
|
Chris@42
|
154 static inline V VBYI(V x)
|
Chris@42
|
155 {
|
Chris@42
|
156 /* Complicated low-level stuff. vpermxor is really a cryptographic instruction that is only
|
Chris@42
|
157 * available in the low-level inteface both for GCC and XLC. However, on little-endian
|
Chris@42
|
158 * platforms there is also the complicated swapping going on. XLC does this here too, but
|
Chris@42
|
159 * not GCC, so we need different permute constants.
|
Chris@42
|
160 */
|
Chris@42
|
161 #if defined(__POWER8_VECTOR__) && defined(__GNUC__) && defined(__LITTLE_ENDIAN__)
|
Chris@42
|
162 # ifdef FFTW_SINGLE
|
Chris@42
|
163 const vector unsigned char perm = { 0xbb, 0xaa, 0x99, 0x88, 0xff, 0xee, 0xdd, 0xcc, 0x33, 0x22, 0x11, 0x00, 0x77, 0x66, 0x55, 0x44 };
|
Chris@42
|
164 # else
|
Chris@42
|
165 const vector unsigned char perm = { 0x77, 0x66, 0x55, 0x44, 0x33, 0x22, 0x11, 0x00, 0xff, 0xee, 0xdd, 0xcc, 0xbb, 0xaa, 0x99, 0x88 };
|
Chris@42
|
166 # endif
|
Chris@42
|
167 const V pmpm = vec_mergel(vec_splats((R)0.0),-(vec_splats((R)0.0)));
|
Chris@42
|
168 return (V)__builtin_crypto_vpermxor((vector unsigned char)x,(vector unsigned char)pmpm,perm);
|
Chris@42
|
169 #elif defined(__POWER8_VECTOR__) && (defined(__ibmxl__) || (defined(__GNUC__) && !defined(__LITTLE_ENDIAN__)))
|
Chris@42
|
170 # ifdef FFTW_SINGLE
|
Chris@42
|
171 const vector unsigned char perm = { 0x44, 0x55, 0x66, 0x77, 0x00, 0x11, 0x22, 0x33, 0xCC, 0xDD, 0xEE, 0xFF, 0x88, 0x99, 0xAA, 0xBB };
|
Chris@42
|
172 # else
|
Chris@42
|
173 const vector unsigned char perm = { 0x88, 0x99, 0xAA, 0xBB, 0xCC, 0xDD, 0xEE, 0xFF, 0x00, 0x11, 0x22, 0x33, 0x44, 0x55, 0x66, 0x77 };
|
Chris@42
|
174 # endif
|
Chris@42
|
175 const V pmpm = vec_mergel(vec_splats((R)0.0),-(vec_splats((R)0.0)));
|
Chris@42
|
176 return (V)__vpermxor((vector unsigned char)x,(vector unsigned char)pmpm,perm);
|
Chris@42
|
177 #else
|
Chris@42
|
178 /* The safe option */
|
Chris@42
|
179 return FLIP_RI(VCONJ(x));
|
Chris@42
|
180 #endif
|
Chris@42
|
181 }
|
Chris@42
|
182
|
Chris@42
|
183 /* FMA support */
|
Chris@42
|
184 #define VFMA(a, b, c) vec_madd(a,b,c)
|
Chris@42
|
185 #define VFNMS(a, b, c) vec_nmsub(a,b,c)
|
Chris@42
|
186 #define VFMS(a, b, c) vec_msub(a,b,c)
|
Chris@42
|
187 #define VFMAI(b, c) VADD(c, VBYI(b))
|
Chris@42
|
188 #define VFNMSI(b, c) VSUB(c, VBYI(b))
|
Chris@42
|
189 #define VFMACONJ(b,c) VADD(VCONJ(b),c)
|
Chris@42
|
190 #define VFMSCONJ(b,c) VSUB(VCONJ(b),c)
|
Chris@42
|
191 #define VFNMSCONJ(b,c) VSUB(c, VCONJ(b))
|
Chris@42
|
192
|
Chris@42
|
193 static inline V VZMUL(V tx, V sr)
|
Chris@42
|
194 {
|
Chris@42
|
195 V tr = VDUPL(tx);
|
Chris@42
|
196 V ti = VDUPH(tx);
|
Chris@42
|
197 tr = VMUL(sr, tr);
|
Chris@42
|
198 sr = VBYI(sr);
|
Chris@42
|
199 return VFMA(ti, sr, tr);
|
Chris@42
|
200 }
|
Chris@42
|
201
|
Chris@42
|
202 static inline V VZMULJ(V tx, V sr)
|
Chris@42
|
203 {
|
Chris@42
|
204 V tr = VDUPL(tx);
|
Chris@42
|
205 V ti = VDUPH(tx);
|
Chris@42
|
206 tr = VMUL(sr, tr);
|
Chris@42
|
207 sr = VBYI(sr);
|
Chris@42
|
208 return VFNMS(ti, sr, tr);
|
Chris@42
|
209 }
|
Chris@42
|
210
|
Chris@42
|
211 static inline V VZMULI(V tx, V sr)
|
Chris@42
|
212 {
|
Chris@42
|
213 V tr = VDUPL(tx);
|
Chris@42
|
214 V ti = VDUPH(tx);
|
Chris@42
|
215 ti = VMUL(ti, sr);
|
Chris@42
|
216 sr = VBYI(sr);
|
Chris@42
|
217 return VFMS(tr, sr, ti);
|
Chris@42
|
218 }
|
Chris@42
|
219
|
Chris@42
|
220 static inline V VZMULIJ(V tx, V sr)
|
Chris@42
|
221 {
|
Chris@42
|
222 V tr = VDUPL(tx);
|
Chris@42
|
223 V ti = VDUPH(tx);
|
Chris@42
|
224 ti = VMUL(ti, sr);
|
Chris@42
|
225 sr = VBYI(sr);
|
Chris@42
|
226 return VFMA(tr, sr, ti);
|
Chris@42
|
227 }
|
Chris@42
|
228
|
Chris@42
|
229 /* twiddle storage #1: compact, slower */
|
Chris@42
|
230 #ifdef FFTW_SINGLE
|
Chris@42
|
231 # define VTW1(v,x) \
|
Chris@42
|
232 {TW_COS, v, x}, {TW_COS, v+1, x}, {TW_SIN, v, x}, {TW_SIN, v+1, x}
|
Chris@42
|
233 static inline V BYTW1(const R *t, V sr)
|
Chris@42
|
234 {
|
Chris@42
|
235 V tx = LDA(t,0,t);
|
Chris@42
|
236 V tr = UNPCKH(tx, tx);
|
Chris@42
|
237 V ti = UNPCKL(tx, tx);
|
Chris@42
|
238 tr = VMUL(tr, sr);
|
Chris@42
|
239 sr = VBYI(sr);
|
Chris@42
|
240 return VFMA(ti, sr, tr);
|
Chris@42
|
241 }
|
Chris@42
|
242 static inline V BYTWJ1(const R *t, V sr)
|
Chris@42
|
243 {
|
Chris@42
|
244 V tx = LDA(t,0,t);
|
Chris@42
|
245 V tr = UNPCKH(tx, tx);
|
Chris@42
|
246 V ti = UNPCKL(tx, tx);
|
Chris@42
|
247 tr = VMUL(tr, sr);
|
Chris@42
|
248 sr = VBYI(sr);
|
Chris@42
|
249 return VFNMS(ti, sr, tr);
|
Chris@42
|
250 }
|
Chris@42
|
251 #else /* !FFTW_SINGLE */
|
Chris@42
|
252 # define VTW1(v,x) {TW_CEXP, v, x}
|
Chris@42
|
253 static inline V BYTW1(const R *t, V sr)
|
Chris@42
|
254 {
|
Chris@42
|
255 V tx = LD(t, 1, t);
|
Chris@42
|
256 return VZMUL(tx, sr);
|
Chris@42
|
257 }
|
Chris@42
|
258 static inline V BYTWJ1(const R *t, V sr)
|
Chris@42
|
259 {
|
Chris@42
|
260 V tx = LD(t, 1, t);
|
Chris@42
|
261 return VZMULJ(tx, sr);
|
Chris@42
|
262 }
|
Chris@42
|
263 #endif
|
Chris@42
|
264 #define TWVL1 (VL)
|
Chris@42
|
265
|
Chris@42
|
266 /* twiddle storage #2: twice the space, faster (when in cache) */
|
Chris@42
|
267 #ifdef FFTW_SINGLE
|
Chris@42
|
268 # define VTW2(v,x) \
|
Chris@42
|
269 {TW_COS, v, x}, {TW_COS, v, x}, {TW_COS, v+1, x}, {TW_COS, v+1, x}, \
|
Chris@42
|
270 {TW_SIN, v, -x}, {TW_SIN, v, x}, {TW_SIN, v+1, -x}, {TW_SIN, v+1, x}
|
Chris@42
|
271 #else /* !FFTW_SINGLE */
|
Chris@42
|
272 # define VTW2(v,x) \
|
Chris@42
|
273 {TW_COS, v, x}, {TW_COS, v, x}, {TW_SIN, v, -x}, {TW_SIN, v, x}
|
Chris@42
|
274 #endif
|
Chris@42
|
275 #define TWVL2 (2 * VL)
|
Chris@42
|
276 static inline V BYTW2(const R *t, V sr)
|
Chris@42
|
277 {
|
Chris@42
|
278 V si = FLIP_RI(sr);
|
Chris@42
|
279 V ti = LDA(t+2*VL,0,t);
|
Chris@42
|
280 V tt = VMUL(ti, si);
|
Chris@42
|
281 V tr = LDA(t,0,t);
|
Chris@42
|
282 return VFMA(tr, sr, tt);
|
Chris@42
|
283 }
|
Chris@42
|
284 static inline V BYTWJ2(const R *t, V sr)
|
Chris@42
|
285 {
|
Chris@42
|
286 V si = FLIP_RI(sr);
|
Chris@42
|
287 V tr = LDA(t,0,t);
|
Chris@42
|
288 V tt = VMUL(tr, sr);
|
Chris@42
|
289 V ti = LDA(t+2*VL,0,t);
|
Chris@42
|
290 return VFNMS(ti, si, tt);
|
Chris@42
|
291 }
|
Chris@42
|
292
|
Chris@42
|
293 /* twiddle storage #3 */
|
Chris@42
|
294 #ifdef FFTW_SINGLE
|
Chris@42
|
295 # define VTW3(v,x) {TW_CEXP, v, x}, {TW_CEXP, v+1, x}
|
Chris@42
|
296 # define TWVL3 (VL)
|
Chris@42
|
297 #else
|
Chris@42
|
298 # define VTW3(v,x) VTW1(v,x)
|
Chris@42
|
299 # define TWVL3 TWVL1
|
Chris@42
|
300 #endif
|
Chris@42
|
301
|
Chris@42
|
302 /* twiddle storage for split arrays */
|
Chris@42
|
303 #ifdef FFTW_SINGLE
|
Chris@42
|
304 # define VTWS(v,x) \
|
Chris@42
|
305 {TW_COS, v, x}, {TW_COS, v+1, x}, {TW_COS, v+2, x}, {TW_COS, v+3, x}, \
|
Chris@42
|
306 {TW_SIN, v, x}, {TW_SIN, v+1, x}, {TW_SIN, v+2, x}, {TW_SIN, v+3, x}
|
Chris@42
|
307 #else
|
Chris@42
|
308 # define VTWS(v,x) \
|
Chris@42
|
309 {TW_COS, v, x}, {TW_COS, v+1, x}, {TW_SIN, v, x}, {TW_SIN, v+1, x}
|
Chris@42
|
310 #endif
|
Chris@42
|
311 #define TWVLS (2 * VL)
|
Chris@42
|
312
|
Chris@42
|
313 #define VLEAVE() /* nothing */
|
Chris@42
|
314
|
Chris@42
|
315 #include "simd-common.h"
|