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
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
|
// 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 <xmmintrin.h>
#ifdef __SSE4_1__
// Intel's SSE 4.1
# include <smmintrin.h>
#endif
#ifdef __SSE4A__
// AMD's SSE 4a
# include <ammintrin.h>
#endif
#ifdef __FMA4__
# include <fma4intrin.h>
#endif
#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
// Integer and boolean types corresponding to this real type
#define CCTK_INTEGER4 CCTK_REAL4
#define CCTK_BOOLEAN4 CCTK_REAL4
#define CCTK_INTEGER4_VEC CCTK_REAL4_VEC
#define CCTK_BOOLEAN4_VEC CCTK_REAL4_VEC
union k4const_t {
unsigned i[4];
float f[4];
__m128i vi;
__m128 vf;
};
#define K4_ZERO 0x00000000UL
#define K4_IMIN 0x80000000UL
#define K4_IMAX 0x7fffffffUL
// 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 x__=(x_); \
CCTK_REAL4_VEC const x=x__; \
_mm_shuffle_ps(x,x, _MM_SHUFFLE(2,3,0,1)); \
})
#define vec4_swap2301(x_) \
({ \
CCTK_REAL4_VEC const x__=(x_); \
CCTK_REAL4_VEC const x=x__; \
_mm_shuffle_ps(x,x, _MM_SHUFFLE(1,0,3,2)); \
})
#define vec4_swap3210(x_) \
({ \
CCTK_REAL4_VEC const x__=(x_); \
CCTK_REAL4_VEC const x=x__; \
_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 x__=(x_); \
CCTK_REAL4_VEC const x=x__; \
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 x__=(x_); \
CCTK_REAL4_VEC const x=x__; \
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& p__=(p_); \
CCTK_REAL4 const& p=p__; \
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_shuffle_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& p__=(p_); \
CCTK_REAL4 const& p=p__; \
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& p__=(p_); \
CCTK_REAL4 const& p=p__; \
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_shuffle_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& p__=(p_); \
CCTK_REAL4 const& p=p__; \
(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 vector to memory (aligned and non-temporal); this stores to
// a reference to a scalar
#define vec4_store(p,x) (_mm_store_ps(&(p),x))
#define vec4_storeu(p,x) (_mm_storeu_ps(&(p),x))
#if ! VECTORISE_STREAMING_STORES
# define vec4_store_nta(p,x) vec4_store(p,x)
#else
# define vec4_store_nta(p,x) (_mm_stream_ps(&(p),x))
#endif
// Store a partial vector (aligned and non-temporal)
#define vec4_store_partial_prepare(i,imin,imax) \
int v4stp_lo_skip = (imin)-(i); \
int v4stp_hi_skip = (i)+CCTK_REAL_VEC_SIZE-(imax); \
if (CCTK_BUILTIN_EXPECT(v4stp_lo_skip < 0, true)) v4stp_lo_skip = 0; \
if (CCTK_BUILTIN_EXPECT(v4stp_hi_skip < 0, true)) v4stp_hi_skip = 0;
// Ignoring VECTORISE_STREAMING_STORES for partial stores
#define vec4_store_nta_partial(p_,x_) \
({ \
CCTK_REAL4& p__=(p_); \
CCTK_REAL4& p=p__; \
CCTK_REAL4_VEC const x__=(x_); \
CCTK_REAL4_VEC const x=x__; \
if (CCTK_BUILTIN_EXPECT(v4stp_lo_skip==0 and v4stp_hi_skip==0, true)) { \
vec4_store_nta(p,x); \
} else { \
/* these cases fall through */ \
switch (v4stp_lo_skip) { \
case 0: \
(&p)[0] = vec4_elt0(x); \
case 1: \
if (v4stp_hi_skip>=3) break; \
(&p)[1] = vec4_elt1(x); \
case 2: \
if (v4stp_hi_skip>=2) break; \
(&p)[2] = vec4_elt2(x); \
case 3: \
if (v4stp_hi_skip>=1) break; \
(&p)[3] = vec4_elt3(x); \
} \
} \
})
// Ignoring VECTORISE_STREAMING_STORES for partial stores
#define vec4_store_nta_partial_lo(p_,x_,n) \
({ \
CCTK_REAL4 & p__=(p_); \
CCTK_REAL4_VEC const x__=(x_); \
CCTK_REAL4 & p=p__; \
CCTK_REAL4_VEC const x=x__; \
/* these cases fall through */ \
switch (n) { \
case 3: (&p)[2] = vec4_elt2(x); \
case 2: (&p)[1] = vec4_elt1(x); \
case 1: (&p)[0] = vec4_elt0(x); \
} \
})
#define vec4_store_nta_partial_hi(p_,x_,n) \
({ \
CCTK_REAL4 & p__=(p_); \
CCTK_REAL4_VEC const x__=(x_); \
CCTK_REAL4 & p=p__; \
CCTK_REAL4_VEC const x=x__; \
/* these cases fall through */ \
switch (n) { \
case 3: (&p)[1]=vec4_elt1(x); \
case 2: (&p)[2]=vec4_elt2(x); \
case 1: (&p)[3]=vec4_elt3(x); \
} \
})
#define vec4_store_nta_partial_mid(p_,x_,nlo,nhi) \
({ \
CCTK_REAL4 & p__=(p_); \
CCTK_REAL4_VEC const x__=(x_); \
CCTK_REAL4 & p=p__; \
CCTK_REAL4_VEC const x=x__; \
/* these cases fall through */ \
switch (nhi) { \
case 3: if (nlo<2) break; (&p)[1] = vec4_elt1(x); \
case 2: if (nlo<3) break; (&p)[2] = vec4_elt2(x); \
} \
})
// Functions and operators
static const k4const_t k4sign_mask = {{ K4_IMIN, K4_IMIN, K4_IMIN, K4_IMIN, }};
// Operators
#define k4neg(x) (_mm_xor_ps(k4sign_mask.vf,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
#ifdef __FMA4__
# define k4madd(x,y,z) (_mm_macc_ps(x,y,z))
# define k4msub(x,y,z) (_mm_msub_ps(x,y,z))
# define k4nmadd(x,y,z) (_mm_nmsub_ps(x,y,z))
# define k4nmsub(x,y,z) (_mm_nmacc_ps(x,y,z))
#else
# 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)))
#endif
// Cheap functions
#define k4copysign(x,y) \
(_mm_or_ps(_mm_andnot_ps(k4sign_mask.vf,x), \
_mm_and_ps(k4sign_mask.vf,y)))
#define k4fabs(x) (_mm_andnot_ps(k4sign_mask.vf,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.vf,x))
static const k4const_t k4zero = { f: { 0.0f, 0.0f, 0.0f, 0.0f, }};
static const k4const_t k4one = { f: { 1.0f, 1.0f, 1.0f, 1.0f, }};
#define k4sgn(x_) \
({ \
CCTK_REAL_VEC const x__=(x_); \
CCTK_REAL_VEC const x=x__; \
CCTK_REAL_VEC const iszero = _mm_cmpeq_ps(k4zero.vf, x); \
CCTK_REAL_VEC const sign = _mm_and_ps(k4sign_mask.vf, x); \
CCTK_REAL_VEC const signedone = _mm_or_ps(k4one.vf, sign); \
k4ifthen(iszero, k4zero.vf, signedone); \
})
// TODO: maybe use rsqrt and Newton-Raphson
#define k4sqrt(x) (_mm_sqrt_ps(x))
// Expensive functions
#define K4REPL(f,x_) \
({ \
CCTK_REAL4_VEC const x__=(x_); \
CCTK_REAL4_VEC const x=x__; \
vec4_set(f(vec4_elt0(x)), \
f(vec4_elt1(x)), \
f(vec4_elt2(x)), \
f(vec4_elt3(x))); \
})
#define K4REPL2S(f,x_,a_) \
({ \
CCTK_REAL4_VEC const x__=(x_); \
CCTK_REAL4 const a__=(a_); \
CCTK_REAL4_VEC const x=x__; \
CCTK_REAL4 const a=a__; \
vec4_set(f(vec4_elt0(x),a), \
f(vec4_elt1(x),a), \
f(vec4_elt2(x),a), \
f(vec4_elt3(x),a)); \
})
#define K4REPL2(f,x_,y_) \
({ \
CCTK_REAL4_VEC const x__=(x_); \
CCTK_REAL4_VEC const y__=(y_); \
CCTK_REAL4_VEC const x=x__; \
CCTK_REAL4_VEC const y=y__; \
vec4_set(f(vec4_elt0(x),vec4_elt0(y)), \
f(vec4_elt1(x),vec4_elt1(y)), \
f(vec4_elt2(x),vec4_elt2(y)), \
f(vec4_elt3(x),vec4_elt3(y))); \
})
#define k4acos(x) K4REPL(acosf,x)
#define k4acosh(x) K4REPL(acoshf,x)
#define k4asin(x) K4REPL(asinf,x)
#define k4asinh(x) K4REPL(asinhf,x)
#define k4atan(x) K4REPL(atanf,x)
#define k4atan2(x,y) K4REPL2(atan2f,x,y)
#define k4atanh(x) K4REPL(atanhf,x)
#define k4cos(x) K4REPL(cosf,x)
#define k4cosh(x) K4REPL(coshf,x)
#define k4exp(x) K4REPL(expf,x)
#define k4log(x) K4REPL(logf,x)
#define k4pow(x,a) K4REPL2S(powf,x,a)
#define k4sin(x) K4REPL(sinf,x)
#define k4sinh(x) K4REPL(sinhf,x)
#define k4tan(x) K4REPL(tanf,x)
#define k4tanh(x) K4REPL(tanhf,x)
static const k4const_t k4lfalse_ = {{ +0U, +0U, +0U, +0U, }};
static const k4const_t k4ltrue_ = {{ -1U, -1U, -1U, -1U, }};
#define k4lfalse (k4lfalse_.vf)
#define k4ltrue (k4ltrue_.vf)
#define k4lnot(x) (_mm_xor_ps(k4ltrue,x))
#define k4land(x,y) (_mm_and_ps(x,y))
#define k4lor(x,y) (_mm_or_ps(x,y))
#define k4lxor(x,y) (_mm_xor_ps(x,y))
#ifdef __SSE4_1__
# define k4ifthen(x,y,z) (_mm_blendv_ps(z,y,x))
#elif 0
# ifdef __cplusplus
# define k4signbit(x) ({ using namespace std; signbit(x); })
# else
# define k4signbit(x) (signbitf(x))
# endif
# define k4ifthen(x,y,z) \
({ \
CCTK_REAL4_VEC const x__=(x_); \
CCTK_REAL4_VEC const y__=(y_); \
CCTK_REAL4_VEC const z__=(z_); \
CCTK_REAL4_VEC const x=x__; \
CCTK_REAL4_VEC const y=y__; \
CCTK_REAL4_VEC const z=z__; \
vec4_set(k4signbit(vec4_elt0(x)) ? vec4_elt0(y) : vec4_elt0(z), \
k4signbit(vec4_elt1(x)) ? vec4_elt1(y) : vec4_elt1(z), \
k4signbit(vec4_elt2(x)) ? vec4_elt2(y) : vec4_elt2(z), \
k4signbit(vec4_elt3(x)) ? vec4_elt3(y) : vec4_elt3(z)); \
})
#elif 0
// We don't need to shift -- the condition (mask) will be either all
// zeros or all ones
# define k4ifthen(x_,y_,z_) \
({ \
CCTK_REAL4_VEC const x__=(x_); \
CCTK_REAL4_VEC const y__=(y_); \
CCTK_REAL4_VEC const z__=(z_); \
CCTK_REAL4_VEC const x=x__; \
CCTK_REAL4_VEC const y=y__; \
CCTK_REAL4_VEC const z=z__; \
CCTK_REAL4_VEC const mask = \
(__m128)_mm_srai_epi32((__m128i)x, 31); \
/* (z & ~mask) | (y & mask) */ \
_mm_or_ps(_mm_andnot_ps(mask, z), _mm_and_ps(mask, y)); \
})
#else
# define k4ifthen(x_,y_,z_) \
({ \
CCTK_REAL4_VEC const x__=(x_); \
CCTK_REAL4_VEC const y__=(y_); \
CCTK_REAL4_VEC const z__=(z_); \
CCTK_REAL4_VEC const x=x__; \
CCTK_REAL4_VEC const y=y__; \
CCTK_REAL4_VEC const z=z__; \
/* (z & ~mask) | (y & mask) where imask = ~mask */ \
_mm_or_ps(_mm_and_ps(x, y), _mm_andnot_ps(x, z)); \
})
#endif
#define k4cmpeq(x,y) (_mm_cmpeq_ps(x,y))
#define k4cmpne(x,y) (_mm_cmpneq_ps(x,y))
#define k4cmpgt(x,y) (_mm_cmpgt_ps(x,y))
#define k4cmpge(x,y) (_mm_cmpge_ps(x,y))
#define k4cmplt(x,y) (_mm_cmplt_ps(x,y))
#define k4cmple(x,y) (_mm_cmple_ps(x,y))
|
