third_party/highway/hwy/nanobenchmark.cc

   1 // Copyright 2019 Google LLC
   2 // SPDX-License-Identifier: Apache-2.0
   3 //
   4 // Licensed under the Apache License, Version 2.0 (the "License");
   5 // you may not use this file except in compliance with the License.
   6 // You may obtain a copy of the License at
   7 //
   8 //      http://www.apache.org/licenses/LICENSE-2.0
   9 //
  10 // Unless required by applicable law or agreed to in writing, software
  11 // distributed under the License is distributed on an "AS IS" BASIS,
  12 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  13 // See the License for the specific language governing permissions and
  14 // limitations under the License.
  15
  16 #include "hwy/nanobenchmark.h"
  17
  18 #include <stdio.h>
  19 #include <stdlib.h>
  20 #include <time.h>  // clock_gettime
  21
  22 #include <algorithm>  // std::sort, std::find_if
  23 #include <numeric>    // std::iota
  24 #include <random>
  25 #include <vector>
  26
  27 #include "hwy/robust_statistics.h"
  28 #include "hwy/timer-inl.h"
  29 #include "hwy/timer.h"
  30
  31 namespace hwy {
  32 namespace {
  33 namespace timer = hwy::HWY_NAMESPACE::timer;
  34
  35 static const timer::Ticks timer_resolution = platform::TimerResolution();
  36
  37 // Estimates the expected value of "lambda" values with a variable number of
  38 // samples until the variability "rel_mad" is less than "max_rel_mad".
  39 template <class Lambda>
  40 timer::Ticks SampleUntilStable(const double max_rel_mad, double* rel_mad,
  41                                const Params& p, const Lambda& lambda) {
  42   // Choose initial samples_per_eval based on a single estimated duration.
  43   timer::Ticks t0 = timer::Start();
  44   lambda();
  45   timer::Ticks t1 = timer::Stop();  // Caller checks HaveTimerStop
  46   timer::Ticks est = t1 - t0;
  47   static const double ticks_per_second = platform::InvariantTicksPerSecond();
  48   const size_t ticks_per_eval =
  49       static_cast<size_t>(ticks_per_second * p.seconds_per_eval);
  50   size_t samples_per_eval = est == 0
  51                                 ? p.min_samples_per_eval
  52                                 : static_cast<size_t>(ticks_per_eval / est);
  53   samples_per_eval = HWY_MAX(samples_per_eval, p.min_samples_per_eval);
  54
  55   std::vector<timer::Ticks> samples;
  56   samples.reserve(1 + samples_per_eval);
  57   samples.push_back(est);
  58
  59   // Percentage is too strict for tiny differences, so also allow a small
  60   // absolute "median absolute deviation".
  61   const timer::Ticks max_abs_mad = (timer_resolution + 99) / 100;
  62   *rel_mad = 0.0;  // ensure initialized
  63
  64   for (size_t eval = 0; eval < p.max_evals; ++eval, samples_per_eval *= 2) {
  65     samples.reserve(samples.size() + samples_per_eval);
  66     for (size_t i = 0; i < samples_per_eval; ++i) {
  67       t0 = timer::Start();
  68       lambda();
  69       t1 = timer::Stop();  // Caller checks HaveTimerStop
  70       samples.push_back(t1 - t0);
  71     }
  72
  73     if (samples.size() >= p.min_mode_samples) {
  74       est = robust_statistics::Mode(samples.data(), samples.size());
  75     } else {
  76       // For "few" (depends also on the variance) samples, Median is safer.
  77       est = robust_statistics::Median(samples.data(), samples.size());
  78     }
  79     NANOBENCHMARK_CHECK(est != 0);
  80
  81     // Median absolute deviation (mad) is a robust measure of 'variability'.
  82     const timer::Ticks abs_mad = robust_statistics::MedianAbsoluteDeviation(
  83         samples.data(), samples.size(), est);
  84     *rel_mad = static_cast<double>(abs_mad) / static_cast<double>(est);
  85
  86     if (*rel_mad <= max_rel_mad || abs_mad <= max_abs_mad) {
  87       if (p.verbose) {
  88         printf("%6d samples => %5d (abs_mad=%4d, rel_mad=%4.2f%%)\n",
  89                static_cast<int>(samples.size()), static_cast<int>(est),
  90                static_cast<int>(abs_mad), *rel_mad * 100.0);
  91       }
  92       return est;
  93     }
  94   }
  95
  96   if (p.verbose) {
  97     printf("WARNING: rel_mad=%4.2f%% still exceeds %4.2f%% after %6d samples\n",
  98            *rel_mad * 100.0, max_rel_mad * 100.0,
  99            static_cast<int>(samples.size()));
 100   }
 101   return est;
 102 }
 103
 104 using InputVec = std::vector<FuncInput>;
 105
 106 // Returns vector of unique input values.
 107 InputVec UniqueInputs(const FuncInput* inputs, const size_t num_inputs) {
 108   InputVec unique(inputs, inputs + num_inputs);
 109   std::sort(unique.begin(), unique.end());
 110   unique.erase(std::unique(unique.begin(), unique.end()), unique.end());
 111   return unique;
 112 }
 113
 114 // Returns how often we need to call func for sufficient precision.
 115 size_t NumSkip(const Func func, const uint8_t* arg, const InputVec& unique,
 116                const Params& p) {
 117   // Min elapsed ticks for any input.
 118   timer::Ticks min_duration = ~timer::Ticks(0);
 119
 120   for (const FuncInput input : unique) {
 121     double rel_mad;
 122     const timer::Ticks total = SampleUntilStable(
 123         p.target_rel_mad, &rel_mad, p,
 124         [func, arg, input]() { PreventElision(func(arg, input)); });
 125     min_duration = HWY_MIN(min_duration, total - timer_resolution);
 126   }
 127
 128   // Number of repetitions required to reach the target resolution.
 129   const size_t max_skip = p.precision_divisor;
 130   // Number of repetitions given the estimated duration.
 131   const size_t num_skip =
 132       min_duration == 0
 133           ? 0
 134           : static_cast<size_t>((max_skip + min_duration - 1) / min_duration);
 135   if (p.verbose) {
 136     printf("res=%d max_skip=%d min_dur=%d num_skip=%d\n",
 137            static_cast<int>(timer_resolution), static_cast<int>(max_skip),
 138            static_cast<int>(min_duration), static_cast<int>(num_skip));
 139   }
 140   return num_skip;
 141 }
 142
 143 // Replicates inputs until we can omit "num_skip" occurrences of an input.
 144 InputVec ReplicateInputs(const FuncInput* inputs, const size_t num_inputs,
 145                          const size_t num_unique, const size_t num_skip,
 146                          const Params& p) {
 147   InputVec full;
 148   if (num_unique == 1) {
 149     full.assign(p.subset_ratio * num_skip, inputs[0]);
 150     return full;
 151   }
 152
 153   full.reserve(p.subset_ratio * num_skip * num_inputs);
 154   for (size_t i = 0; i < p.subset_ratio * num_skip; ++i) {
 155     full.insert(full.end(), inputs, inputs + num_inputs);
 156   }
 157   std::mt19937 rng;
 158   std::shuffle(full.begin(), full.end(), rng);
 159   return full;
 160 }
 161
 162 // Copies the "full" to "subset" in the same order, but with "num_skip"
 163 // randomly selected occurrences of "input_to_skip" removed.
 164 void FillSubset(const InputVec& full, const FuncInput input_to_skip,
 165                 const size_t num_skip, InputVec* subset) {
 166   const size_t count =
 167       static_cast<size_t>(std::count(full.begin(), full.end(), input_to_skip));
 168   // Generate num_skip random indices: which occurrence to skip.
 169   std::vector<uint32_t> omit(count);
 170   std::iota(omit.begin(), omit.end(), 0);
 171   // omit[] is the same on every call, but that's OK because they identify the
 172   // Nth instance of input_to_skip, so the position within full[] differs.
 173   std::mt19937 rng;
 174   std::shuffle(omit.begin(), omit.end(), rng);
 175   omit.resize(num_skip);
 176   std::sort(omit.begin(), omit.end());
 177
 178   uint32_t occurrence = ~0u;  // 0 after preincrement
 179   size_t idx_omit = 0;        // cursor within omit[]
 180   size_t idx_subset = 0;      // cursor within *subset
 181   for (const FuncInput next : full) {
 182     if (next == input_to_skip) {
 183       ++occurrence;
 184       // Haven't removed enough already
 185       if (idx_omit < num_skip) {
 186         // This one is up for removal
 187         if (occurrence == omit[idx_omit]) {
 188           ++idx_omit;
 189           continue;
 190         }
 191       }
 192     }
 193     if (idx_subset < subset->size()) {
 194       (*subset)[idx_subset++] = next;
 195     }
 196   }
 197   NANOBENCHMARK_CHECK(idx_subset == subset->size());
 198   NANOBENCHMARK_CHECK(idx_omit == omit.size());
 199   NANOBENCHMARK_CHECK(occurrence == count - 1);
 200 }
 201
 202 // Returns total ticks elapsed for all inputs.
 203 timer::Ticks TotalDuration(const Func func, const uint8_t* arg,
 204                            const InputVec* inputs, const Params& p,
 205                            double* max_rel_mad) {
 206   double rel_mad;
 207   const timer::Ticks duration =
 208       SampleUntilStable(p.target_rel_mad, &rel_mad, p, [func, arg, inputs]() {
 209         for (const FuncInput input : *inputs) {
 210           PreventElision(func(arg, input));
 211         }
 212       });
 213   *max_rel_mad = HWY_MAX(*max_rel_mad, rel_mad);
 214   return duration;
 215 }
 216
 217 // (Nearly) empty Func for measuring timer overhead/resolution.
 218 HWY_NOINLINE FuncOutput EmptyFunc(const void* /*arg*/, const FuncInput input) {
 219   return input;
 220 }
 221
 222 // Returns overhead of accessing inputs[] and calling a function; this will
 223 // be deducted from future TotalDuration return values.
 224 timer::Ticks Overhead(const uint8_t* arg, const InputVec* inputs,
 225                       const Params& p) {
 226   double rel_mad;
 227   // Zero tolerance because repeatability is crucial and EmptyFunc is fast.
 228   return SampleUntilStable(0.0, &rel_mad, p, [arg, inputs]() {
 229     for (const FuncInput input : *inputs) {
 230       PreventElision(EmptyFunc(arg, input));
 231     }
 232   });
 233 }
 234
 235 }  // namespace
 236
 237 HWY_DLLEXPORT int Unpredictable1() { return timer::Start() != ~0ULL; }
 238
 239 HWY_DLLEXPORT size_t Measure(const Func func, const uint8_t* arg,
 240                              const FuncInput* inputs, const size_t num_inputs,
 241                              Result* results, const Params& p) {
 242   NANOBENCHMARK_CHECK(num_inputs != 0);
 243
 244   char cpu100[100];
 245   if (!platform::HaveTimerStop(cpu100)) {
 246     fprintf(stderr, "CPU '%s' does not support RDTSCP, skipping benchmark.\n",
 247             cpu100);
 248     return 0;
 249   }
 250
 251   const InputVec& unique = UniqueInputs(inputs, num_inputs);
 252
 253   const size_t num_skip = NumSkip(func, arg, unique, p);  // never 0
 254   if (num_skip == 0) return 0;  // NumSkip already printed error message
 255   // (slightly less work on x86 to cast from signed integer)
 256   const float mul = 1.0f / static_cast<float>(static_cast<int>(num_skip));
 257
 258   const InputVec& full =
 259       ReplicateInputs(inputs, num_inputs, unique.size(), num_skip, p);
 260   InputVec subset(full.size() - num_skip);
 261
 262   const timer::Ticks overhead = Overhead(arg, &full, p);
 263   const timer::Ticks overhead_skip = Overhead(arg, &subset, p);
 264   if (overhead < overhead_skip) {
 265     fprintf(stderr, "Measurement failed: overhead %d < %d\n",
 266             static_cast<int>(overhead), static_cast<int>(overhead_skip));
 267     return 0;
 268   }
 269
 270   if (p.verbose) {
 271     printf("#inputs=%5d,%5d overhead=%5d,%5d\n", static_cast<int>(full.size()),
 272            static_cast<int>(subset.size()), static_cast<int>(overhead),
 273            static_cast<int>(overhead_skip));
 274   }
 275
 276   double max_rel_mad = 0.0;
 277   const timer::Ticks total = TotalDuration(func, arg, &full, p, &max_rel_mad);
 278
 279   for (size_t i = 0; i < unique.size(); ++i) {
 280     FillSubset(full, unique[i], num_skip, &subset);
 281     const timer::Ticks total_skip =
 282         TotalDuration(func, arg, &subset, p, &max_rel_mad);
 283
 284     if (total < total_skip) {
 285       fprintf(stderr, "Measurement failed: total %f < %f\n",
 286               static_cast<double>(total), static_cast<double>(total_skip));
 287       return 0;
 288     }
 289
 290     const timer::Ticks duration =
 291         (total - overhead) - (total_skip - overhead_skip);
 292     results[i].input = unique[i];
 293     results[i].ticks = static_cast<float>(duration) * mul;
 294     results[i].variability = static_cast<float>(max_rel_mad);
 295   }
 296
 297   return unique.size();
 298 }
 299
 300 }  // namespace hwy