NDFT blocked recurrence

[jenskeiner]

This PR seeks to improve the efficiency as well as the accuracy of the 1D direct NDFT forward/adjoint transforms for serial and threaded cases.

I was originally just trying to improve the performance but then noticed that the forward transform's error (max abs node-wise error divided by 1-norm of Fourier coefficients vector) grows like O(sqrt(N)) and that of the adjoint transform like O(N) while the error bound we use in our tests is just O(eps), independent of N.

After some rounding error analysis, the two sources of errors were

1. the calculation of the argument to sin/cos, and

2. the final summation over the frequencies.

3. can be reduced by using an FMA to reduce the argument to [-pi,pi] before feeding into sin/cos. This change has removed the N-dependency of the errors entirely.

After confirming the accuracy improvements, I applied another optimization that avoids the repeated costly sin/cos evaluations, by pulling out a constant factor added to the phase from one loop iteration to the next. This has to be done block-wise to not remove the accuracy improvements from the first set of changes again. There's now a tunable parameter B=32, the block size, that trades accuracy (B smaller) for speed (B larger).

I have attached a HTML write-up generated by my coding assistant that describes the changes and the rationale:
direct-1d-recurrence-optimization.html

This was tested on ARM64 (macOS) as well as in CI on x64 (Linux) architectures. If a platform does not support the single-rounding FMA semantics, the error bound could still be N-dependent again.

[codspeed-hq]

Merging this PR will improve performance by ×3.8

⚠️ Different runtime environments detected

Some benchmarks with significant performance changes were compared across different runtime environments,
which may affect the accuracy of the results.

Open the report in CodSpeed to investigate

⚡ 56 improved benchmarks
✅ 76 untouched benchmarks

Performance Changes

	Mode	Benchmark	BASE	HEAD	Efficiency
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_adjoint_direct_1d[512/1600]	41.5 ms	5.1 ms	×8.2
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_adjoint_direct_1d[256/800]	10.5 ms	1.3 ms	×8.2
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_adjoint_direct_1d[64/200]	660.1 µs	81.8 µs	×8.1
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_adjoint_direct_1d[128/400]	2,641.1 µs	329.3 µs	×8
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_forward_direct_1d[512/1600]	41.3 ms	5.4 ms	×7.6
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_forward_direct_1d[256/800]	10.4 ms	1.4 ms	×7.6
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_forward_direct_1d[64/200]	656.5 µs	87 µs	×7.6
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_forward_direct_1d_omp[512/1600]	2,998.4 µs	398 µs	×7.5
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_forward_direct_1d[128/400]	2,634.1 µs	352 µs	×7.5
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_forward_direct_1d_omp[256/800]	778.1 µs	110.1 µs	×7.1
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_adjoint_direct_1d[32/100]	160.3 µs	23.6 µs	×6.8
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_forward_direct_1d[32/100]	160.5 µs	25.3 µs	×6.3
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_forward_direct_1d_omp[128/400]	207.5 µs	34.5 µs	×6
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_adjoint_direct_1d_omp[512/1600]	3,299.8 µs	609.4 µs	×5.4
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_float/nfft_adjoint_direct_1d[512/1600]	21.6 ms	4.4 ms	×5
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_adjoint_direct_1d_omp[256/800]	859.3 µs	184 µs	×4.7
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_float/nfft_adjoint_direct_1d[256/800]	5 ms	1.1 ms	×4.5
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_float/nfft_forward_direct_1d_omp[512/1600]	1,556.6 µs	353.1 µs	×4.4
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_float/nfft_forward_direct_1d[512/1600]	21.2 ms	4.8 ms	×4.4
⚡	WallTime	codspeed-macro_gcc_kaiserbessel_double/nfft_forward_direct_1d_omp[64/200]	57.8 µs	14.1 µs	×4.1
...	...	...	...	...	...

ℹ️ Only the first 20 benchmarks are displayed. Go to the app to view all benchmarks.

Tip

Curious why this is faster? Comment @codspeedbot explain why this is faster on this PR, or directly use the CodSpeed MCP with your agent.

---

_{Comparing feature/nfft-direct-blocked-recurrence (83c6882) with develop (1e1cd06)}

[github-actions]

Accuracy report

2219 unchanged · 141 improved · 0 regressed

module	improved	regressed
nfft	141	0

📊 Full accuracy report

[jenskeiner]

By the way, there's some noise on the CodSpeed report regarding changes to the runtime environment. But upon closer inspection, the warning alerts on a change to the NFFT library that gets linked to the benchmark binaries. This will of course be different as changes are made to the code. Looking for a way to remove this noise...