1 // Copyright 2019 Google LLC
2 //
3 // Licensed under the Apache License, Version 2.0 (the "License");
4 // you may not use this file except in compliance with the License.
5 // You may obtain a copy of the License at
6 //
7 // http://www.apache.org/licenses/LICENSE-2.0
8 //
9 // Unless required by applicable law or agreed to in writing, software
10 // distributed under the License is distributed on an "AS IS" BASIS,
11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 // See the License for the specific language governing permissions and
13 // limitations under the License.
14
15 #include "hwy/nanobenchmark.h"
16
17 #include <stddef.h>
18 #include <stdio.h>
19 #include <stdlib.h> // abort
20 #include <string.h> // memcpy
21 #include <time.h> // clock_gettime
22
23 #include <algorithm> // sort
24 #include <array>
25 #include <atomic>
26 #include <limits>
27 #include <numeric> // iota
28 #include <random>
29 #include <string>
30 #include <vector>
31
32 #if defined(_WIN32) || defined(_WIN64)
33 #ifndef NOMINMAX
34 #define NOMINMAX
35 #endif // NOMINMAX
36 #include <windows.h>
37 #endif
38
39 #if defined(__MACH__)
40 #include <mach/mach.h>
41 #include <mach/mach_time.h>
42 #endif
43
44 #if defined(__HAIKU__)
45 #include <OS.h>
46 #endif
47
48 #include "hwy/base.h"
49 #if HWY_ARCH_PPC
50 #include <sys/platform/ppc.h> // NOLINT __ppc_get_timebase_freq
51 #elif HWY_ARCH_X86
52
53 #if HWY_COMPILER_MSVC
54 #include <intrin.h>
55 #else
56 #include <cpuid.h> // NOLINT
57 #endif // HWY_COMPILER_MSVC
58
59 #endif // HWY_ARCH_X86
60
61 namespace hwy {
62 namespace {
63 namespace timer {
64
65 // Ticks := platform-specific timer values (CPU cycles on x86). Must be
66 // unsigned to guarantee wraparound on overflow.
67 using Ticks = uint64_t;
68
69 // Start/Stop return absolute timestamps and must be placed immediately before
70 // and after the region to measure. We provide separate Start/Stop functions
71 // because they use different fences.
72 //
73 // Background: RDTSC is not 'serializing'; earlier instructions may complete
74 // after it, and/or later instructions may complete before it. 'Fences' ensure
75 // regions' elapsed times are independent of such reordering. The only
76 // documented unprivileged serializing instruction is CPUID, which acts as a
77 // full fence (no reordering across it in either direction). Unfortunately
78 // the latency of CPUID varies wildly (perhaps made worse by not initializing
79 // its EAX input). Because it cannot reliably be deducted from the region's
80 // elapsed time, it must not be included in the region to measure (i.e.
81 // between the two RDTSC).
82 //
83 // The newer RDTSCP is sometimes described as serializing, but it actually
84 // only serves as a half-fence with release semantics. Although all
85 // instructions in the region will complete before the final timestamp is
86 // captured, subsequent instructions may leak into the region and increase the
87 // elapsed time. Inserting another fence after the final RDTSCP would prevent
88 // such reordering without affecting the measured region.
89 //
90 // Fortunately, such a fence exists. The LFENCE instruction is only documented
91 // to delay later loads until earlier loads are visible. However, Intel's
92 // reference manual says it acts as a full fence (waiting until all earlier
93 // instructions have completed, and delaying later instructions until it
94 // completes). AMD assigns the same behavior to MFENCE.
95 //
96 // We need a fence before the initial RDTSC to prevent earlier instructions
97 // from leaking into the region, and arguably another after RDTSC to avoid
98 // region instructions from completing before the timestamp is recorded.
99 // When surrounded by fences, the additional RDTSCP half-fence provides no
100 // benefit, so the initial timestamp can be recorded via RDTSC, which has
101 // lower overhead than RDTSCP because it does not read TSC_AUX. In summary,
102 // we define Start = LFENCE/RDTSC/LFENCE; Stop = RDTSCP/LFENCE.
103 //
104 // Using Start+Start leads to higher variance and overhead than Stop+Stop.
105 // However, Stop+Stop includes an LFENCE in the region measurements, which
106 // adds a delay dependent on earlier loads. The combination of Start+Stop
107 // is faster than Start+Start and more consistent than Stop+Stop because
108 // the first LFENCE already delayed subsequent loads before the measured
109 // region. This combination seems not to have been considered in prior work:
110 // http://akaros.cs.berkeley.edu/lxr/akaros/kern/arch/x86/rdtsc_test.c
111 //
112 // Note: performance counters can measure 'exact' instructions-retired or
113 // (unhalted) cycle counts. The RDPMC instruction is not serializing and also
114 // requires fences. Unfortunately, it is not accessible on all OSes and we
115 // prefer to avoid kernel-mode drivers. Performance counters are also affected
116 // by several under/over-count errata, so we use the TSC instead.
117
118 // Returns a 64-bit timestamp in unit of 'ticks'; to convert to seconds,
119 // divide by InvariantTicksPerSecond.
Start()120 inline Ticks Start() {
121 Ticks t;
122 #if HWY_ARCH_PPC
123 asm volatile("mfspr %0, %1" : "=r"(t) : "i"(268));
124 #elif HWY_ARCH_X86 && HWY_COMPILER_MSVC
125 _ReadWriteBarrier();
126 _mm_lfence();
127 _ReadWriteBarrier();
128 t = __rdtsc();
129 _ReadWriteBarrier();
130 _mm_lfence();
131 _ReadWriteBarrier();
132 #elif HWY_ARCH_X86_64
133 asm volatile(
134 "lfence\n\t"
135 "rdtsc\n\t"
136 "shl $32, %%rdx\n\t"
137 "or %%rdx, %0\n\t"
138 "lfence"
139 : "=a"(t)
140 :
141 // "memory" avoids reordering. rdx = TSC >> 32.
142 // "cc" = flags modified by SHL.
143 : "rdx", "memory", "cc");
144 #elif HWY_ARCH_RVV
145 asm volatile("rdcycle %0" : "=r"(t));
146 #elif defined(_WIN32) || defined(_WIN64)
147 LARGE_INTEGER counter;
148 (void)QueryPerformanceCounter(&counter);
149 t = counter.QuadPart;
150 #elif defined(__MACH__)
151 t = mach_absolute_time();
152 #elif defined(__HAIKU__)
153 t = system_time_nsecs(); // since boot
154 #else // POSIX
155 timespec ts;
156 clock_gettime(CLOCK_MONOTONIC, &ts);
157 t = ts.tv_sec * 1000000000LL + ts.tv_nsec;
158 #endif
159 return t;
160 }
161
Stop()162 inline Ticks Stop() {
163 uint64_t t;
164 #if HWY_ARCH_PPC
165 asm volatile("mfspr %0, %1" : "=r"(t) : "i"(268));
166 #elif HWY_ARCH_X86 && HWY_COMPILER_MSVC
167 _ReadWriteBarrier();
168 unsigned aux;
169 t = __rdtscp(&aux);
170 _ReadWriteBarrier();
171 _mm_lfence();
172 _ReadWriteBarrier();
173 #elif HWY_ARCH_X86_64
174 // Use inline asm because __rdtscp generates code to store TSC_AUX (ecx).
175 asm volatile(
176 "rdtscp\n\t"
177 "shl $32, %%rdx\n\t"
178 "or %%rdx, %0\n\t"
179 "lfence"
180 : "=a"(t)
181 :
182 // "memory" avoids reordering. rcx = TSC_AUX. rdx = TSC >> 32.
183 // "cc" = flags modified by SHL.
184 : "rcx", "rdx", "memory", "cc");
185 #else
186 t = Start();
187 #endif
188 return t;
189 }
190
191 } // namespace timer
192
193 namespace robust_statistics {
194
195 // Sorts integral values in ascending order (e.g. for Mode). About 3x faster
196 // than std::sort for input distributions with very few unique values.
197 template <class T>
CountingSort(T * values,size_t num_values)198 void CountingSort(T* values, size_t num_values) {
199 // Unique values and their frequency (similar to flat_map).
200 using Unique = std::pair<T, int>;
201 std::vector<Unique> unique;
202 for (size_t i = 0; i < num_values; ++i) {
203 const T value = values[i];
204 const auto pos =
205 std::find_if(unique.begin(), unique.end(),
206 [value](const Unique u) { return u.first == value; });
207 if (pos == unique.end()) {
208 unique.push_back(std::make_pair(value, 1));
209 } else {
210 ++pos->second;
211 }
212 }
213
214 // Sort in ascending order of value (pair.first).
215 std::sort(unique.begin(), unique.end());
216
217 // Write that many copies of each unique value to the array.
218 T* HWY_RESTRICT p = values;
219 for (const auto& value_count : unique) {
220 std::fill(p, p + value_count.second, value_count.first);
221 p += value_count.second;
222 }
223 NANOBENCHMARK_CHECK(p == values + num_values);
224 }
225
226 // @return i in [idx_begin, idx_begin + half_count) that minimizes
227 // sorted[i + half_count] - sorted[i].
228 template <typename T>
MinRange(const T * const HWY_RESTRICT sorted,const size_t idx_begin,const size_t half_count)229 size_t MinRange(const T* const HWY_RESTRICT sorted, const size_t idx_begin,
230 const size_t half_count) {
231 T min_range = std::numeric_limits<T>::max();
232 size_t min_idx = 0;
233
234 for (size_t idx = idx_begin; idx < idx_begin + half_count; ++idx) {
235 NANOBENCHMARK_CHECK(sorted[idx] <= sorted[idx + half_count]);
236 const T range = sorted[idx + half_count] - sorted[idx];
237 if (range < min_range) {
238 min_range = range;
239 min_idx = idx;
240 }
241 }
242
243 return min_idx;
244 }
245
246 // Returns an estimate of the mode by calling MinRange on successively
247 // halved intervals. "sorted" must be in ascending order. This is the
248 // Half Sample Mode estimator proposed by Bickel in "On a fast, robust
249 // estimator of the mode", with complexity O(N log N). The mode is less
250 // affected by outliers in highly-skewed distributions than the median.
251 // The averaging operation below assumes "T" is an unsigned integer type.
252 template <typename T>
ModeOfSorted(const T * const HWY_RESTRICT sorted,const size_t num_values)253 T ModeOfSorted(const T* const HWY_RESTRICT sorted, const size_t num_values) {
254 size_t idx_begin = 0;
255 size_t half_count = num_values / 2;
256 while (half_count > 1) {
257 idx_begin = MinRange(sorted, idx_begin, half_count);
258 half_count >>= 1;
259 }
260
261 const T x = sorted[idx_begin + 0];
262 if (half_count == 0) {
263 return x;
264 }
265 NANOBENCHMARK_CHECK(half_count == 1);
266 const T average = (x + sorted[idx_begin + 1] + 1) / 2;
267 return average;
268 }
269
270 // Returns the mode. Side effect: sorts "values".
271 template <typename T>
Mode(T * values,const size_t num_values)272 T Mode(T* values, const size_t num_values) {
273 CountingSort(values, num_values);
274 return ModeOfSorted(values, num_values);
275 }
276
277 template <typename T, size_t N>
278 T Mode(T (&values)[N]) {
279 return Mode(&values[0], N);
280 }
281
282 // Returns the median value. Side effect: sorts "values".
283 template <typename T>
Median(T * values,const size_t num_values)284 T Median(T* values, const size_t num_values) {
285 NANOBENCHMARK_CHECK(!values->empty());
286 std::sort(values, values + num_values);
287 const size_t half = num_values / 2;
288 // Odd count: return middle
289 if (num_values % 2) {
290 return values[half];
291 }
292 // Even count: return average of middle two.
293 return (values[half] + values[half - 1] + 1) / 2;
294 }
295
296 // Returns a robust measure of variability.
297 template <typename T>
MedianAbsoluteDeviation(const T * values,const size_t num_values,const T median)298 T MedianAbsoluteDeviation(const T* values, const size_t num_values,
299 const T median) {
300 NANOBENCHMARK_CHECK(num_values != 0);
301 std::vector<T> abs_deviations;
302 abs_deviations.reserve(num_values);
303 for (size_t i = 0; i < num_values; ++i) {
304 const int64_t abs = std::abs(int64_t(values[i]) - int64_t(median));
305 abs_deviations.push_back(static_cast<T>(abs));
306 }
307 return Median(abs_deviations.data(), num_values);
308 }
309
310 } // namespace robust_statistics
311 } // namespace
312 namespace platform {
313 namespace {
314
315 // Prevents the compiler from eliding the computations that led to "output".
316 template <class T>
PreventElision(T && output)317 inline void PreventElision(T&& output) {
318 #if HWY_COMPILER_MSVC == 0
319 // Works by indicating to the compiler that "output" is being read and
320 // modified. The +r constraint avoids unnecessary writes to memory, but only
321 // works for built-in types (typically FuncOutput).
322 asm volatile("" : "+r"(output) : : "memory");
323 #else
324 // MSVC does not support inline assembly anymore (and never supported GCC's
325 // RTL constraints). Self-assignment with #pragma optimize("off") might be
326 // expected to prevent elision, but it does not with MSVC 2015. Type-punning
327 // with volatile pointers generates inefficient code on MSVC 2017.
328 static std::atomic<T> dummy(T{});
329 dummy.store(output, std::memory_order_relaxed);
330 #endif
331 }
332
333 #if HWY_ARCH_X86
334
Cpuid(const uint32_t level,const uint32_t count,uint32_t * HWY_RESTRICT abcd)335 void Cpuid(const uint32_t level, const uint32_t count,
336 uint32_t* HWY_RESTRICT abcd) {
337 #if HWY_COMPILER_MSVC
338 int regs[4];
339 __cpuidex(regs, level, count);
340 for (int i = 0; i < 4; ++i) {
341 abcd[i] = regs[i];
342 }
343 #else
344 uint32_t a;
345 uint32_t b;
346 uint32_t c;
347 uint32_t d;
348 __cpuid_count(level, count, a, b, c, d);
349 abcd[0] = a;
350 abcd[1] = b;
351 abcd[2] = c;
352 abcd[3] = d;
353 #endif
354 }
355
BrandString()356 std::string BrandString() {
357 char brand_string[49];
358 std::array<uint32_t, 4> abcd;
359
360 // Check if brand string is supported (it is on all reasonable Intel/AMD)
361 Cpuid(0x80000000U, 0, abcd.data());
362 if (abcd[0] < 0x80000004U) {
363 return std::string();
364 }
365
366 for (size_t i = 0; i < 3; ++i) {
367 Cpuid(static_cast<uint32_t>(0x80000002U + i), 0, abcd.data());
368 memcpy(brand_string + i * 16, abcd.data(), sizeof(abcd));
369 }
370 brand_string[48] = 0;
371 return brand_string;
372 }
373
374 // Returns the frequency quoted inside the brand string. This does not
375 // account for throttling nor Turbo Boost.
NominalClockRate()376 double NominalClockRate() {
377 const std::string& brand_string = BrandString();
378 // Brand strings include the maximum configured frequency. These prefixes are
379 // defined by Intel CPUID documentation.
380 const char* prefixes[3] = {"MHz", "GHz", "THz"};
381 const double multipliers[3] = {1E6, 1E9, 1E12};
382 for (size_t i = 0; i < 3; ++i) {
383 const size_t pos_prefix = brand_string.find(prefixes[i]);
384 if (pos_prefix != std::string::npos) {
385 const size_t pos_space = brand_string.rfind(' ', pos_prefix - 1);
386 if (pos_space != std::string::npos) {
387 const std::string digits =
388 brand_string.substr(pos_space + 1, pos_prefix - pos_space - 1);
389 return std::stod(digits) * multipliers[i];
390 }
391 }
392 }
393
394 return 0.0;
395 }
396
397 #endif // HWY_ARCH_X86
398
399 } // namespace
400
InvariantTicksPerSecond()401 double InvariantTicksPerSecond() {
402 #if HWY_ARCH_PPC
403 return __ppc_get_timebase_freq();
404 #elif HWY_ARCH_X86
405 // We assume the TSC is invariant; it is on all recent Intel/AMD CPUs.
406 return NominalClockRate();
407 #elif defined(_WIN32) || defined(_WIN64)
408 LARGE_INTEGER freq;
409 (void)QueryPerformanceFrequency(&freq);
410 return double(freq.QuadPart);
411 #elif defined(__MACH__)
412 // https://developer.apple.com/library/mac/qa/qa1398/_index.html
413 mach_timebase_info_data_t timebase;
414 (void)mach_timebase_info(&timebase);
415 return double(timebase.denom) / timebase.numer * 1E9;
416 #else
417 // TODO(janwas): ARM? Unclear how to reliably query cntvct_el0 frequency.
418 return 1E9; // Haiku and clock_gettime return nanoseconds.
419 #endif
420 }
421
Now()422 double Now() {
423 static const double mul = 1.0 / InvariantTicksPerSecond();
424 return static_cast<double>(timer::Start()) * mul;
425 }
426
TimerResolution()427 uint64_t TimerResolution() {
428 // Nested loop avoids exceeding stack/L1 capacity.
429 timer::Ticks repetitions[Params::kTimerSamples];
430 for (size_t rep = 0; rep < Params::kTimerSamples; ++rep) {
431 timer::Ticks samples[Params::kTimerSamples];
432 for (size_t i = 0; i < Params::kTimerSamples; ++i) {
433 const timer::Ticks t0 = timer::Start();
434 const timer::Ticks t1 = timer::Stop();
435 samples[i] = t1 - t0;
436 }
437 repetitions[rep] = robust_statistics::Mode(samples);
438 }
439 return robust_statistics::Mode(repetitions);
440 }
441
442 } // namespace platform
443 namespace {
444
445 static const timer::Ticks timer_resolution = platform::TimerResolution();
446
447 // Estimates the expected value of "lambda" values with a variable number of
448 // samples until the variability "rel_mad" is less than "max_rel_mad".
449 template <class Lambda>
SampleUntilStable(const double max_rel_mad,double * rel_mad,const Params & p,const Lambda & lambda)450 timer::Ticks SampleUntilStable(const double max_rel_mad, double* rel_mad,
451 const Params& p, const Lambda& lambda) {
452 // Choose initial samples_per_eval based on a single estimated duration.
453 timer::Ticks t0 = timer::Start();
454 lambda();
455 timer::Ticks t1 = timer::Stop();
456 timer::Ticks est = t1 - t0;
457 static const double ticks_per_second = platform::InvariantTicksPerSecond();
458 const size_t ticks_per_eval =
459 static_cast<size_t>(ticks_per_second * p.seconds_per_eval);
460 size_t samples_per_eval =
461 est == 0 ? p.min_samples_per_eval : ticks_per_eval / est;
462 samples_per_eval = std::max(samples_per_eval, p.min_samples_per_eval);
463
464 std::vector<timer::Ticks> samples;
465 samples.reserve(1 + samples_per_eval);
466 samples.push_back(est);
467
468 // Percentage is too strict for tiny differences, so also allow a small
469 // absolute "median absolute deviation".
470 const timer::Ticks max_abs_mad = (timer_resolution + 99) / 100;
471 *rel_mad = 0.0; // ensure initialized
472
473 for (size_t eval = 0; eval < p.max_evals; ++eval, samples_per_eval *= 2) {
474 samples.reserve(samples.size() + samples_per_eval);
475 for (size_t i = 0; i < samples_per_eval; ++i) {
476 t0 = timer::Start();
477 lambda();
478 t1 = timer::Stop();
479 samples.push_back(t1 - t0);
480 }
481
482 if (samples.size() >= p.min_mode_samples) {
483 est = robust_statistics::Mode(samples.data(), samples.size());
484 } else {
485 // For "few" (depends also on the variance) samples, Median is safer.
486 est = robust_statistics::Median(samples.data(), samples.size());
487 }
488 NANOBENCHMARK_CHECK(est != 0);
489
490 // Median absolute deviation (mad) is a robust measure of 'variability'.
491 const timer::Ticks abs_mad = robust_statistics::MedianAbsoluteDeviation(
492 samples.data(), samples.size(), est);
493 *rel_mad = static_cast<double>(abs_mad) / static_cast<double>(est);
494
495 if (*rel_mad <= max_rel_mad || abs_mad <= max_abs_mad) {
496 if (p.verbose) {
497 printf("%6zu samples => %5zu (abs_mad=%4zu, rel_mad=%4.2f%%)\n",
498 samples.size(), size_t(est), size_t(abs_mad), *rel_mad * 100.0);
499 }
500 return est;
501 }
502 }
503
504 if (p.verbose) {
505 printf(
506 "WARNING: rel_mad=%4.2f%% still exceeds %4.2f%% after %6zu samples.\n",
507 *rel_mad * 100.0, max_rel_mad * 100.0, samples.size());
508 }
509 return est;
510 }
511
512 using InputVec = std::vector<FuncInput>;
513
514 // Returns vector of unique input values.
UniqueInputs(const FuncInput * inputs,const size_t num_inputs)515 InputVec UniqueInputs(const FuncInput* inputs, const size_t num_inputs) {
516 InputVec unique(inputs, inputs + num_inputs);
517 std::sort(unique.begin(), unique.end());
518 unique.erase(std::unique(unique.begin(), unique.end()), unique.end());
519 return unique;
520 }
521
522 // Returns how often we need to call func for sufficient precision.
NumSkip(const Func func,const uint8_t * arg,const InputVec & unique,const Params & p)523 size_t NumSkip(const Func func, const uint8_t* arg, const InputVec& unique,
524 const Params& p) {
525 // Min elapsed ticks for any input.
526 timer::Ticks min_duration = ~timer::Ticks(0);
527
528 for (const FuncInput input : unique) {
529 double rel_mad;
530 const timer::Ticks total = SampleUntilStable(
531 p.target_rel_mad, &rel_mad, p,
532 [func, arg, input]() { platform::PreventElision(func(arg, input)); });
533 min_duration = std::min(min_duration, total - timer_resolution);
534 }
535
536 // Number of repetitions required to reach the target resolution.
537 const size_t max_skip = p.precision_divisor;
538 // Number of repetitions given the estimated duration.
539 const size_t num_skip =
540 min_duration == 0 ? 0 : (max_skip + min_duration - 1) / min_duration;
541 if (p.verbose) {
542 printf("res=%zu max_skip=%zu min_dur=%zu num_skip=%zu\n",
543 size_t(timer_resolution), max_skip, size_t(min_duration), num_skip);
544 }
545 return num_skip;
546 }
547
548 // Replicates inputs until we can omit "num_skip" occurrences of an input.
ReplicateInputs(const FuncInput * inputs,const size_t num_inputs,const size_t num_unique,const size_t num_skip,const Params & p)549 InputVec ReplicateInputs(const FuncInput* inputs, const size_t num_inputs,
550 const size_t num_unique, const size_t num_skip,
551 const Params& p) {
552 InputVec full;
553 if (num_unique == 1) {
554 full.assign(p.subset_ratio * num_skip, inputs[0]);
555 return full;
556 }
557
558 full.reserve(p.subset_ratio * num_skip * num_inputs);
559 for (size_t i = 0; i < p.subset_ratio * num_skip; ++i) {
560 full.insert(full.end(), inputs, inputs + num_inputs);
561 }
562 std::mt19937 rng;
563 std::shuffle(full.begin(), full.end(), rng);
564 return full;
565 }
566
567 // Copies the "full" to "subset" in the same order, but with "num_skip"
568 // randomly selected occurrences of "input_to_skip" removed.
FillSubset(const InputVec & full,const FuncInput input_to_skip,const size_t num_skip,InputVec * subset)569 void FillSubset(const InputVec& full, const FuncInput input_to_skip,
570 const size_t num_skip, InputVec* subset) {
571 const size_t count =
572 static_cast<size_t>(std::count(full.begin(), full.end(), input_to_skip));
573 // Generate num_skip random indices: which occurrence to skip.
574 std::vector<uint32_t> omit(count);
575 std::iota(omit.begin(), omit.end(), 0);
576 // omit[] is the same on every call, but that's OK because they identify the
577 // Nth instance of input_to_skip, so the position within full[] differs.
578 std::mt19937 rng;
579 std::shuffle(omit.begin(), omit.end(), rng);
580 omit.resize(num_skip);
581 std::sort(omit.begin(), omit.end());
582
583 uint32_t occurrence = ~0u; // 0 after preincrement
584 size_t idx_omit = 0; // cursor within omit[]
585 size_t idx_subset = 0; // cursor within *subset
586 for (const FuncInput next : full) {
587 if (next == input_to_skip) {
588 ++occurrence;
589 // Haven't removed enough already
590 if (idx_omit < num_skip) {
591 // This one is up for removal
592 if (occurrence == omit[idx_omit]) {
593 ++idx_omit;
594 continue;
595 }
596 }
597 }
598 if (idx_subset < subset->size()) {
599 (*subset)[idx_subset++] = next;
600 }
601 }
602 NANOBENCHMARK_CHECK(idx_subset == subset->size());
603 NANOBENCHMARK_CHECK(idx_omit == omit.size());
604 NANOBENCHMARK_CHECK(occurrence == count - 1);
605 }
606
607 // Returns total ticks elapsed for all inputs.
TotalDuration(const Func func,const uint8_t * arg,const InputVec * inputs,const Params & p,double * max_rel_mad)608 timer::Ticks TotalDuration(const Func func, const uint8_t* arg,
609 const InputVec* inputs, const Params& p,
610 double* max_rel_mad) {
611 double rel_mad;
612 const timer::Ticks duration =
613 SampleUntilStable(p.target_rel_mad, &rel_mad, p, [func, arg, inputs]() {
614 for (const FuncInput input : *inputs) {
615 platform::PreventElision(func(arg, input));
616 }
617 });
618 *max_rel_mad = std::max(*max_rel_mad, rel_mad);
619 return duration;
620 }
621
622 // (Nearly) empty Func for measuring timer overhead/resolution.
EmptyFunc(const void *,const FuncInput input)623 HWY_NOINLINE FuncOutput EmptyFunc(const void* /*arg*/, const FuncInput input) {
624 return input;
625 }
626
627 // Returns overhead of accessing inputs[] and calling a function; this will
628 // be deducted from future TotalDuration return values.
Overhead(const uint8_t * arg,const InputVec * inputs,const Params & p)629 timer::Ticks Overhead(const uint8_t* arg, const InputVec* inputs,
630 const Params& p) {
631 double rel_mad;
632 // Zero tolerance because repeatability is crucial and EmptyFunc is fast.
633 return SampleUntilStable(0.0, &rel_mad, p, [arg, inputs]() {
634 for (const FuncInput input : *inputs) {
635 platform::PreventElision(EmptyFunc(arg, input));
636 }
637 });
638 }
639
640 } // namespace
641
Unpredictable1()642 int Unpredictable1() { return timer::Start() != ~0ULL; }
643
Measure(const Func func,const uint8_t * arg,const FuncInput * inputs,const size_t num_inputs,Result * results,const Params & p)644 size_t Measure(const Func func, const uint8_t* arg, const FuncInput* inputs,
645 const size_t num_inputs, Result* results, const Params& p) {
646 NANOBENCHMARK_CHECK(num_inputs != 0);
647 const InputVec& unique = UniqueInputs(inputs, num_inputs);
648
649 const size_t num_skip = NumSkip(func, arg, unique, p); // never 0
650 if (num_skip == 0) return 0; // NumSkip already printed error message
651 // (slightly less work on x86 to cast from signed integer)
652 const float mul = 1.0f / static_cast<float>(static_cast<int>(num_skip));
653
654 const InputVec& full =
655 ReplicateInputs(inputs, num_inputs, unique.size(), num_skip, p);
656 InputVec subset(full.size() - num_skip);
657
658 const timer::Ticks overhead = Overhead(arg, &full, p);
659 const timer::Ticks overhead_skip = Overhead(arg, &subset, p);
660 if (overhead < overhead_skip) {
661 fprintf(stderr, "Measurement failed: overhead %zu < %zu\n",
662 size_t(overhead), size_t(overhead_skip));
663 return 0;
664 }
665
666 if (p.verbose) {
667 printf("#inputs=%5zu,%5zu overhead=%5zu,%5zu\n", full.size(), subset.size(),
668 size_t(overhead), size_t(overhead_skip));
669 }
670
671 double max_rel_mad = 0.0;
672 const timer::Ticks total = TotalDuration(func, arg, &full, p, &max_rel_mad);
673
674 for (size_t i = 0; i < unique.size(); ++i) {
675 FillSubset(full, unique[i], num_skip, &subset);
676 const timer::Ticks total_skip =
677 TotalDuration(func, arg, &subset, p, &max_rel_mad);
678
679 if (total < total_skip) {
680 fprintf(stderr, "Measurement failed: total %zu < %zu\n", size_t(total),
681 size_t(total_skip));
682 return 0;
683 }
684
685 const timer::Ticks duration =
686 (total - overhead) - (total_skip - overhead_skip);
687 results[i].input = unique[i];
688 results[i].ticks = static_cast<float>(duration) * mul;
689 results[i].variability = static_cast<float>(max_rel_mad);
690 }
691
692 return unique.size();
693 }
694
695 } // namespace hwy
696