add FFT size-sweep benchmark vs FFTW with gnuplot output

- benchmark/algorithms_fft.cpp: sweeps 1-D (16..2^21), 2-D (16^2..2048^2),
  3-D (8^3..256^3); plan-recycled execution timing, cache flushed per rep,
  FFTW_ESTIMATE, single thread; writes gnuplot-friendly fft_bench_{1d,2d,3d}.dat
  (n, N, ms, benchFFT MFLOPS, mine/FFTW ratio)
- benchmark/algorithms_fft_plots.gp: renders the three speed-vs-size plots
- register algorithms_fft.cpp.x in the standalone benchmark project
- fft.NOTES.md: add CUDA adaptation plan (section 8)

Measured on this machine: faster than FFTW for 1-D n <= 128, 2-D >= 256^2
(1.4-1.7x), and 3-D 256^3 (2.2x); FFTW leads mid-range 1-D (SIMD codelets).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: alfC

benchmark/CMakeLists.txt

@@ -52,6 +52,11 @@ target_link_libraries(matrix_multiplication.cpp.x benchmark::benchmark)
 target_link_libraries(matrix_multiplication.cpp.x -ltbb)
 add_test(NAME matrix_multiplication.cpp.x COMMAND matrix_multiplication.cpp.x)
 
+add_executable(algorithms_fft.cpp.x algorithms_fft.cpp)
+target_link_libraries(algorithms_fft.cpp.x fftw3)
+target_include_directories(algorithms_fft.cpp.x PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../include)
+add_test(NAME algorithms_fft.cpp.x COMMAND algorithms_fft.cpp.x)
+
 add_executable(algorithms_gemm.cpp.x algorithms_gemm.cpp)
 target_link_libraries(algorithms_gemm.cpp.x benchmark::benchmark Boost::unit_test_framework)
 target_link_libraries(algorithms_gemm.cpp.x -ltbb)

benchmark/algorithms_fft.cpp

@@ -0,0 +1,127 @@
+// Copyright 2024-2026 Alfredo A. Correa
+// Distributed under the Boost Software License, Version 1.0.
+// https://www.boost.org/LICENSE_1_0.txt
+
+#ifdef COMPILATION_INSTRUCTIONS
+g++ - std = c++ 17 - O3 - march = native - DNDEBUG - I../ include $0 - o $0.x - lfftw3 &&./ $0.x && gnuplot algorithms_fft_plots.gp;
+exit
+#endif
+
+// Size sweep of multi::fft_plan vs FFTW, emitting gnuplot-friendly .dat files
+// (one per dimensionality). Methodology: both libraries build their plan once
+// (untimed) and recycle it; only execution is timed; CPU cache flushed before
+// every repetition; FFTW_ESTIMATE, no wisdom; single thread.
+#include <boost/multi/algorithms/fft.hpp>
+#include <boost/multi/array.hpp>
+
+#include <fftw3.h>
+
+#include <algorithm>
+#include <chrono>
+#include <cmath>
+#include <complex>
+#include <cstdio>
+#include <random>
+#include <vector>
+
+	namespace multi = boost::multi;
+using complex       = std::complex<double>;
+
+struct watch {
+	std::chrono::high_resolution_clock::time_point s = std::chrono::high_resolution_clock::now();
+	auto                                           sec() const { return std::chrono::duration<double>(std::chrono::high_resolution_clock::now() - s).count(); }
+};
+static std::vector<char> g_thrash(64 << 20);
+void                     flush_cache() {
+    for(std::size_t i = 0; i < g_thrash.size(); i += 64) {
+        g_thrash[i]++;
+    }
+    char volatile x = g_thrash[0];
+    (void)x;
+}
+template<class F> double time_it(long reps, F f) {
+	double t = 0;
+	for(long r = 0; r != reps; ++r) {
+		flush_cache();
+		watch w;
+		f();
+		t += w.sec();
+	}
+	return t / reps;
+}
+
+auto reps_for(double n_total) -> long {
+	double const work = 5.0 * n_total * std::log2(std::max(n_total, 2.0));
+	return std::clamp<long>(static_cast<long>(2e8 / work), 5, 300);
+}
+
+template<std::ptrdiff_t D>
+void sweep(std::vector<int> const& sides, char const* fname, char const* label) {
+	std::FILE* out = std::fopen(fname, "w");
+	std::fprintf(out, "# %s: multi::fft_plan vs FFTW 3 (plan recycled, exec only, cache flushed/rep, single thread, FFTW_ESTIMATE)\n", label);
+	std::fprintf(out, "# mflops = 5*N*log2(N)/time_us (benchFFT convention), N = total points\n");
+	std::fprintf(out, "# n_side  N_total  mine_ms  fftw_ms  mine_mflops  fftw_mflops  ratio_mine_over_fftw\n");
+	for(int n : sides) {
+		long N = 1;
+		for(int d = 0; d != D; ++d) {
+			N *= n;
+		}
+		std::vector<complex>                   base(N);
+		std::mt19937                           gen(42);
+		std::uniform_real_distribution<double> dist(-1.0, 1.0);
+		for(auto& e : base) {
+			e = complex{dist(gen), dist(gen)};
+		}
+		long const reps = reps_for(static_cast<double>(N));
+
+		multi::array<complex, 1> flat(multi::extensions_t<1>{N});
+		auto                     load = [&] { std::copy(base.begin(), base.end(), flat.begin()); };
+		double                   mine = 0.0;
+		{
+			std::array<multi::ssize_t, D> ext_arr{};
+			ext_arr.fill(n);
+			multi::fft_plan<complex, D> const plan{ext_arr, multi::fft_forward};
+			if constexpr(D == 1) {
+				load();
+				plan(flat);
+				mine = time_it(reps, [&] { load(); plan(flat); });
+			} else if constexpr(D == 2) {
+				multi::array_ref<complex, 2> v(flat.data_elements(), {n, n});
+				load();
+				plan(v);
+				mine = time_it(reps, [&] { load(); plan(v); });
+			} else {
+				multi::array_ref<complex, 3> v(flat.data_elements(), {n, n, n});
+				load();
+				plan(v);
+				mine = time_it(reps, [&] { load(); plan(v); });
+			}
+		}
+
+		auto* in    = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * N));
+		auto* fo    = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * N));
+		auto  loadf = [&] { for(long i = 0; i != N; ++i) { in[i][0] = base[i].real(); in[i][1] = base[i].imag(); } };
+		fftw_forget_wisdom();
+		fftw_plan p = D == 1 ? fftw_plan_dft_1d(n, in, fo, FFTW_FORWARD, FFTW_ESTIMATE)
+					: D == 2 ? fftw_plan_dft_2d(n, n, in, fo, FFTW_FORWARD, FFTW_ESTIMATE)
+							 : fftw_plan_dft_3d(n, n, n, in, fo, FFTW_FORWARD, FFTW_ESTIMATE);
+		loadf();
+		double ffw = time_it(reps, [&] { loadf(); fftw_execute(p); });
+		fftw_destroy_plan(p);
+		fftw_free(in);
+		fftw_free(fo);
+
+		double const work = 5.0 * static_cast<double>(N) * std::log2(static_cast<double>(N));
+		std::fprintf(out, "%8d %10ld %12.5f %12.5f %12.1f %12.1f %8.3f\n", n, N, mine * 1e3, ffw * 1e3, work / (mine * 1e6), work / (ffw * 1e6), mine / ffw);
+		std::fflush(out);
+		std::fprintf(stderr, "%s n=%d done (reps=%ld)\n", label, n, reps);
+	}
+	std::fclose(out);
+}
+
+int main() {
+	sweep<1>({16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768, 65536, 131072, 262144, 524288, 1048576, 2097152}, "fft_bench_1d.dat", "1D n");
+	sweep<2>({16, 32, 64, 128, 256, 512, 1024, 2048}, "fft_bench_2d.dat", "2D n x n");
+	sweep<3>({8, 16, 32, 64, 128, 256}, "fft_bench_3d.dat", "3D n x n x n");
+	return 0;
+}

benchmark/algorithms_fft_plots.gp

@@ -0,0 +1,32 @@
+# Plots multi::fft_plan vs FFTW 3 from fft_bench_{1d,2d,3d}.dat
+# usage: run algorithms_fft.cpp.x first (writes the .dat files to the
+# current directory), then: gnuplot algorithms_fft_plots.gp
+# data columns: n_side N_total mine_ms fftw_ms mine_mflops fftw_mflops ratio
+
+set terminal pngcairo size 900,600 font ",11"
+set grid xtics ytics mytics
+set logscale x 2
+set format x "2^{%L}"                 # x tics as powers of two
+set xlabel "transform size n (per axis)"
+set ylabel "speed  5·N·log_2(N) / time   [MFLOPS]  (higher is better)"
+set key top left
+set style line 1 lw 2 pt 7  ps 1.1 lc rgb "#c0392b"   # multi
+set style line 2 lw 2 pt 5  ps 1.1 lc rgb "#2c3e50"   # fftw
+
+set format x "%.0f"
+set xtics rotate by -30
+
+set output "fft_bench_1d.png"
+set title "1-D complex FFT — multi::fft\\_plan vs FFTW 3 (plan recycled, single thread)"
+plot "fft_bench_1d.dat" using 1:5 with linespoints ls 1 title "multi::fft\\_plan", \
+     "fft_bench_1d.dat" using 1:6 with linespoints ls 2 title "FFTW 3 (estimate)"
+
+set output "fft_bench_2d.png"
+set title "2-D complex FFT (n × n) — multi::fft\\_plan vs FFTW 3"
+plot "fft_bench_2d.dat" using 1:5 with linespoints ls 1 title "multi::fft\\_plan", \
+     "fft_bench_2d.dat" using 1:6 with linespoints ls 2 title "FFTW 3 (estimate)"
+
+set output "fft_bench_3d.png"
+set title "3-D complex FFT (n × n × n) — multi::fft\\_plan vs FFTW 3"
+plot "fft_bench_3d.dat" using 1:5 with linespoints ls 1 title "multi::fft\\_plan", \
+     "fft_bench_3d.dat" using 1:6 with linespoints ls 2 title "FFTW 3 (estimate)"

include/boost/multi/algorithms/fft.NOTES.md

@@ -396,3 +396,41 @@ expressing butterflies through strided iterators would block the tight,
 auto-vectorizable batched loops that make this competitive with FFTW. The
 boundary between the two is the gather/scatter (or, on the fast paths, a raw
 `base()` pointer handed across) — FFTW draws the same line.
+
+---
+
+## 8. CUDA adaptation plan (not yet implemented)
+
+The architecture is already CUDA-shaped: Stockham (no bit-reversal, uniform
+stages) is what GPU FFT libraries use; the batched layout `buf[k*m + j]`
+(batch contiguous) is exactly warp-coalesced when `j` maps to `threadIdx.x`;
+the fused strided stages' (pointer, stride, width) interface is
+layout-agnostic; and plans already separate "build tables once" (host) from
+"execute many" (upload once, launch many). Incremental steps:
+
+1. **Extract butterfly bodies** into `BOOST_MULTI_HD` inline functions (the
+   per-(block, r, j) work) shared by the host loops and device kernels. Pure
+   refactor; CPU codegen must not regress (re-benchmark).
+2. **`fft_memory_space<Ptr>` trait** dispatching host (raw `T*`, current
+   paths) vs device (thrust/Multi-CUDA fancy pointers), specialized in an
+   opt-in adaptor header so the core stays CUDA-free. The current
+   host-iterator gather fallback must *never* run on device pointers.
+3. **Device table mirror**: keep host `std::vector` tables; add a lazily
+   uploaded POD `engine_view` (raw device pointers + sizes) per engine.
+4. **Device stage launchers**: one `__global__` kernel per stage kind around
+   the shared butterfly body; thread → (j, r, block); stream-ordered launches
+   replace the host stage loop. Grids are natively 3-D, so an extra
+   (outer-slab, stride) index runs a whole 3-D middle axis in one launch.
+5. **Granularity flip**: same rotation/rank-descent orchestration, but the
+   device branch batches the whole axis (m = all fibers; `mb_` is a host
+   cache concern) and skips slab pairing (a CPU cache optimization). Needs a
+   cap-and-tile fallback for the 2·n·m device scratch on huge arrays.
+6. **Bluestein**: pointwise chirp/kernel multiplies as elementwise kernels;
+   sub-engine recursion stays host-side (it only orchestrates launches).
+   **Six-step**: skip on device (its cache-blocking rationale is CPU-only).
+7. **Tests**: mirror `algorithms_fft.cpp` on Multi-thrust device arrays (the
+   test already compiles as CUDA in the nvcc lane); `thrust::complex` should
+   work unchanged through the `T{re,im}`/`+`/`-`/`*` contract and `fft_ops`.
+8. **Benchmark vs cuFFT** with the `adaptors/cufft` harness — the goal is
+   *generic* device FFT (custom element types, arbitrary strided layouts,
+   header-only), not beating cuFFT's tuned codelets for `float/double`.