1 /*
2 * Copyright Nick Thompson, 2020
3 * Use, modification and distribution are subject to the
4 * Boost Software License, Version 1.0. (See accompanying file
5 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6 */
7
8 #include <cmath>
9 #include <random>
10 #include <benchmark/benchmark.h>
11 #include <boost/random/uniform_real_distribution.hpp>
12 #include <boost/math/special_functions/daubechies_scaling.hpp>
13 #include <boost/math/interpolators/cubic_hermite.hpp>
14 #include <boost/math/interpolators/detail/quintic_hermite_detail.hpp>
15 #include <boost/math/interpolators/detail/septic_hermite_detail.hpp>
16
exponential(benchmark::IterationCount j)17 double exponential(benchmark::IterationCount j)
18 {
19 return std::pow(2, j);
20 }
21
22
23 template<typename Real, int p>
DyadicGrid(benchmark::State & state)24 void DyadicGrid(benchmark::State & state)
25 {
26 int j = state.range(0);
27 size_t s = 0;
28 for (auto _ : state)
29 {
30 auto v = boost::math::daubechies_scaling_dyadic_grid<Real, 4, 0>(j);
31 benchmark::DoNotOptimize(v[0]);
32 s = v.size();
33 }
34
35 state.counters["RAM"] = s*sizeof(Real);
36 state.SetComplexityN(state.range(0));
37 }
38
39 BENCHMARK_TEMPLATE(DyadicGrid, double, 4)->DenseRange(3, 22, 1)->Unit(benchmark::kMillisecond)->Complexity(exponential);
40 //BENCHMARK_TEMPLATE(DyadicGrid, double, 8)->DenseRange(3, 22, 1)->Unit(benchmark::kMillisecond)->Complexity(exponential);
41 //BENCHMARK_TEMPLATE(DyadicGrid, double, 11)->DenseRange(3,22,1)->Unit(benchmark::kMillisecond)->Complexity(exponential);
42
43 uint64_t s[2] = { 0x41, 0x29837592 };
44
rotl(const uint64_t x,int k)45 static inline uint64_t rotl(const uint64_t x, int k) {
46 return (x << k) | (x >> (64 - k));
47 }
48
next(void)49 uint64_t next(void) {
50 const uint64_t s0 = s[0];
51 uint64_t s1 = s[1];
52 const uint64_t result = s0 + s1;
53
54 s1 ^= s0;
55 s[0] = rotl(s0, 55) ^ s1 ^ (s1 << 14); // a, b
56 s[1] = rotl(s1, 36); // c
57
58 return result;
59 }
60
uniform()61 double uniform() {
62 return next()*(1.0/18446744073709551616.0);
63 }
64
65 template<typename Real, int p>
ScalingEvaluation(benchmark::State & state)66 void ScalingEvaluation(benchmark::State & state)
67 {
68 auto phi = boost::math::daubechies_scaling<Real, p>();
69 Real xmax = phi.support().second;
70 Real x = 0;
71 Real step = uniform()/2048;
72 for (auto _ : state)
73 {
74 benchmark::DoNotOptimize(phi(x));
75 x += step;
76 if (x > xmax) {
77 x = 0;
78 step = uniform()/2048;
79 }
80 }
81 }
82
83
84 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 2);
85 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 3);
86 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 4);
87 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 5);
88 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 6);
89 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 7);
90 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 8);
91 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 9);
92 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 10);
93 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 11);
94 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 12);
95 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 13);
96 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 14);
97 BENCHMARK_TEMPLATE(ScalingEvaluation, double, 15);
98
99
100 template<typename Real, int p>
ScalingConstructor(benchmark::State & state)101 void ScalingConstructor(benchmark::State & state)
102 {
103 for (auto _ : state)
104 {
105 benchmark::DoNotOptimize(boost::math::daubechies_scaling<Real, p>());
106 }
107 }
108
109 BENCHMARK_TEMPLATE(ScalingConstructor, float, 2)->Unit(benchmark::kMillisecond);
110 BENCHMARK_TEMPLATE(ScalingConstructor, double, 2)->Unit(benchmark::kMillisecond);
111 BENCHMARK_TEMPLATE(ScalingConstructor, long double, 2)->Unit(benchmark::kMillisecond);
112
113 BENCHMARK_TEMPLATE(ScalingConstructor, float, 3)->Unit(benchmark::kMillisecond);
114 BENCHMARK_TEMPLATE(ScalingConstructor, double, 3)->Unit(benchmark::kMillisecond);
115 BENCHMARK_TEMPLATE(ScalingConstructor, long double, 3)->Unit(benchmark::kMillisecond);
116
117 BENCHMARK_TEMPLATE(ScalingConstructor, float, 4)->Unit(benchmark::kMillisecond);
118 BENCHMARK_TEMPLATE(ScalingConstructor, double, 4)->Unit(benchmark::kMillisecond);
119 BENCHMARK_TEMPLATE(ScalingConstructor, long double, 4)->Unit(benchmark::kMillisecond);
120
121 BENCHMARK_TEMPLATE(ScalingConstructor, float, 5)->Unit(benchmark::kMillisecond);
122 BENCHMARK_TEMPLATE(ScalingConstructor, double, 5)->Unit(benchmark::kMillisecond);
123 BENCHMARK_TEMPLATE(ScalingConstructor, long double, 5)->Unit(benchmark::kMillisecond);
124
125 BENCHMARK_TEMPLATE(ScalingConstructor, float, 11)->Unit(benchmark::kMillisecond);
126 BENCHMARK_TEMPLATE(ScalingConstructor, double, 11)->Unit(benchmark::kMillisecond);
127 BENCHMARK_TEMPLATE(ScalingConstructor, long double, 11)->Unit(benchmark::kMillisecond);
128
129 template<typename Real>
CubicHermite(benchmark::State & state)130 void CubicHermite(benchmark::State & state)
131 {
132 using boost::math::interpolators::cubic_hermite;
133 auto n = state.range(0);
134 std::vector<Real> x(n);
135 std::vector<Real> y(n);
136 std::vector<Real> dydx(n);
137 std::random_device rd;
138 boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
139 x[0] = dis(rd);
140 y[0] = dis(rd);
141 dydx[0] = dis(rd);
142 for (size_t i = 1; i < y.size(); ++i)
143 {
144 x[i] = x[i-1] + dis(rd);
145 y[i] = dis(rd);
146 dydx[i] = dis(rd);
147 }
148 Real x0 = x.front();
149 Real xf = x.back();
150
151 auto qh = cubic_hermite(std::move(x), std::move(y), std::move(dydx));
152 Real t = x0;
153 Real step = uniform()*(xf-x0)/2048;
154 for (auto _ : state)
155 {
156 benchmark::DoNotOptimize(qh(t));
157 t += step;
158 if (t >= xf)
159 {
160 t = x0;
161 step = uniform()*(xf-x0)/2048;
162 }
163 }
164 state.SetComplexityN(state.range(0));
165 }
166
167 BENCHMARK_TEMPLATE(CubicHermite, double)->RangeMultiplier(2)->Range(1<<8, 1<<20)->Complexity(benchmark::oLogN);
168
169 template<typename Real>
CardinalCubicHermite(benchmark::State & state)170 void CardinalCubicHermite(benchmark::State & state)
171 {
172 using boost::math::interpolators::detail::cardinal_cubic_hermite_detail;
173 auto n = state.range(0);
174 std::vector<Real> y(n);
175 std::vector<Real> dydx(n);
176 std::random_device rd;
177 boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
178 for (size_t i = 0; i < y.size(); ++i)
179 {
180 y[i] = uniform();
181 dydx[i] = uniform();
182 }
183
184 Real dx = Real(1)/Real(8);
185 Real x0 = 0;
186 Real xf = x0 + (y.size()-1)*dx;
187
188 auto qh = cardinal_cubic_hermite_detail(std::move(y), std::move(dydx), x0, dx);
189 Real x = x0;
190 Real step = uniform()*(xf-x0)/2048;
191 for (auto _ : state)
192 {
193 benchmark::DoNotOptimize(qh.unchecked_evaluation(x));
194 x += step;
195 if (x >= xf)
196 {
197 x = x0;
198 step = uniform()*(xf-x0)/2048;
199 }
200 }
201 state.SetComplexityN(state.range(0));
202 }
203
204 template<typename Real>
CardinalCubicHermiteAOS(benchmark::State & state)205 void CardinalCubicHermiteAOS(benchmark::State & state)
206 {
207 auto n = state.range(0);
208 std::vector<std::array<Real, 2>> dat(n);
209 std::random_device rd;
210 boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
211 for (size_t i = 0; i < dat.size(); ++i)
212 {
213 dat[i][0] = uniform();
214 dat[i][1] = uniform();
215 }
216
217 using boost::math::interpolators::detail::cardinal_cubic_hermite_detail_aos;
218 Real dx = Real(1)/Real(8);
219 Real x0 = 0;
220 Real xf = x0 + (dat.size()-1)*dx;
221 auto qh = cardinal_cubic_hermite_detail_aos(std::move(dat), x0, dx);
222 Real x = 0;
223 Real step = uniform()*(xf-x0)/2048;
224 for (auto _ : state)
225 {
226 benchmark::DoNotOptimize(qh.unchecked_evaluation(x));
227 x += step;
228 if (x >= xf)
229 {
230 x = x0;
231 step = uniform()*(xf-x0)/2048;
232 }
233 }
234 state.SetComplexityN(state.range(0));
235 }
236
237 BENCHMARK_TEMPLATE(CardinalCubicHermiteAOS, double)->RangeMultiplier(2)->Range(1<<8, 1<<21)->Complexity(benchmark::o1);
238 BENCHMARK_TEMPLATE(CardinalCubicHermite, double)->RangeMultiplier(2)->Range(1<<8, 1<<21)->Complexity(benchmark::o1);
239
240 template<class Real>
SineEvaluation(benchmark::State & state)241 void SineEvaluation(benchmark::State& state)
242 {
243 std::default_random_engine gen;
244 std::uniform_real_distribution<Real> x_dis(0, 3.14159);
245
246 Real x = x_dis(gen);
247 for (auto _ : state)
248 {
249 benchmark::DoNotOptimize(std::sin(x));
250 x += std::numeric_limits<Real>::epsilon();
251 }
252 }
253
254 BENCHMARK_TEMPLATE(SineEvaluation, float);
255 BENCHMARK_TEMPLATE(SineEvaluation, double);
256 BENCHMARK_TEMPLATE(SineEvaluation, long double);
257
258 template<class Real>
ExpEvaluation(benchmark::State & state)259 void ExpEvaluation(benchmark::State& state)
260 {
261 std::default_random_engine gen;
262 std::uniform_real_distribution<Real> x_dis(0, 3.14159);
263
264 Real x = x_dis(gen);
265 for (auto _ : state)
266 {
267 benchmark::DoNotOptimize(std::exp(x));
268 x += std::numeric_limits<Real>::epsilon();
269 }
270 }
271
272 BENCHMARK_TEMPLATE(ExpEvaluation, float);
273 BENCHMARK_TEMPLATE(ExpEvaluation, double);
274 BENCHMARK_TEMPLATE(ExpEvaluation, long double);
275
276 template<class Real>
PowEvaluation(benchmark::State & state)277 void PowEvaluation(benchmark::State& state)
278 {
279 std::default_random_engine gen;
280 std::uniform_real_distribution<Real> x_dis(0, 3.14159);
281
282 Real x = x_dis(gen);
283 for (auto _ : state)
284 {
285 benchmark::DoNotOptimize(std::pow(x, x+1));
286 x += std::numeric_limits<Real>::epsilon();
287 }
288 }
289
290 BENCHMARK_TEMPLATE(PowEvaluation, float);
291 BENCHMARK_TEMPLATE(PowEvaluation, double);
292 BENCHMARK_TEMPLATE(PowEvaluation, long double);
293
294
295 template<typename Real>
CardinalQuinticHermite(benchmark::State & state)296 void CardinalQuinticHermite(benchmark::State & state)
297 {
298 using boost::math::interpolators::detail::cardinal_quintic_hermite_detail;
299 auto n = state.range(0);
300 std::vector<Real> y(n);
301 std::vector<Real> dydx(n);
302 std::vector<Real> d2ydx2(n);
303 std::random_device rd;
304 boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
305 for (size_t i = 0; i < y.size(); ++i)
306 {
307 y[i] = uniform();
308 dydx[i] = uniform();
309 d2ydx2[i] = uniform();
310 }
311
312 Real dx = Real(1)/Real(8);
313 Real x0 = 0;
314 Real xf = x0 + (y.size()-1)*dx;
315
316 auto qh = cardinal_quintic_hermite_detail(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx);
317 Real x = 0;
318 Real step = uniform()*(xf-x0)/2048;
319 for (auto _ : state)
320 {
321 benchmark::DoNotOptimize(qh.unchecked_evaluation(x));
322 x += step;
323 if (x >= xf)
324 {
325 x = x0;
326 step = uniform()*(xf-x0)/2048;
327 }
328 }
329 state.SetComplexityN(state.range(0));
330 }
331
332 template<typename Real>
CardinalQuinticHermiteAOS(benchmark::State & state)333 void CardinalQuinticHermiteAOS(benchmark::State & state)
334 {
335 auto n = state.range(0);
336 std::vector<std::array<Real, 3>> dat(n);
337 std::random_device rd;
338 boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
339 for (size_t i = 0; i < dat.size(); ++i)
340 {
341 dat[i][0] = uniform();
342 dat[i][1] = uniform();
343 dat[i][2] = uniform();
344 }
345
346 using boost::math::interpolators::detail::cardinal_quintic_hermite_detail_aos;
347 Real dx = Real(1)/Real(8);
348 Real x0 = 0;
349 Real xf = x0 + (dat.size()-1)*dx;
350 auto qh = cardinal_quintic_hermite_detail_aos(std::move(dat), x0, dx);
351 Real x = x0;
352 Real step = uniform()*(xf-x0)/2048;
353 for (auto _ : state) {
354 benchmark::DoNotOptimize(qh.unchecked_evaluation(x));
355 x += step;
356 if (x >= xf)
357 {
358 x = x0;
359 step = uniform()*(xf-x0)/2048;
360 }
361 }
362 state.SetComplexityN(state.range(0));
363 }
364
365 BENCHMARK_TEMPLATE(CardinalQuinticHermiteAOS, double)->RangeMultiplier(2)->Range(1<<8, 1<<22)->Complexity(benchmark::o1);
366 BENCHMARK_TEMPLATE(CardinalQuinticHermite, double)->RangeMultiplier(2)->Range(1<<8, 1<<22)->Complexity(benchmark::o1);
367
368 template<typename Real>
SepticHermite(benchmark::State & state)369 void SepticHermite(benchmark::State & state)
370 {
371 using boost::math::interpolators::detail::septic_hermite_detail;
372 auto n = state.range(0);
373 std::vector<Real> x(n);
374 std::vector<Real> y(n);
375 std::vector<Real> dydx(n);
376 std::vector<Real> d2ydx2(n);
377 std::vector<Real> d3ydx3(n);
378 std::random_device rd;
379 boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
380 Real x0 = dis(rd);
381 x[0] = x0;
382 for (size_t i = 1; i < n; ++i)
383 {
384 x[i] = x[i-1] + dis(rd);
385 }
386 for (size_t i = 0; i < y.size(); ++i)
387 {
388 y[i] = dis(rd);
389 dydx[i] = dis(rd);
390 d2ydx2[i] = dis(rd);
391 d3ydx3[i] = dis(rd);
392 }
393
394 Real xf = x.back();
395
396 auto sh = septic_hermite_detail(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
397 Real t = x0;
398 for (auto _ : state)
399 {
400 benchmark::DoNotOptimize(sh(t));
401 t += xf/128;
402 if (t >= xf)
403 {
404 t = x0;
405 }
406 }
407 state.SetComplexityN(state.range(0));
408 }
409
410 BENCHMARK_TEMPLATE(SepticHermite, double)->RangeMultiplier(2)->Range(1<<8, 1<<20)->Complexity();
411
412
413 template<typename Real>
CardinalSepticHermite(benchmark::State & state)414 void CardinalSepticHermite(benchmark::State & state)
415 {
416 using boost::math::interpolators::detail::cardinal_septic_hermite_detail;
417 auto n = state.range(0);
418 std::vector<Real> y(n);
419 std::vector<Real> dydx(n);
420 std::vector<Real> d2ydx2(n);
421 std::vector<Real> d3ydx3(n);
422 std::random_device rd;
423 boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
424 for (size_t i = 0; i < y.size(); ++i)
425 {
426 y[i] = dis(rd);
427 dydx[i] = dis(rd);
428 d2ydx2[i] = dis(rd);
429 d3ydx3[i] = dis(rd);
430 }
431
432 Real dx = Real(1)/Real(8);
433 Real x0 = 0;
434 Real xf = x0 + (y.size()-1)*dx;
435
436 auto sh = cardinal_septic_hermite_detail(std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), x0, dx);
437 Real x = 0;
438 for (auto _ : state)
439 {
440 benchmark::DoNotOptimize(sh.unchecked_evaluation(x));
441 x += xf/128;
442 if (x >= xf)
443 {
444 x = x0;
445 }
446 }
447 state.SetComplexityN(state.range(0));
448 }
449
450 BENCHMARK_TEMPLATE(CardinalSepticHermite, double)->RangeMultiplier(2)->Range(1<<8, 1<<20)->Complexity();
451
452 template<typename Real>
CardinalSepticHermiteAOS(benchmark::State & state)453 void CardinalSepticHermiteAOS(benchmark::State & state)
454 {
455 using boost::math::interpolators::detail::cardinal_septic_hermite_detail_aos;
456 auto n = state.range(0);
457 std::vector<std::array<Real, 4>> data(n);
458 std::random_device rd;
459 boost::random::uniform_real_distribution<Real> dis(Real(0), Real(1));
460 for (size_t i = 0; i < data.size(); ++i)
461 {
462 for (size_t j = 0; j < 4; ++j)
463 {
464 data[i][j] = dis(rd);
465 }
466 }
467
468 Real dx = Real(1)/Real(8);
469 Real x0 = 0;
470 Real xf = x0 + (data.size()-1)*dx;
471
472 auto sh = cardinal_septic_hermite_detail_aos(std::move(data), x0, dx);
473 Real x = 0;
474 for (auto _ : state)
475 {
476 benchmark::DoNotOptimize(sh.unchecked_evaluation(x));
477 x += xf/128;
478 if (x >= xf)
479 {
480 x = x0;
481 }
482 }
483 state.SetComplexityN(state.range(0));
484 }
485
486 BENCHMARK_TEMPLATE(CardinalSepticHermiteAOS, double)->RangeMultiplier(2)->Range(1<<8, 1<<20)->Complexity();
487
488
489 BENCHMARK_MAIN();
490