|
| 1 | +#include "../src/detail/bruun_dif_kernel.hpp" |
| 2 | +#include "../src/detail/bruun_dip_kernel.hpp" |
| 3 | + |
| 4 | +#include <algorithm> |
| 5 | +#include <chrono> |
| 6 | +#include <cmath> |
| 7 | +#include <cstdio> |
| 8 | +#include <cstdlib> |
| 9 | +#include <exception> |
| 10 | +#include <random> |
| 11 | +#include <stdexcept> |
| 12 | +#include <vector> |
| 13 | + |
| 14 | +namespace { |
| 15 | + |
| 16 | +using clock_type = std::chrono::steady_clock; |
| 17 | + |
| 18 | +constexpr double pi = 3.141592653589793238462643383279502884; |
| 19 | + |
| 20 | +struct timing { |
| 21 | + double best_ns = 0.0; |
| 22 | + double sink = 0.0; |
| 23 | +}; |
| 24 | + |
| 25 | +struct result { |
| 26 | + std::size_t n = 0; |
| 27 | + int iters = 0; |
| 28 | + double dif_forward_ns = 0.0; |
| 29 | + double dip_forward_ns = 0.0; |
| 30 | + double dif_inverse_ns = 0.0; |
| 31 | + double dip_inverse_ns = 0.0; |
| 32 | + double dif_roundtrip_ns = 0.0; |
| 33 | + double dip_roundtrip_ns = 0.0; |
| 34 | + double forward_maxerr = 0.0; |
| 35 | + double dif_roundtrip_maxerr = 0.0; |
| 36 | + double dip_roundtrip_maxerr = 0.0; |
| 37 | + double sink = 0.0; |
| 38 | +}; |
| 39 | + |
| 40 | +bool is_power2(std::size_t n) { |
| 41 | + return n > 0 && (n & (n - 1)) == 0; |
| 42 | +} |
| 43 | + |
| 44 | +std::size_t parse_size(const char* text) { |
| 45 | + char* end = nullptr; |
| 46 | + const unsigned long long value = std::strtoull(text, &end, 0); |
| 47 | + if (!end || *end != '\0' || value == 0) { |
| 48 | + throw std::invalid_argument("invalid size"); |
| 49 | + } |
| 50 | + return static_cast<std::size_t>(value); |
| 51 | +} |
| 52 | + |
| 53 | +int parse_iters(const char* text) { |
| 54 | + char* end = nullptr; |
| 55 | + const long value = std::strtol(text, &end, 0); |
| 56 | + if (!end || *end != '\0' || value <= 0) { |
| 57 | + throw std::invalid_argument("invalid iteration count"); |
| 58 | + } |
| 59 | + return static_cast<int>(value); |
| 60 | +} |
| 61 | + |
| 62 | +int default_iters(std::size_t n) { |
| 63 | + const std::size_t target_samples = 1u << 20; |
| 64 | + std::size_t iters = target_samples / n; |
| 65 | + if (iters < 3) iters = 3; |
| 66 | + if (iters > 4096) iters = 4096; |
| 67 | + return static_cast<int>(iters); |
| 68 | +} |
| 69 | + |
| 70 | +std::vector<double> make_signal(std::size_t n) { |
| 71 | + std::vector<double> input(n); |
| 72 | + std::mt19937_64 rng(0xD1F0D1F0ULL + static_cast<unsigned long long>(n)); |
| 73 | + std::uniform_real_distribution<double> noise(-0.05, 0.05); |
| 74 | + for (std::size_t i = 0; i < n; ++i) { |
| 75 | + const double t = static_cast<double>(i) / static_cast<double>(n); |
| 76 | + input[i] = std::sin(2.0 * pi * 13.0 * t) |
| 77 | + + 0.5 * std::cos(2.0 * pi * 37.0 * t) |
| 78 | + + 0.25 * std::sin(2.0 * pi * 89.0 * t) |
| 79 | + + noise(rng); |
| 80 | + } |
| 81 | + return input; |
| 82 | +} |
| 83 | + |
| 84 | +template <typename Func> |
| 85 | +timing bench_ns(int iters, Func&& func) { |
| 86 | + timing best; |
| 87 | + best.best_ns = 1.0e300; |
| 88 | + for (int pass = 0; pass < 7; ++pass) { |
| 89 | + double sink = 0.0; |
| 90 | + const auto start = clock_type::now(); |
| 91 | + for (int i = 0; i < iters; ++i) { |
| 92 | + sink += func(i, sink); |
| 93 | + } |
| 94 | + const auto stop = clock_type::now(); |
| 95 | + const double ns = std::chrono::duration<double, std::nano>(stop - start).count() |
| 96 | + / static_cast<double>(iters); |
| 97 | + if (ns < best.best_ns) { |
| 98 | + best.best_ns = ns; |
| 99 | + best.sink = sink; |
| 100 | + } |
| 101 | + } |
| 102 | + return best; |
| 103 | +} |
| 104 | + |
| 105 | +template <typename A, typename B> |
| 106 | +double max_abs_complex( const std::vector<A>& a, const std::vector<B>& b) { |
| 107 | + double err = 0.0; |
| 108 | + for (std::size_t i = 0; i < a.size(); ++i) { |
| 109 | + err = std::max(err, std::abs(a[i].re - b[i].re)); |
| 110 | + err = std::max(err, std::abs(a[i].im - b[i].im)); |
| 111 | + } |
| 112 | + return err; |
| 113 | +} |
| 114 | + |
| 115 | +double max_abs_real(const std::vector<double>& a, const std::vector<double>& b) { |
| 116 | + double err = 0.0; |
| 117 | + for (std::size_t i = 0; i < a.size(); ++i) { |
| 118 | + err = std::max(err, std::abs(a[i] - b[i])); |
| 119 | + } |
| 120 | + return err; |
| 121 | +} |
| 122 | + |
| 123 | +result run_one(std::size_t n, int forced_iters) { |
| 124 | + if (n > static_cast<std::size_t>(2147483647)) { |
| 125 | + throw std::invalid_argument("N is too large for the experimental kernels"); |
| 126 | + } |
| 127 | + |
| 128 | + bruun::DIF_RFFT_kernel dif; |
| 129 | + bruun::DIP_RFFT_kernel dip; |
| 130 | + if (!dif.init(static_cast<int>(n)) || !dip.init(static_cast<int>(n))) { |
| 131 | + throw std::runtime_error("kernel setup failed"); |
| 132 | + } |
| 133 | + |
| 134 | + const std::size_t nb = n / 2 + 1; |
| 135 | + const int iters = forced_iters > 0 ? forced_iters : default_iters(n); |
| 136 | + const std::vector<double> original = make_signal(n); |
| 137 | + std::vector<double> input(original); |
| 138 | + |
| 139 | + std::vector<bruun::complex_t> dif_bins(nb); |
| 140 | + std::vector<bruun::complex_t> dip_bins(nb); |
| 141 | + std::vector<bruun::complex_t> dif_scratch(static_cast<std::size_t>(dif.native_scratch_size())); |
| 142 | + std::vector<double> dif_work(static_cast<std::size_t>(dif.work_size())); |
| 143 | + std::vector<double> dip_work(static_cast<std::size_t>(dip.work_size())); |
| 144 | + std::vector<double> dif_out(n); |
| 145 | + std::vector<double> dip_out(n); |
| 146 | + |
| 147 | + dif.forward_standard(original.data(), dif_bins.data(), dif_work.data(), dif_scratch.data()); |
| 148 | + dip.forward_standard(original.data(), dip_bins.data(), dip_work.data()); |
| 149 | + dif.inverse(dif_bins.data(), dif_out.data()); |
| 150 | + dip.inverse_standard(dip_bins.data(), dip_out.data(), dip_work.data()); |
| 151 | + |
| 152 | + result r; |
| 153 | + r.n = n; |
| 154 | + r.iters = iters; |
| 155 | + r.forward_maxerr = max_abs_complex(dif_bins, dip_bins); |
| 156 | + r.dif_roundtrip_maxerr = max_abs_real(original, dif_out); |
| 157 | + r.dip_roundtrip_maxerr = max_abs_real(original, dip_out); |
| 158 | + |
| 159 | + input = original; |
| 160 | + timing t = bench_ns(iters, [&](int i, double sink) { |
| 161 | + input[(static_cast<std::size_t>(i) * 131u + static_cast<std::size_t>(sink)) & (n - 1)] += 1e-12; |
| 162 | + dif.forward_standard(input.data(), dif_bins.data(), dif_work.data(), dif_scratch.data()); |
| 163 | + return dif_bins[(static_cast<std::size_t>(i) * 17u) % nb].re; |
| 164 | + }); |
| 165 | + r.dif_forward_ns = t.best_ns; |
| 166 | + r.sink += t.sink; |
| 167 | + |
| 168 | + input = original; |
| 169 | + t = bench_ns(iters, [&](int i, double sink) { |
| 170 | + input[(static_cast<std::size_t>(i) * 131u + static_cast<std::size_t>(sink)) & (n - 1)] += 1e-12; |
| 171 | + dip.forward_standard(input.data(), dip_bins.data(), dip_work.data()); |
| 172 | + return dip_bins[(static_cast<std::size_t>(i) * 17u) % nb].re; |
| 173 | + }); |
| 174 | + r.dip_forward_ns = t.best_ns; |
| 175 | + r.sink += t.sink; |
| 176 | + |
| 177 | + dif.forward_standard(original.data(), dif_bins.data(), dif_work.data(), dif_scratch.data()); |
| 178 | + t = bench_ns(iters, [&](int i, double) { |
| 179 | + dif_bins[(static_cast<std::size_t>(i) * 17u) % nb].re += 1e-12; |
| 180 | + dif.inverse(dif_bins.data(), dif_out.data()); |
| 181 | + return dif_out[(static_cast<std::size_t>(i) * 31u) & (n - 1)]; |
| 182 | + }); |
| 183 | + r.dif_inverse_ns = t.best_ns; |
| 184 | + r.sink += t.sink; |
| 185 | + |
| 186 | + dip.forward_standard(original.data(), dip_bins.data(), dip_work.data()); |
| 187 | + t = bench_ns(iters, [&](int i, double) { |
| 188 | + dip_bins[(static_cast<std::size_t>(i) * 17u) % nb].re += 1e-12; |
| 189 | + dip.inverse_standard(dip_bins.data(), dip_out.data(), dip_work.data()); |
| 190 | + return dip_out[(static_cast<std::size_t>(i) * 31u) & (n - 1)]; |
| 191 | + }); |
| 192 | + r.dip_inverse_ns = t.best_ns; |
| 193 | + r.sink += t.sink; |
| 194 | + |
| 195 | + input = original; |
| 196 | + t = bench_ns(iters, [&](int i, double sink) { |
| 197 | + input[(static_cast<std::size_t>(i) * 131u + static_cast<std::size_t>(sink)) & (n - 1)] += 1e-12; |
| 198 | + dif.forward_standard(input.data(), dif_bins.data(), dif_work.data(), dif_scratch.data()); |
| 199 | + dif.inverse(dif_bins.data(), dif_out.data()); |
| 200 | + return dif_out[(static_cast<std::size_t>(i) * 31u) & (n - 1)]; |
| 201 | + }); |
| 202 | + r.dif_roundtrip_ns = t.best_ns; |
| 203 | + r.sink += t.sink; |
| 204 | + |
| 205 | + input = original; |
| 206 | + t = bench_ns(iters, [&](int i, double sink) { |
| 207 | + input[(static_cast<std::size_t>(i) * 131u + static_cast<std::size_t>(sink)) & (n - 1)] += 1e-12; |
| 208 | + dip.forward_standard(input.data(), dip_bins.data(), dip_work.data()); |
| 209 | + dip.inverse_standard(dip_bins.data(), dip_out.data(), dip_work.data()); |
| 210 | + return dip_out[(static_cast<std::size_t>(i) * 31u) & (n - 1)]; |
| 211 | + }); |
| 212 | + r.dip_roundtrip_ns = t.best_ns; |
| 213 | + r.sink += t.sink; |
| 214 | + |
| 215 | + return r; |
| 216 | +} |
| 217 | + |
| 218 | +void print_header() { |
| 219 | + std::printf("%9s %8s %13s %13s %13s %13s %13s %13s %8s %8s %8s %12s\n", |
| 220 | + "N", "iters", "DIF_fwd_ns", "DIP_fwd_ns", "DIF_inv_ns", |
| 221 | + "DIP_inv_ns", "DIF_rt_ns", "DIP_rt_ns", "fwd_x", |
| 222 | + "inv_x", "rt_x", "checks"); |
| 223 | + std::fflush(stdout); |
| 224 | +} |
| 225 | + |
| 226 | +void print_result(const result& r) { |
| 227 | + const double fwd_ratio = r.dip_forward_ns / r.dif_forward_ns; |
| 228 | + const double inv_ratio = r.dip_inverse_ns / r.dif_inverse_ns; |
| 229 | + const double rt_ratio = r.dip_roundtrip_ns / r.dif_roundtrip_ns; |
| 230 | + std::printf("%9zu %8d %13.2f %13.2f %13.2f %13.2f %13.2f %13.2f %8.3f %8.3f %8.3f " |
| 231 | + "ferr %.2e dif_rt %.2e dip_rt %.2e sink %.3e\n", |
| 232 | + r.n, |
| 233 | + r.iters, |
| 234 | + r.dif_forward_ns, |
| 235 | + r.dip_forward_ns, |
| 236 | + r.dif_inverse_ns, |
| 237 | + r.dip_inverse_ns, |
| 238 | + r.dif_roundtrip_ns, |
| 239 | + r.dip_roundtrip_ns, |
| 240 | + fwd_ratio, |
| 241 | + inv_ratio, |
| 242 | + rt_ratio, |
| 243 | + r.forward_maxerr, |
| 244 | + r.dif_roundtrip_maxerr, |
| 245 | + r.dip_roundtrip_maxerr, |
| 246 | + r.sink); |
| 247 | + std::fflush(stdout); |
| 248 | +} |
| 249 | + |
| 250 | +} // namespace |
| 251 | + |
| 252 | +int main(int argc, char** argv) { |
| 253 | + try { |
| 254 | + if (argc > 3) { |
| 255 | + std::fprintf(stderr, "usage: %s [N | max_pow] [iters]\n", argv[0]); |
| 256 | + return 2; |
| 257 | + } |
| 258 | + |
| 259 | + int forced_iters = 0; |
| 260 | + if (argc >= 3) { |
| 261 | + forced_iters = parse_iters(argv[2]); |
| 262 | + } |
| 263 | + |
| 264 | + print_header(); |
| 265 | + if (argc >= 2) { |
| 266 | + const std::size_t arg = parse_size(argv[1]); |
| 267 | + if (is_power2(arg) && arg >= 4) { |
| 268 | + print_result(run_one(arg, forced_iters)); |
| 269 | + } else if (arg >= 2 && arg <= 30) { |
| 270 | + for (int p = 2; p <= static_cast<int>(arg); ++p) { |
| 271 | + print_result(run_one(static_cast<std::size_t>(1) << p, forced_iters)); |
| 272 | + } |
| 273 | + } else { |
| 274 | + throw std::invalid_argument("argument must be a power-of-two N or max_pow in [2, 30]"); |
| 275 | + } |
| 276 | + return 0; |
| 277 | + } |
| 278 | + |
| 279 | + for (int p = 4; p <= 20; ++p) { |
| 280 | + print_result(run_one(static_cast<std::size_t>(1) << p, forced_iters)); |
| 281 | + } |
| 282 | + return 0; |
| 283 | + } catch (const std::exception& e) { |
| 284 | + std::fprintf(stderr, "dif_vs_dip_benchmark failed: %s\n", e.what()); |
| 285 | + return 1; |
| 286 | + } |
| 287 | +} |
0 commit comments