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
|
// Instantiante some functions to be able to inspect the generated
// machine code
#define VML_NODEBUG
#define restrict __restrict__
#include "vecmathlib.h"
namespace vecmathlib {
template<typename realvec_t, int n>
typename realvec_t::real_t get_elt(realvec_t x)
{
return x[n];
}
template<typename realvec_t, int n>
realvec_t set_elt(realvec_t x, typename realvec_t::real_t a)
{
return x.set_elt(n, a);
}
// template realbuiltinvec<float,1> fabs(realbuiltinvec<float,1> x);
// template realbuiltinvec<float,1> fmin(realbuiltinvec<float,1> x, realbuiltinvec<float,1> y);
// template intbuiltinvec<float,1> lsr(intbuiltinvec<float,1> x, intbuiltinvec<float,1>::int_t n);
// template intbuiltinvec<double,1> lsr(intbuiltinvec<double,1> x, intbuiltinvec<double,1>::int_t n);
// template intbuiltinvec<double,2> lsr(intbuiltinvec<double,2> x, intbuiltinvec<double,2>::int_t n);
// template intbuiltinvec<double,2> lsr(intbuiltinvec<double,2> x, intbuiltinvec<double,2> n);
// template realbuiltinvec<float,1> ifthen(realbuiltinvec<float,1>::boolvec_t c, realbuiltinvec<float,1> x, realbuiltinvec<float,1> y);
// template realbuiltinvec<double,1> ifthen(realbuiltinvec<double,1>::boolvec_t c, realbuiltinvec<double,1> x, realbuiltinvec<double,1> y);
// template realbuiltinvec<float,4> ifthen(realbuiltinvec<float,4>::boolvec_t c, realbuiltinvec<float,4> x, realbuiltinvec<float,4> y);
// template realbuiltinvec<double,2> ifthen(realbuiltinvec<double,2>::boolvec_t c, realbuiltinvec<double,2> x, realbuiltinvec<double,2> y);
#ifdef VECMATHLIB_HAVE_VEC_FLOAT_1
template realvec<float,1> round(realvec<float,1> x);
#endif
#ifdef VECMATHLIB_HAVE_VEC_FLOAT_8
template intvec<float,8> popcount(intvec<float,8>);
#endif
#ifdef VECMATHLIB_HAVE_VEC_DOUBLE_1
template realvec<double,1> exp(realvec<double,1> x);
template realvec<double,1> log(realvec<double,1> x);
template realvec<double,1> sin(realvec<double,1> x);
template realvec<double,1> sqrt(realvec<double,1> x);
template realvec<double,1>::real_t get_elt<realvec<double,1>,0>(realvec<double,1> x);
template realvec<double,1> set_elt<realvec<double,1>,0>(realvec<double,1> x, realvec<double,1>::real_t a);
#endif
#ifdef VECMATHLIB_HAVE_VEC_DOUBLE_2
template realvec<double,2> exp(realvec<double,2> x);
template realvec<double,2> log(realvec<double,2> x);
template realvec<double,2> sin(realvec<double,2> x);
template realvec<double,2> sqrt(realvec<double,2> x);
template realvec<double,2>::real_t get_elt<realvec<double,2>,0>(realvec<double,2>);
template realvec<double,2>::real_t get_elt<realvec<double,2>,1>(realvec<double,2>);
template realvec<double,2> set_elt<realvec<double,2>,0>(realvec<double,2> x, realvec<double,2>::real_t a);
template realvec<double,2> set_elt<realvec<double,2>,1>(realvec<double,2> x, realvec<double,2>::real_t a);
#endif
#ifdef VECMATHLIB_HAVE_VEC_DOUBLE_4
template realvec<double,4> exp(realvec<double,4> x);
template realvec<double,4> log(realvec<double,4> x);
template realvec<double,4> sin(realvec<double,4> x);
template realvec<double,4> sqrt(realvec<double,4> x);
template realvec<double,4>::real_t get_elt<realvec<double,4>,0>(realvec<double,4>);
template realvec<double,4>::real_t get_elt<realvec<double,4>,1>(realvec<double,4>);
template realvec<double,4>::real_t get_elt<realvec<double,4>,2>(realvec<double,4>);
template realvec<double,4>::real_t get_elt<realvec<double,4>,3>(realvec<double,4>);
template realvec<double,4> set_elt<realvec<double,4>,0>(realvec<double,4> x, realvec<double,4>::real_t a);
template realvec<double,4> set_elt<realvec<double,4>,1>(realvec<double,4> x, realvec<double,4>::real_t a);
template realvec<double,4> set_elt<realvec<double,4>,2>(realvec<double,4> x, realvec<double,4>::real_t a);
template realvec<double,4> set_elt<realvec<double,4>,3>(realvec<double,4> x, realvec<double,4>::real_t a);
template intvec<double,4> popcount(intvec<double,4>);
#endif
}
// Various tests to detect auto-vectorization features
#include <cassert>
#include <cstdlib>
using namespace std;
using namespace vecmathlib;
#if defined VECMATHLIB_HAVE_VEC_DOUBLE_4
typedef realvec<double,4> realV;
#elif defined VECMATHLIB_HAVE_VEC_DOUBLE_2
typedef realvec<double,2> realV;
#elif defined VECMATHLIB_HAVE_VEC_FLOAT_8
typedef realvec<float,8> realV;
#elif defined VECMATHLIB_HAVE_VEC_FLOAT_4
typedef realvec<float,4> realV;
#elif defined VECMATHLIB_HAVE_VEC_FLOAT_2
typedef realvec<float,2> realV;
#else
# error "There are no vector types"
#endif
typedef realV::scalar_t real;
const int vecsize = realV::size;
#ifdef __clang__
# define __builtin_assume_aligned(exp, align) (exp)
#endif
// Simple, naive loop adding two arrays
extern "C"
void loop_add(real* a,
real* b,
real* c,
ptrdiff_t n)
{
for (ptrdiff_t i=0; i<n; i+=vecsize) {
realV tmpb = realV::loadu(&b[i]);
realV tmpc = realV::loadu(&c[i]);
realV tmpa = tmpb + tmpc;
storeu(tmpa, &a[i]);
}
}
// Declare pointers as restrict
extern "C"
void loop_add_restrict(real *restrict a,
real *restrict b,
real *restrict c,
ptrdiff_t n)
{
for (ptrdiff_t i=0; i<n; i+=vecsize) {
realV tmpb = realV::loadu(&b[i]);
realV tmpc = realV::loadu(&c[i]);
realV tmpa = tmpb + tmpc;
storeu(tmpa, &a[i]);
}
}
// Declare pointers as restrict and aligned
extern "C"
void loop_add_aligned(real *restrict a,
real *restrict b,
real *restrict c,
ptrdiff_t n)
{
for (ptrdiff_t i=0; i<n; i+=vecsize) {
realV tmpb = realV::loada(&b[i]);
realV tmpc = realV::loada(&c[i]);
realV tmpa = tmpb + tmpc;
storea(tmpa, &a[i]);
}
}
// Reduction loop
extern "C"
real loop_dot_reduce(real *restrict a,
real *restrict b,
ptrdiff_t n)
{
realV sumV = 0.0;
for (ptrdiff_t i=0; i<n; i+=vecsize) {
realV tmpa = realV::loada(&a[i]);
realV tmpb = realV::loada(&b[i]);
sumV += tmpa * tmpb;
}
return sum(sumV);
}
// Loop with a simple if condition (fmax)
extern "C"
void loop_if_simple(real *restrict a,
real *restrict b,
real *restrict c,
ptrdiff_t n)
{
for (ptrdiff_t i=0; i<n; i+=vecsize) {
realV tmpb = realV::loada(&b[i]);
realV tmpc = realV::loada(&c[i]);
realV tmpa = ifthen(tmpb > tmpc, tmpb, tmpc);
storea(tmpa, &a[i]);
}
}
// Loop with a complex if condition (select)
extern "C"
void loop_if(real *restrict a,
real *restrict b,
real *restrict c,
ptrdiff_t n)
{
for (ptrdiff_t i=0; i<n; i+=vecsize) {
realV tmpb = realV::loada(&b[i]);
realV tmpc = realV::loada(&c[i]);
realV tmpa = ifthen(tmpb > realV(0.0), tmpb * tmpc, realV(1.0));
storea(tmpa, &a[i]);
}
}
// Skip ghost points
extern "C"
void loop_add_masked(real *restrict a,
real *restrict b,
real *restrict c,
ptrdiff_t n)
{
for (realV::mask_t mask(1, n-1, 0); mask; ++mask) {
ptrdiff_t i = mask.index();
realV tmpb = realV::loada(&b[i]);
realV tmpc = realV::loada(&c[i]);
realV tmpa = tmpb + tmpc;
storea(tmpa, &a[i], mask);
}
}
|
