Skip to content

Commit 2340ed8

Browse files
falseywinchnetfalseywinchnet
committed
finished DIP
1 parent 63d6d52 commit 2340ed8

6 files changed

Lines changed: 668 additions & 79 deletions

File tree

‎Makefile‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -84,6 +84,7 @@ BENCH := $(BUILD_DIR)/examples/benchmark
8484
BODFT_BENCH := $(BUILD_DIR)/examples/bodft_benchmark
8585
DIP_BENCH := $(BUILD_DIR)/examples/dip_benchmark
8686
DIF_DIP_BENCH := $(BUILD_DIR)/examples/dif_vs_dip_benchmark
87+
DIF_DIT_DIP_BENCH := $(BUILD_DIR)/examples/dif_dit_dip_benchmark
8788
DIF_DIT_BENCH := $(BUILD_DIR)/examples/dif_vs_dit_benchmark
8889
DIF_DIT_F32_BENCH := $(BUILD_DIR)/examples/dif_vs_dit_f32_benchmark
8990
APPLE_BENCH := $(BUILD_DIR)/examples/apple_benchmark
@@ -138,9 +139,9 @@ $(PC_FILE): pkgconfig/bfft.pc.in | $(BUILD_DIR)
138139
-e 's|@PROJECT_VERSION@|$(VERSION)|g' \
139140
$< > $@
140141

141-
examples: $(BENCH) $(BODFT_BENCH) $(DIP_BENCH) $(DIF_DIP_BENCH) $(APPLE_EXAMPLES) $(LOCALITY_PROBE) $(C_DEMO) $(CPP_DEMO)
142+
examples: $(BENCH) $(BODFT_BENCH) $(DIP_BENCH) $(DIF_DIP_BENCH) $(DIF_DIT_DIP_BENCH) $(APPLE_EXAMPLES) $(LOCALITY_PROBE) $(C_DEMO) $(CPP_DEMO)
142143

143-
benchmarks: $(BENCH) $(BODFT_BENCH) $(DIP_BENCH) $(DIF_DIP_BENCH) $(DIF_DIT_BENCH) $(DIF_DIT_F32_BENCH) $(ATAN2_BENCH)
144+
benchmarks: $(BENCH) $(BODFT_BENCH) $(DIP_BENCH) $(DIF_DIP_BENCH) $(DIF_DIT_DIP_BENCH) $(DIF_DIT_BENCH) $(DIF_DIT_F32_BENCH) $(ATAN2_BENCH)
144145

145146
asm-check: $(ASM_OUTPUTS)
146147
@if [ -z "$(ASM_OUTPUTS)" ]; then echo "No x86 assembly variants supported by $(CXX)."; fi
@@ -166,6 +167,9 @@ $(DIP_BENCH): examples/dip_benchmark.cpp src/detail/bruun_dip_kernel.hpp src/det
166167
$(DIF_DIP_BENCH): examples/dif_vs_dip_benchmark.cpp src/detail/bruun_dip_kernel.hpp src/detail/bruun_dif_kernel.hpp src/detail/bruun_simd_backend.hpp $(STATIC_LIB) | $(BUILD_DIR)
167168
$(CXX) $(CPPFLAGS) $(INCLUDES) $(CXXFLAGS) $(AUTO_SIMD_FLAGS) $< $(STATIC_LIB) $(LDLIBS) -o $@
168169

170+
$(DIF_DIT_DIP_BENCH): examples/dif_dit_dip_benchmark.cpp src/detail/bruun_dip_kernel.hpp src/detail/bruun_dit_kernel.hpp src/detail/bruun_dif_kernel.hpp src/detail/bruun_simd_backend.hpp src/detail/MAG_REPRESENT_KERNEL.hpp include/bfft/bfft.hpp $(STATIC_LIB) | $(BUILD_DIR)
171+
$(CXX) $(CPPFLAGS) $(INCLUDES) $(CXXFLAGS) $(AUTO_SIMD_FLAGS) $< $(STATIC_LIB) $(LDLIBS) -o $@
172+
169173
$(DIF_DIT_BENCH): examples/dif_vs_dit_benchmark.cpp src/detail/bruun_dit_kernel.hpp src/detail/bruun_dif_kernel.hpp src/detail/MAG_REPRESENT_KERNEL.hpp include/bfft/bfft.hpp $(STATIC_LIB) | $(BUILD_DIR)
170174
$(CXX) $(CPPFLAGS) $(INCLUDES) $(CXXFLAGS) $(AUTO_SIMD_FLAGS) $< $(STATIC_LIB) $(LDLIBS) -o $@
171175

‎examples/dif_dit_dip_benchmark.cpp‎

Lines changed: 333 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,333 @@
1+
// Three-way real-FFT microbenchmark: DIF vs DIT vs DIP, forward / inverse /
2+
// roundtrip, f64. DIF is the ratio baseline (its residue kernel is the tuned
3+
// reference the other two are chasing).
4+
//
5+
// Hardening (carried over from dif_vs_dip_benchmark, extended to 3 engines):
6+
// - INTERLEAVED passes: every pass times a chunk of DIF, then DIT, then DIP,
7+
// so all three ride the same thermal / frequency-scaling trajectory and
8+
// the "measure one fully then the next" ordering bias is cancelled.
9+
// - MEDIAN of many passes (not best-of): robust to transient spikes.
10+
// - CALIBRATED iteration counts: each timed chunk runs >= kTargetChunkNs, so
11+
// sub-microsecond transforms at small N get thousands of iterations and
12+
// become measurable instead of quantized to the clock.
13+
14+
#include "../src/detail/bruun_dif_kernel.hpp"
15+
#include "../src/detail/bruun_dit_kernel.hpp"
16+
#include "../src/detail/bruun_dip_kernel.hpp"
17+
18+
#include <algorithm>
19+
#include <chrono>
20+
#include <cmath>
21+
#include <cstdio>
22+
#include <cstdlib>
23+
#include <exception>
24+
#include <random>
25+
#include <stdexcept>
26+
#include <vector>
27+
28+
namespace {
29+
30+
using clock_type = std::chrono::steady_clock;
31+
32+
constexpr double pi = 3.141592653589793238462643383279502884;
33+
constexpr double kTargetChunkNs = 2.0e6;
34+
constexpr int kPasses = 11;
35+
36+
struct triple_timing {
37+
double a_ns = 0.0; // DIF
38+
double b_ns = 0.0; // DIT
39+
double c_ns = 0.0; // DIP
40+
double sink = 0.0;
41+
};
42+
43+
struct result {
44+
std::size_t n = 0;
45+
int iters = 0;
46+
double fwd[3] = {0, 0, 0}; // DIF, DIT, DIP
47+
double inv[3] = {0, 0, 0};
48+
double rt[3] = {0, 0, 0};
49+
double dit_fwd_err = 0.0;
50+
double dip_fwd_err = 0.0;
51+
double dit_rt_err = 0.0;
52+
double dip_rt_err = 0.0;
53+
double sink = 0.0;
54+
};
55+
56+
bool is_power2(std::size_t n) {
57+
return n > 0 && (n & (n - 1)) == 0;
58+
}
59+
60+
std::size_t parse_size(const char* text) {
61+
char* end = nullptr;
62+
const unsigned long long value = std::strtoull(text, &end, 0);
63+
if (!end || *end != '\0' || value == 0) {
64+
throw std::invalid_argument("invalid size");
65+
}
66+
return static_cast<std::size_t>(value);
67+
}
68+
69+
int parse_iters(const char* text) {
70+
char* end = nullptr;
71+
const long value = std::strtol(text, &end, 0);
72+
if (!end || *end != '\0' || value <= 0) {
73+
throw std::invalid_argument("invalid iteration count");
74+
}
75+
return static_cast<int>(value);
76+
}
77+
78+
template <typename Func>
79+
int calibrate_iters(Func&& func) {
80+
const auto start = clock_type::now();
81+
double sink = 0.0;
82+
for (int i = 0; i < 4; ++i) {
83+
sink += func(i, sink);
84+
}
85+
const auto stop = clock_type::now();
86+
const double per_op =
87+
std::chrono::duration<double, std::nano>(stop - start).count() / 4.0 + 1.0;
88+
double iters = kTargetChunkNs / per_op;
89+
if (iters < 3.0) iters = 3.0;
90+
if (iters > 262144.0) iters = 262144.0;
91+
return static_cast<int>(iters) + (sink == 0.12345 ? 1 : 0);
92+
}
93+
94+
double median_of(double* v, int count) {
95+
std::sort(v, v + count);
96+
return (count & 1) ? v[count / 2] : 0.5 * (v[count / 2 - 1] + v[count / 2]);
97+
}
98+
99+
// Interleaved three-engine measurement: every pass times a chunk of A, then B,
100+
// then C, so all ride the same thermal trajectory. Reported value is the
101+
// per-engine median of the per-pass times.
102+
template <typename FA, typename FB, typename FC>
103+
triple_timing bench_triple(int iters, FA&& fa, FB&& fb, FC&& fc) {
104+
double as[kPasses], bs[kPasses], cs[kPasses];
105+
triple_timing out;
106+
for (int pass = 0; pass < kPasses; ++pass) {
107+
double sink = 0.0;
108+
auto t0 = clock_type::now();
109+
for (int i = 0; i < iters; ++i) sink += fa(i, sink);
110+
auto t1 = clock_type::now();
111+
for (int i = 0; i < iters; ++i) sink += fb(i, sink);
112+
auto t2 = clock_type::now();
113+
for (int i = 0; i < iters; ++i) sink += fc(i, sink);
114+
auto t3 = clock_type::now();
115+
as[pass] = std::chrono::duration<double, std::nano>(t1 - t0).count() / iters;
116+
bs[pass] = std::chrono::duration<double, std::nano>(t2 - t1).count() / iters;
117+
cs[pass] = std::chrono::duration<double, std::nano>(t3 - t2).count() / iters;
118+
out.sink += sink;
119+
}
120+
out.a_ns = median_of(as, kPasses);
121+
out.b_ns = median_of(bs, kPasses);
122+
out.c_ns = median_of(cs, kPasses);
123+
return out;
124+
}
125+
126+
template <typename A, typename B>
127+
double max_abs_complex(const std::vector<A>& a, const std::vector<B>& b) {
128+
double err = 0.0;
129+
for (std::size_t i = 0; i < a.size(); ++i) {
130+
err = std::max(err, std::abs(a[i].re - b[i].re));
131+
err = std::max(err, std::abs(a[i].im - b[i].im));
132+
}
133+
return err;
134+
}
135+
136+
double max_abs_real(const std::vector<double>& a, const std::vector<double>& b) {
137+
double err = 0.0;
138+
for (std::size_t i = 0; i < a.size(); ++i) {
139+
err = std::max(err, std::abs(a[i] - b[i]));
140+
}
141+
return err;
142+
}
143+
144+
std::vector<double> make_signal(std::size_t n) {
145+
std::vector<double> input(n);
146+
std::mt19937_64 rng(0xD1F0D17DULL + static_cast<unsigned long long>(n));
147+
std::uniform_real_distribution<double> noise(-0.05, 0.05);
148+
for (std::size_t i = 0; i < n; ++i) {
149+
const double t = static_cast<double>(i) / static_cast<double>(n);
150+
input[i] = std::sin(2.0 * pi * 13.0 * t)
151+
+ 0.5 * std::cos(2.0 * pi * 37.0 * t)
152+
+ 0.25 * std::sin(2.0 * pi * 89.0 * t)
153+
+ noise(rng);
154+
}
155+
return input;
156+
}
157+
158+
result run_one(std::size_t n, int forced_iters) {
159+
if (n > static_cast<std::size_t>(2147483647)) {
160+
throw std::invalid_argument("N is too large for the experimental kernels");
161+
}
162+
163+
bruun::DIF_RFFT_kernel dif;
164+
bruun::DIT_RFFT_kernel dit;
165+
bruun::DIP_RFFT_kernel dip;
166+
if (!dif.init(static_cast<int>(n)) || !dit.init(static_cast<int>(n)) ||
167+
!dip.init(static_cast<int>(n))) {
168+
throw std::runtime_error("kernel setup failed");
169+
}
170+
171+
const std::size_t nb = n / 2 + 1;
172+
const std::vector<double> original = make_signal(n);
173+
std::vector<double> input(original);
174+
175+
std::vector<bruun::complex_t> dif_bins(nb), dit_bins(nb), dip_bins(nb);
176+
std::vector<bruun::complex_t> dif_scratch(static_cast<std::size_t>(dif.native_scratch_size()));
177+
std::vector<double> dif_work(static_cast<std::size_t>(dif.work_size()));
178+
std::vector<double> dit_work(static_cast<std::size_t>(dit.work_size()));
179+
std::vector<double> dip_work(static_cast<std::size_t>(dip.work_size()));
180+
std::vector<double> dif_out(n), dit_out(n), dip_out(n);
181+
182+
// warm + correctness reference
183+
dif.forward_standard(original.data(), dif_bins.data(), dif_work.data(), dif_scratch.data());
184+
dit.forward_simd(original.data(), dit_bins.data(), dit_work.data());
185+
dip.forward_standard(original.data(), dip_bins.data(), dip_work.data());
186+
dif.inverse(dif_bins.data(), dif_out.data());
187+
dit.inverse_simd(dit_bins.data(), dit_out.data(), dit_work.data());
188+
dip.inverse_standard(dip_bins.data(), dip_out.data(), dip_work.data());
189+
190+
result r;
191+
r.n = n;
192+
r.dit_fwd_err = max_abs_complex(dif_bins, dit_bins);
193+
r.dip_fwd_err = max_abs_complex(dif_bins, dip_bins);
194+
r.dit_rt_err = max_abs_real(original, dit_out);
195+
r.dip_rt_err = max_abs_real(original, dip_out);
196+
197+
auto dif_fwd = [&](int i, double sink) {
198+
input[(static_cast<std::size_t>(i) * 131u + static_cast<std::size_t>(sink)) & (n - 1)] += 1e-12;
199+
dif.forward_standard(input.data(), dif_bins.data(), dif_work.data(), dif_scratch.data());
200+
return dif_bins[(static_cast<std::size_t>(i) * 17u) % nb].re;
201+
};
202+
auto dit_fwd = [&](int i, double sink) {
203+
input[(static_cast<std::size_t>(i) * 131u + static_cast<std::size_t>(sink)) & (n - 1)] += 1e-12;
204+
dit.forward_simd(input.data(), dit_bins.data(), dit_work.data());
205+
return dit_bins[(static_cast<std::size_t>(i) * 17u) % nb].re;
206+
};
207+
auto dip_fwd = [&](int i, double sink) {
208+
input[(static_cast<std::size_t>(i) * 131u + static_cast<std::size_t>(sink)) & (n - 1)] += 1e-12;
209+
dip.forward_standard(input.data(), dip_bins.data(), dip_work.data());
210+
return dip_bins[(static_cast<std::size_t>(i) * 17u) % nb].re;
211+
};
212+
auto dif_inv = [&](int i, double) {
213+
dif_bins[(static_cast<std::size_t>(i) * 17u) % nb].re += 1e-12;
214+
dif.inverse(dif_bins.data(), dif_out.data());
215+
return dif_out[(static_cast<std::size_t>(i) * 31u) & (n - 1)];
216+
};
217+
auto dit_inv = [&](int i, double) {
218+
dit_bins[(static_cast<std::size_t>(i) * 17u) % nb].re += 1e-12;
219+
dit.inverse_simd(dit_bins.data(), dit_out.data(), dit_work.data());
220+
return dit_out[(static_cast<std::size_t>(i) * 31u) & (n - 1)];
221+
};
222+
auto dip_inv = [&](int i, double) {
223+
dip_bins[(static_cast<std::size_t>(i) * 17u) % nb].re += 1e-12;
224+
dip.inverse_standard(dip_bins.data(), dip_out.data(), dip_work.data());
225+
return dip_out[(static_cast<std::size_t>(i) * 31u) & (n - 1)];
226+
};
227+
auto dif_rt = [&](int i, double sink) {
228+
input[(static_cast<std::size_t>(i) * 131u + static_cast<std::size_t>(sink)) & (n - 1)] += 1e-12;
229+
dif.forward_standard(input.data(), dif_bins.data(), dif_work.data(), dif_scratch.data());
230+
dif.inverse(dif_bins.data(), dif_out.data());
231+
return dif_out[(static_cast<std::size_t>(i) * 31u) & (n - 1)];
232+
};
233+
auto dit_rt = [&](int i, double sink) {
234+
input[(static_cast<std::size_t>(i) * 131u + static_cast<std::size_t>(sink)) & (n - 1)] += 1e-12;
235+
dit.forward_simd(input.data(), dit_bins.data(), dit_work.data());
236+
dit.inverse_simd(dit_bins.data(), dit_out.data(), dit_work.data());
237+
return dit_out[(static_cast<std::size_t>(i) * 31u) & (n - 1)];
238+
};
239+
auto dip_rt = [&](int i, double sink) {
240+
input[(static_cast<std::size_t>(i) * 131u + static_cast<std::size_t>(sink)) & (n - 1)] += 1e-12;
241+
dip.forward_standard(input.data(), dip_bins.data(), dip_work.data());
242+
dip.inverse_standard(dip_bins.data(), dip_out.data(), dip_work.data());
243+
return dip_out[(static_cast<std::size_t>(i) * 31u) & (n - 1)];
244+
};
245+
246+
input = original;
247+
const int iters = forced_iters > 0 ? forced_iters : calibrate_iters(dif_fwd);
248+
r.iters = iters;
249+
250+
input = original;
251+
triple_timing t = bench_triple(iters, dif_fwd, dit_fwd, dip_fwd);
252+
r.fwd[0] = t.a_ns; r.fwd[1] = t.b_ns; r.fwd[2] = t.c_ns;
253+
r.sink += t.sink;
254+
255+
dif.forward_standard(original.data(), dif_bins.data(), dif_work.data(), dif_scratch.data());
256+
dit.forward_simd(original.data(), dit_bins.data(), dit_work.data());
257+
dip.forward_standard(original.data(), dip_bins.data(), dip_work.data());
258+
t = bench_triple(iters, dif_inv, dit_inv, dip_inv);
259+
r.inv[0] = t.a_ns; r.inv[1] = t.b_ns; r.inv[2] = t.c_ns;
260+
r.sink += t.sink;
261+
262+
input = original;
263+
t = bench_triple(iters, dif_rt, dit_rt, dip_rt);
264+
r.rt[0] = t.a_ns; r.rt[1] = t.b_ns; r.rt[2] = t.c_ns;
265+
r.sink += t.sink;
266+
267+
return r;
268+
}
269+
270+
void print_header() {
271+
std::printf("%9s %8s %11s %11s %11s %11s %11s %11s %11s %11s %11s %8s %8s %8s %8s %8s %8s %s\n",
272+
"N", "iters",
273+
"DIF_fwd", "DIT_fwd", "DIP_fwd",
274+
"DIF_inv", "DIT_inv", "DIP_inv",
275+
"DIF_rt", "DIT_rt", "DIP_rt",
276+
"ditf_x", "dipf_x", "diti_x", "dipi_x", "ditr_x", "dipr_x",
277+
"checks");
278+
std::fflush(stdout);
279+
}
280+
281+
void print_result(const result& r) {
282+
std::printf("%9zu %8d %11.2f %11.2f %11.2f %11.2f %11.2f %11.2f %11.2f %11.2f %11.2f "
283+
"%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f "
284+
"dit_fe %.1e dip_fe %.1e dit_re %.1e dip_re %.1e\n",
285+
r.n, r.iters,
286+
r.fwd[0], r.fwd[1], r.fwd[2],
287+
r.inv[0], r.inv[1], r.inv[2],
288+
r.rt[0], r.rt[1], r.rt[2],
289+
r.fwd[1] / r.fwd[0], r.fwd[2] / r.fwd[0],
290+
r.inv[1] / r.inv[0], r.inv[2] / r.inv[0],
291+
r.rt[1] / r.rt[0], r.rt[2] / r.rt[0],
292+
r.dit_fwd_err, r.dip_fwd_err, r.dit_rt_err, r.dip_rt_err);
293+
std::fflush(stdout);
294+
}
295+
296+
} // namespace
297+
298+
int main(int argc, char** argv) {
299+
try {
300+
if (argc > 3) {
301+
std::fprintf(stderr, "usage: %s [N | max_pow] [iters]\n", argv[0]);
302+
return 2;
303+
}
304+
305+
int forced_iters = 0;
306+
if (argc >= 3) {
307+
forced_iters = parse_iters(argv[2]);
308+
}
309+
310+
print_header();
311+
if (argc >= 2) {
312+
const std::size_t arg = parse_size(argv[1]);
313+
if (is_power2(arg) && arg >= 4) {
314+
print_result(run_one(arg, forced_iters));
315+
} else if (arg >= 2 && arg <= 30) {
316+
for (int p = 2; p <= static_cast<int>(arg); ++p) {
317+
print_result(run_one(static_cast<std::size_t>(1) << p, forced_iters));
318+
}
319+
} else {
320+
throw std::invalid_argument("argument must be a power-of-two N or max_pow in [2, 30]");
321+
}
322+
return 0;
323+
}
324+
325+
for (int p = 4; p <= 20; ++p) {
326+
print_result(run_one(static_cast<std::size_t>(1) << p, forced_iters));
327+
}
328+
return 0;
329+
} catch (const std::exception& e) {
330+
std::fprintf(stderr, "dif_dit_dip_benchmark failed: %s\n", e.what());
331+
return 1;
332+
}
333+
}

0 commit comments

Comments
 (0)