1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
|
// Vectorise using Intel's or AMD's SSE
// Use the type __m128 directly, without introducing a wrapper class
// Use macros instead of inline functions
#include <assert.h>
#include <math.h>
#include <x86intrin.h>
#ifdef __SSE4_1__
# define vec4_architecture_SSE4_1 "+SSE4.1"
#else
# define vec4_architecture_SSE4_1 ""
#endif
#ifdef __SSE4A__
# define vec4_architecture_SSE4a "+SSE4A"
#else
# define vec4_architecture_SSE4a ""
#endif
#ifdef __FMA4__
# define vec4_architecture_FMA4 "+FMA4"
#else
# define vec4_architecture_FMA4 ""
#endif
#define vec4_architecture "SSE" vec4_architecture_SSE4_1 vec4_architecture_SSE4a vec4_architecture_FMA4 " (32-bit precision)"
// Vector type corresponding to CCTK_REAL
#define CCTK_REAL4_VEC __m128
// Number of vector elements in a CCTK_REAL_VEC
#define CCTK_REAL4_VEC_SIZE 4
// Create vectors, extract vector elements
#define vec4_set1(a) (_mm_set1_ps(a))
#define vec4_set(a,b,c,d) (_mm_set_ps(d,c,b,a)) // note reversed arguments
// original order is 0123
#define vec4_swap1032(x_) \
({ \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
_mm_shuffle_ps(x,x, _MM_SHUFFLE(2,3,0,1)); \
})
#define vec4_swap2301(x_) \
({ \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
_mm_shuffle_ps(x,x, _MM_SHUFFLE(1,0,3,2)); \
})
#define vec4_swap3210(x_) \
({ \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
_mm_shuffle_ps(x,x, _MM_SHUFFLE(0,1,2,3)); \
})
#if defined(__PGI)
// _mm_cvtss_f32 does not exist on PGI compilers
# define vec4_elt0(x) \
({ \
CCTK_REAL4 a; \
asm ("" : "=x" (a) : "0" (x)); \
a; \
})
#else
# define vec4_elt0(x) (_mm_cvtss_f32(x)) // this is a no-op
#endif
#define vec4_elt1(x) vec4_elt0(vec4_swap1032(x))
#define vec4_elt2(x) vec4_elt0(vec4_swap2301(x))
#define vec4_elt3(x) vec4_elt0(vec4_swap3210(x))
#if defined(__PGI)
# define vec4_elt(x_,d) \
({ \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
CCTK_REAL4 a; \
if (d==0) a=vec4_elt0(x); \
else if (d==1) a=vec4_elt1(x); \
else if (d==2) a=vec4_elt2(x); \
else if (d==3) a=vec4_elt3(x); \
a; \
})
#else
# define vec4_elt(x_,d) \
({ \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
CCTK_REAL4 a; \
switch (d) { \
case 0: a=vec4_elt0(x); break; \
case 1: a=vec4_elt1(x); break; \
case 2: a=vec4_elt2(x); break; \
case 3: a=vec4_elt3(x); break; \
} \
a; \
})
#endif
// Load and store vectors
// Load a vector from memory (aligned and unaligned); this loads from
// a reference to a scalar
#define vec4_load(p) (_mm_load_ps(&(p)))
#define vec4_loadu(p) (_mm_loadu_ps(&(p)))
#if ! VECTORISE_ALWAYS_USE_ALIGNED_LOADS
# define vec4_load_off1(p) vec_loadu(p)
# define vec4_load_off2(p) vec_loadu(p)
# define vec4_load_off3(p) vec_loadu(p)
#else
# define vec4_load_off1(p_) \
({ \
CCTK_REAL4 const& pp=(p_); \
CCTK_REAL4 const& p=pp; \
CCTK_REAL4_VEC const lo=vec4_load((&p)[-1]); \
CCTK_REAL4_VEC const hi=vec4_load((&p)[+3]); \
assert(0); \
CCTK_REAL4_VEC const hi2=_mm_suffle_ps(lo,hi, _MM_SHUFFLE(0,1,2,3)); \
_mm_shuffle_ps(lo,hi2, _MM_SHUFFLE(2,1,3,0)); \
})
# define vec4_load_off2(p_) \
({ \
CCTK_REAL4 const& pp=(p_); \
CCTK_REAL4 const& p=pp; \
CCTK_REAL4_VEC const lo=vec4_load((&p)[-2]); \
CCTK_REAL4_VEC const hi=vec4_load((&p)[+2]); \
_mm_shuffle_ps(lo,hi, _MM_SHUFFLE(1,0,3,2)); \
})
# define vec4_load_off1(p_) \
({ \
CCTK_REAL4 const& pp=(p_); \
CCTK_REAL4 const& p=pp; \
CCTK_REAL4_VEC const lo=vec4_load((&p)[-1]); \
CCTK_REAL4_VEC const hi=vec4_load((&p)[+3]); \
assert(0); \
CCTK_REAL4_VEC const lo2=_mm_suffle_ps(lo,hi, _MM_SHUFFLE(0,1,2,3)); \
_mm_shuffle_ps(lo2,hi, _MM_SHUFFLE(3,0,2,1)); \
})
#endif
// Load a vector from memory that may or may not be aligned, as
// decided by the offset off and the vector size
#if VECTORISE_ALWAYS_USE_UNALIGNED_LOADS
// Implementation: Always use unaligned load
# define vec4_loadu_maybe(off,p) vec4_loadu(p)
# define vec4_loadu_maybe3(off1,off2,off3,p) vec4_loadu(p)
#else
# define vec4_loadu_maybe(off,p_) \
({ \
CCTK_REAL4 const& pp=(p_); \
CCTK_REAL4 const& p=pp; \
(off) % CCTK_REAL4_VEC_SIZE == 0 ? \
vec4_load(p) : \
vec4_loadu(p); \
})
# if VECTORISE_ALIGNED_ARRAYS
// Assume all array x sizes are multiples of the vector size
# define vec4_loadu_maybe3(off1,off2,off3,p) \
vec4_loadu_maybe(off1,p)
# else
# define vec4_loadu_maybe3(off1,off2,off3,p) \
vec4_loadu_maybe((off1)|(off2)|(off3),p)
# endif
#endif
// Store a lower or higher partial vector (aligned and non-temporal);
// the non-temporal hint is probably ignored
#if ! VECTORISE_STREAMING_STORES || ! defined(__SSE4A__)
# define vec4_store_nta_partial_lo(p_,x_,n) \
({ \
CCTK_REAL4 const& pp=(p_); \
CCTK_REAL4 const& p=pp; \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
switch (n) { \
case 1: (&p)[0]=vec4_elt0(x); break; \
case 2: _mm_storel_ps(&p,x); break; \
case 3: _mm_storel_ps(&p,x); (&p)[2]=vec4_elt2(x); break; \
} \
})
# define vec4_store_nta_partial_hi(p_,x_,n) \
({ \
CCTK_REAL4 const& pp=(p_); \
CCTK_REAL4 const& p=pp; \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
switch (n) { \
case 1: (&p)[3]=vec4_elt3(x); break; \
case 2: _mm_storeh_ps(&p+2,x); break; \
case 3: _mm_storeh_ps(&p+2,x); (&p)[1]=vec4_elt1(x); break; \
} \
})
#else
# define vec4_store_nta_partial_lo(p_,x_,n) \
({ \
CCTK_REAL4 const& pp=(p_); \
CCTK_REAL4 const& p=pp; \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
switch (n) { \
case 1: \
_mm_stream_ss(&p,x); \
break; \
case 2: \
_mm_storel_ps(&p,x); \
break; \
case 3: \
_mm_storel_ps(&p,x); \
_mm_stream_ss(&p+2, vec4_swap2301(x)); \
break; \
} \
})
# define vec4_store_nta_partial_hi(p_,x_,n) \
({ \
CCTK_REAL4 const& pp=(p_); \
CCTK_REAL4 const& p=pp; \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
switch (n) { \
case 1: \
_mm_stream_ss(&p+3, vec4_swap3210(x)); \
break; \
case 2: \
_mm_storeh_ps(&p+2,x); \
break; \
case 3: \
_mm_storeh_ps(&p+2,x); \
_mm_stream_ss(&p+1, vec4_swap1032(x)); \
break; \
} \
})
#endif
// Functions and operators
static const union {
unsigned i[4];
__m128 v;
} k4sign_mask_union = {{ 0x80000000U, 0x80000000U, 0x80000000U, 0x80000000U }};
#define k4sign_mask (k4sign_mask_union.v)
// Operators
#define k4pos(x) (x)
#define k4neg(x) (_mm_xor_ps(k4sign_mask,x))
// #define k4inv(x)
// TODO: provide k4inv via rcp and Newton-Raphson
// This is described in AMD's publication 47414.
// This should apply for AVX as well.
#define k4add(x,y) (_mm_add_ps(x,y))
#define k4sub(x,y) (_mm_sub_ps(x,y))
#define k4mul(x,y) (_mm_mul_ps(x,y))
// TODO: use k4inv and k4mul instead
#define k4div(x,y) (_mm_div_ps(x,y))
// Fused multiply-add, defined as [+-]x*y[+-]z
#define k4madd(x,y,z) (k4add(k4mul(x,y),z))
#define k4msub(x,y,z) (k4sub(k4mul(x,y),z))
#define k4nmadd(x,y,z) (k4sub(k4neg(z),k4mul(x,y)))
#define k4nmsub(x,y,z) (k4sub(z,k4mul(x,y)))
// Cheap functions
#define k4fabs(x) (_mm_andnot_ps(k4sign_mask,x))
#define k4fmax(x,y) (_mm_max_ps(x,y))
#define k4fmin(x,y) (_mm_min_ps(x,y))
#define k4fnabs(x) (_mm_or_ps(k4sign_mask,x))
// TODO: maybe use rsqrt and Newton-Raphson
#define k4sqrt(x) (_mm_sqrt_ps(x))
// Expensive functions
#define K4REPL(f,x_) \
({ \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
vec4_set(f(vec4_elt0(x)), \
f(vec4_elt1(x)), \
f(vec4_elt2(x)), \
f(vec4_elt3(x))); \
})
#define K4REPL2(f,x_,a_) \
({ \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
CCTK_REAL4 const aa=(a_); \
CCTK_REAL4 const a=aa; \
vec4_set(f(vec4_elt0(x),a), \
f(vec4_elt1(x),a), \
f(vec4_elt2(x),a), \
f(vec4_elt3(x),a)); \
})
#define k4exp(x) K4REPL(exp,x)
#define k4log(x) K4REPL(log,x)
#define k4pow(x,a) K4REPL2(pow,x,a)
// Choice [sign(x)>0 ? y : z]
#ifdef __SSE4_1__
# define k4ifpos(x,y,z) (_mm_blendv_ps(y,z,x))
#elif 0
# ifdef __cplusplus
# define k4sgn(x) ({ using namespace std; signbit(x); })
# else
# define k4sgn(x) (signbit(x))
# endif
# define k4ifpos(x,y,z) \
({ \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
CCTK_REAL4_VEC const yy=(y_); \
CCTK_REAL4_VEC const y=yy; \
CCTK_REAL4_VEC const zz=(z_); \
CCTK_REAL4_VEC const z=zz; \
vec4_set(k4sgn(vec4_elt0(x)) ? vec4_elt0(z) : vec4_elt0(y), \
k4sgn(vec4_elt1(x)) ? vec4_elt1(z) : vec4_elt1(y), \
k4sgn(vec4_elt2(x)) ? vec4_elt2(z) : vec4_elt2(y), \
k4sgn(vec4_elt3(x)) ? vec4_elt3(z) : vec4_elt3(y)); \
})
#else
# define k4ifpos(x_,y_,z_) \
({ \
CCTK_REAL4_VEC const xx=(x_); \
CCTK_REAL4_VEC const x=xx; \
CCTK_REAL4_VEC const yy=(y_); \
CCTK_REAL4_VEC const y=yy; \
CCTK_REAL4_VEC const zz=(z_); \
CCTK_REAL4_VEC const z=zz; \
CCTK_REAL4_VEC const mask = _mm_srai_epi32(x, 31); \
/* (y & ~mask) | (z & mask) */ \
_mm_or_ps(_mm_andnot_ps(mask, y), _mm_and_ps(mask, z)); \
})
#endif
|
