From 229f50cd864d9768324042c10edbb56bfc120f08 Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Fri, 25 Sep 2026 23:35:45 -0400 Subject: [PATCH] Feat: Use FFTW linear convolution for CPU fast-mode PSF resolve Replace the nested spatial global convolution in fast_psf_accumulator_resolve with a reusable double-precision FFTW linear convolution on the CPU PSF backend. - Add the private src/fast_psf_fftw.{c,h} module: zero-padded R2C/C2R plans, cached kernel spectrum and planar scratch, exact 1/(Pwidth*Pheight) and 1/N^2 normalization, (R,R) crop, and additive HDR output. - Keep the previous nested loops as fast_psf_accumulator_resolve_spatial_reference for tests/benchmarks only; it is not a runtime fallback. - Cache the circular row spans on FastPsfAccumulator and report one-time plan, kernel transform, scratch, and per-frame stage timings. - Require fftw3_omp for CPU builds; HIP and dummy builds do not link FFTW. - Namespace test/helper binaries by spacetime and build tag, and reject make test / psf-capture for non-CPU backends. - Add tests/test_fast_psf_fftw.c (FFTW versus spatial), tests/benchmark_fast_psf_fftw.c, an FFTW CLI smoke check, and the 2026-09-25 benchmark record. --- Makefile | 125 ++-- benchmarks/fast_mode_fftw_2026-09-25.md | 280 +++++++++ build.md | 10 +- fast_mode_fftw_optimization_plan.md | 799 ++++++++++++++++++++++++ nr_spacetime_movie_renderer_design.md | 8 +- src/fast_psf_fftw.c | 393 ++++++++++++ src/fast_psf_fftw.h | 67 ++ src/main.c | 2 + src/optics.c | 95 ++- src/optics.h | 21 +- tests/benchmark_fast_psf_fftw.c | 223 +++++++ tests/test_camera_cli.py | 12 +- tests/test_fast_psf_fftw.c | 378 +++++++++++ usage.md | 14 + 14 files changed, 2372 insertions(+), 55 deletions(-) create mode 100644 benchmarks/fast_mode_fftw_2026-09-25.md create mode 100644 fast_mode_fftw_optimization_plan.md create mode 100644 src/fast_psf_fftw.c create mode 100644 src/fast_psf_fftw.h create mode 100644 tests/benchmark_fast_psf_fftw.c create mode 100644 tests/test_fast_psf_fftw.c diff --git a/Makefile b/Makefile index 514efd1..e83769b 100644 --- a/Makefile +++ b/Makefile @@ -29,22 +29,14 @@ else IMAGE_EXT := ppm endif -COMMON_SOURCES := $(filter-out src/main.c src/dummy_psf.c src/spacetime_minkowski.c src/spacetime_schwarzschild.c,$(wildcard src/*.c)) +COMMON_SOURCES := $(filter-out src/main.c src/dummy_psf.c src/fast_psf_fftw.c src/spacetime_minkowski.c src/spacetime_schwarzschild.c,$(wildcard src/*.c)) PROVIDER_SOURCE := src/spacetime_$(SPACETIME).c BUILD_DIR := build/$(BUILD_TYPE) TARGET_BASENAME := $(SPACETIME)_sky OBJECT_DIR := $(BUILD_DIR)/obj/$(SPACETIME) CORE_MINKOWSKI_SOURCES := $(COMMON_SOURCES) src/spacetime_minkowski.c -TEST_TARGET := $(BUILD_DIR)/test_geodesic -FRAME_TEST_TARGET := $(BUILD_DIR)/test_frame -SCHWARZSCHILD_TEST_TARGET := $(BUILD_DIR)/test_schwarzschild -OBSERVER_TRACK_TEST_TARGET := $(BUILD_DIR)/test_observer_track -CATALOG_PREFETCH_TEST_TARGET := $(BUILD_DIR)/test_catalog_prefetch -HIP_PSF_TEST_TARGET := $(BUILD_DIR)/test_hip_psf -HIP_PSF_BENCH_TARGET := $(BUILD_DIR)/benchmark_hip_psf -CAMERA_TEST_TARGETS := $(BUILD_DIR)/test_observer_minkowski $(BUILD_DIR)/test_observer_schwarzschild -.PHONY: all backend clean run test hip-psf-test hip-psf-bench minkowski schwarzschild FORCE +.PHONY: all backend clean run test hip-psf-test hip-psf-bench fast-psf-fftw-bench minkowski schwarzschild FORCE ifneq ($(filter 0 1,$(PSF_EVENT_SINK)),$(PSF_EVENT_SINK)) $(error Unknown PSF_EVENT_SINK '$(PSF_EVENT_SINK)'; choose 0 or 1) @@ -53,6 +45,15 @@ BUILD_CPPFLAGS += -DFRAME_PSF_EVENT_SINK=$(PSF_EVENT_SINK) ifeq ($(PSF_BACKEND),cpu) RENDER_LINKER := $(CC) +ifeq ($(shell pkg-config --exists fftw3_omp && echo yes),yes) +CPU_FFTW_CPPFLAGS := -DFAST_PSF_FFTW $(shell pkg-config --cflags fftw3_omp) +CPU_FFTW_LDLIBS := $(shell pkg-config --libs fftw3_omp) +CPU_FFTW_SOURCES := src/fast_psf_fftw.c +BUILD_CPPFLAGS += $(CPU_FFTW_CPPFLAGS) +LDLIBS += $(CPU_FFTW_LDLIBS) +else +$(error PSF_BACKEND=cpu requires the FFTW OpenMP package 'fftw3_omp' (pkg-config --exists fftw3_omp). On Gentoo enable sci-libs/fftw[openmp].) +endif else ifeq ($(PSF_BACKEND),hip) RENDER_LINKER := $(HIPCC) BUILD_CPPFLAGS += -DPSF_BACKEND_HIP @@ -63,6 +64,23 @@ else $(error Unknown PSF_BACKEND '$(PSF_BACKEND)'; choose cpu, hip, or dummy) endif +# `make test` builds and runs the CPU regression suite. The dummy and hip +# backends do not build those binaries, so reject the goal explicitly instead of +# silently reusing a stale test binary left by an earlier CPU build. +ifneq ($(filter test test-reference-images,$(MAKECMDGOALS)),) +ifneq ($(PSF_BACKEND),cpu) +$(error make test requires PSF_BACKEND=cpu) +endif +endif + +# psf-capture is a CPU producer diagnostic; the dummy and hip backends do not +# provide the producer/frame symbols it links. +ifneq ($(filter psf-capture,$(MAKECMDGOALS)),) +ifneq ($(PSF_BACKEND),cpu) +$(error psf-capture requires PSF_BACKEND=cpu) +endif +endif + ifeq ($(SPACETIME),minkowski) BACKEND_CPPFLAGS := -DSPACETIME_MINKOWSKI else ifeq ($(SPACETIME),schwarzschild) @@ -81,6 +99,21 @@ else $(error Unknown ENABLE_HDR '$(ENABLE_HDR)'; choose 0 or 1) endif +# Test binaries are namespaced by spacetime and build tag so that switching +# PSF_BACKEND/ENABLE_HDR/PSF_EVENT_SINK cannot silently reuse a stale binary +# built for a different backend. +TEST_OUT_DIR := $(OBJECT_DIR)/$(HDR_BUILD_TAG) +TEST_TARGET := $(TEST_OUT_DIR)/test_geodesic +FRAME_TEST_TARGET := $(TEST_OUT_DIR)/test_frame +SCHWARZSCHILD_TEST_TARGET := $(TEST_OUT_DIR)/test_schwarzschild +OBSERVER_TRACK_TEST_TARGET := $(TEST_OUT_DIR)/test_observer_track +CATALOG_PREFETCH_TEST_TARGET := $(TEST_OUT_DIR)/test_catalog_prefetch +HIP_PSF_TEST_TARGET := $(TEST_OUT_DIR)/test_hip_psf +HIP_PSF_BENCH_TARGET := $(TEST_OUT_DIR)/benchmark_hip_psf +CAMERA_TEST_TARGETS := $(TEST_OUT_DIR)/test_observer_minkowski $(TEST_OUT_DIR)/test_observer_schwarzschild +FAST_PSF_FFTW_TEST_TARGET := $(TEST_OUT_DIR)/test_fast_psf_fftw +FAST_PSF_FFTW_BENCH_TARGET := $(TEST_OUT_DIR)/benchmark_fast_psf_fftw + ifeq ($(PSF_BACKEND),hip) TARGET := $(BUILD_DIR)/$(TARGET_BASENAME)_hip else ifeq ($(PSF_BACKEND),dummy) @@ -89,6 +122,7 @@ else TARGET := $(BUILD_DIR)/$(TARGET_BASENAME) endif RENDER_SOURCES := $(COMMON_SOURCES) $(PROVIDER_SOURCE) src/main.c +RENDER_SOURCES += $(CPU_FFTW_SOURCES) ifeq ($(PSF_BACKEND),dummy) RENDER_SOURCES += src/dummy_psf.c endif @@ -140,10 +174,10 @@ ifeq ($(PSF_BACKEND),hip) hip-psf-test: $(HIP_PSF_TEST_TARGET) hip-psf-bench: $(HIP_PSF_BENCH_TARGET) -$(HIP_PSF_BENCH_TARGET): tests/benchmark_hip_psf.hip $(HIP_PSF_OBJECT) $(OBJECT_DIR)/$(HDR_BUILD_TAG)/src/optics.o | $(BUILD_DIR) +$(HIP_PSF_BENCH_TARGET): tests/benchmark_hip_psf.hip $(HIP_PSF_OBJECT) $(OBJECT_DIR)/$(HDR_BUILD_TAG)/src/optics.o | $(TEST_OUT_DIR) $(HIPCC) $(HIP_CXXFLAGS) $(OPENMP_FLAGS) -Isrc -x hip $< -x none $(filter-out $<,$^) $(LDLIBS) $(HDR_LDLIBS) -o $@ -$(HIP_PSF_TEST_TARGET): tests/test_hip_psf.hip $(HIP_PSF_OBJECT) $(OBJECT_DIR)/$(HDR_BUILD_TAG)/src/optics.o | $(BUILD_DIR) +$(HIP_PSF_TEST_TARGET): tests/test_hip_psf.hip $(HIP_PSF_OBJECT) $(OBJECT_DIR)/$(HDR_BUILD_TAG)/src/optics.o | $(TEST_OUT_DIR) $(HIPCC) $(HIP_CXXFLAGS) $(OPENMP_FLAGS) -Isrc -x hip $< -x none $(filter-out $<,$^) $(LDLIBS) $(HDR_LDLIBS) -o $@ else hip-psf-test hip-psf-bench: @@ -154,36 +188,57 @@ run: $(TARGET) mkdir -p output/imgs ./$(TARGET) --catalog assets/sky_grid_5deg.csv --output output/imgs/$(SPACETIME)_sky.$(IMAGE_EXT) -$(TEST_TARGET): tests/test_geodesic.c $(CORE_MINKOWSKI_SOURCES) | $(BUILD_DIR) +$(TEST_OUT_DIR): | $(BUILD_DIR) + mkdir -p $@ + +$(TEST_TARGET): tests/test_geodesic.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ -$(FRAME_TEST_TARGET): tests/test_frame.c $(CORE_MINKOWSKI_SOURCES) | $(BUILD_DIR) +$(FRAME_TEST_TARGET): tests/test_frame.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ -$(SCHWARZSCHILD_TEST_TARGET): tests/test_schwarzschild.c $(COMMON_SOURCES) src/spacetime_schwarzschild.c | $(BUILD_DIR) +$(SCHWARZSCHILD_TEST_TARGET): tests/test_schwarzschild.c $(COMMON_SOURCES) src/spacetime_schwarzschild.c $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ -$(OBSERVER_TRACK_TEST_TARGET): tests/test_observer_track.c $(COMMON_SOURCES) | $(BUILD_DIR) +$(OBSERVER_TRACK_TEST_TARGET): tests/test_observer_track.c $(COMMON_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ -$(CATALOG_PREFETCH_TEST_TARGET): tests/test_catalog_prefetch.c $(CORE_MINKOWSKI_SOURCES) | $(BUILD_DIR) +$(CATALOG_PREFETCH_TEST_TARGET): tests/test_catalog_prefetch.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ -$(BUILD_DIR)/test_observer_minkowski: tests/test_observer.c $(CORE_MINKOWSKI_SOURCES) | $(BUILD_DIR) +$(TEST_OUT_DIR)/test_observer_minkowski: tests/test_observer.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ -$(BUILD_DIR)/test_observer_schwarzschild: tests/test_observer.c $(COMMON_SOURCES) src/spacetime_schwarzschild.c | $(BUILD_DIR) +$(TEST_OUT_DIR)/test_observer_schwarzschild: tests/test_observer.c $(COMMON_SOURCES) src/spacetime_schwarzschild.c $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -DSPACETIME_SCHWARZSCHILD -Isrc $^ $(LDLIBS) -o $@ -test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) - ./$(BUILD_DIR)/test_observer_minkowski - ./$(BUILD_DIR)/test_observer_schwarzschild - ./$(TEST_TARGET) - ./$(FRAME_TEST_TARGET) - ./$(SCHWARZSCHILD_TEST_TARGET) - ./$(OBSERVER_TRACK_TEST_TARGET) - ./$(CATALOG_PREFETCH_TEST_TARGET) - python3 tests/test_camera_cli.py $(BUILD_DIR) +$(FAST_PSF_FFTW_TEST_TARGET): tests/test_fast_psf_fftw.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ + +$(FAST_PSF_FFTW_BENCH_TARGET): tests/benchmark_fast_psf_fftw.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ + +# The FFTW-vs-spatial test is meaningful only in the CPU PSF build. +ifneq ($(CPU_FFTW_SOURCES),) +FAST_PSF_FFTW_TEST_DEP := $(FAST_PSF_FFTW_TEST_TARGET) +FAST_PSF_FFTW_TEST_RUN := $(FAST_PSF_FFTW_TEST_TARGET) +else +FAST_PSF_FFTW_TEST_DEP := +FAST_PSF_FFTW_TEST_RUN := +endif + +test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) + $(TEST_OUT_DIR)/test_observer_minkowski + $(TEST_OUT_DIR)/test_observer_schwarzschild + $(TEST_TARGET) + $(FRAME_TEST_TARGET) + $(SCHWARZSCHILD_TEST_TARGET) + $(OBSERVER_TRACK_TEST_TARGET) + $(CATALOG_PREFETCH_TEST_TARGET) + $(FAST_PSF_FFTW_TEST_RUN) + python3 tests/test_camera_cli.py $(BUILD_DIR) $(TEST_OUT_DIR) + +fast-psf-fftw-bench: $(FAST_PSF_FFTW_BENCH_TARGET) clean: rm -rf build @@ -193,14 +248,14 @@ include mk/reference_images.mk # Test-only producer consumer: never linked into a renderer. .PHONY: psf-capture -psf-capture: $(BUILD_DIR)/capture_psf -$(BUILD_DIR)/capture_psf: tests/capture_psf.c $(CORE_MINKOWSKI_SOURCES) | $(BUILD_DIR) - $(CC) $(CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $< $(filter-out src/frame.c,$(CORE_MINKOWSKI_SOURCES)) $(LDLIBS) -o $@ +psf-capture: $(TEST_OUT_DIR)/capture_psf +$(TEST_OUT_DIR)/capture_psf: tests/capture_psf.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $< $(filter-out src/frame.c,$(CORE_MINKOWSKI_SOURCES)) $(CPU_FFTW_SOURCES) $(LDLIBS) -o $@ .PHONY: hip-psf-replay -hip-psf-replay: $(BUILD_DIR)/replay_psf -$(BUILD_DIR)/replay_psf: tests/replay_psf.hip src/hip_psf.hip src/hip_psf.h src/optics.h $(OBJECT_DIR)/$(HDR_BUILD_TAG)/src/optics.o | $(BUILD_DIR) +hip-psf-replay: $(TEST_OUT_DIR)/replay_psf +$(TEST_OUT_DIR)/replay_psf: tests/replay_psf.hip src/hip_psf.hip src/hip_psf.h src/optics.h $(OBJECT_DIR)/$(HDR_BUILD_TAG)/src/optics.o | $(TEST_OUT_DIR) $(HIPCC) $(HIP_CXXFLAGS) $(OPENMP_FLAGS) -Isrc -x hip $< -x none $(filter %.o,$^) $(LDLIBS) $(HDR_LDLIBS) -o $@ -$(BUILD_DIR)/make_psf_fixture: tests/make_psf_fixture.c src/optics.c src/optics.h | $(BUILD_DIR) - $(CC) $(CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $< src/optics.c $(LDLIBS) -o $@ +$(TEST_OUT_DIR)/make_psf_fixture: tests/make_psf_fixture.c src/optics.c src/optics.h $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $< src/optics.c $(CPU_FFTW_SOURCES) $(LDLIBS) -o $@ diff --git a/benchmarks/fast_mode_fftw_2026-09-25.md b/benchmarks/fast_mode_fftw_2026-09-25.md new file mode 100644 index 0000000..2464b7c --- /dev/null +++ b/benchmarks/fast_mode_fftw_2026-09-25.md @@ -0,0 +1,280 @@ +# Fast-mode FFTW convolution benchmark (2026-09-25) + +Replaces the spatial fast-mode global convolution with a reusable FFTW linear +convolution. This record separates one-time plan/setup from steady-state +execution and repeats the bounded Galactic-center comparison. Raw resolver and +end-to-end output is pasted verbatim from the runs below. + +## Environment + +- Git revision: `3deebfb` plus the uncommitted FFTW work. +- Host: 12th Gen Intel Core i7-12700K, 16 hardware threads, 128 GB. +- FFTW: `pkg-config --modversion fftw3_omp` = `3.3.10`; + `pkg-config --libs fftw3_omp` = `-lfftw3_omp -lfftw3`. +- Build: + `make -j8 BUILD_TYPE=Release PSF_BACKEND=cpu SPACETIME=minkowski backend` + (`-O2 -DNDEBUG -fopenmp -march=native`, libpng, `-DFAST_PSF_FFTW`). +- Binary SHA-256 (this run): + - `build/Release/minkowski_sky`: + `ead2a543b8166abe6f82bc1cd41b2e985a6db517aa45a6416d08d188138b8c47` + - `build/Release/obj/minkowski/standard_sink1_cpu/benchmark_fast_psf_fftw`: + `8da8bbd532104a0c4f524aabf40b66a50c8c715e44479c7db467669ab5a21215` +- Timing wrapper: `/usr/bin/time` is not installed on this host. The + reproducible wrapper used for the end-to-end runs is zsh's built-in `time` + with `TIMEFMT='real=%*E user=%U sys=%S cpu=%P'`, i.e. each command was run as + `TIMEFMT='real=%*E user=%U sys=%S cpu=%P'; time `. The `real=...` + lines below are the verbatim wrapper output. + +## Correctness + +```sh +make -j8 BUILD_TYPE=Debug PSF_BACKEND=cpu SPACETIME=minkowski test +``` + +All C tests pass, including the dedicated FFTW-vs-spatial test (258 comparisons +over sizes `17x13`/`16x16`/`23x31`/`33x17`, `N = 1..4`, nearest and bilinear, +eight scenes; worst `max_abs = 1.8e-14`), the `test_frame` fast-mode +regressions, the non-fast FITS references, and both camera CLI checks. HIP and +dummy renderers compile and link without FFTW; `make PSF_BACKEND=dummy test` +now fails fast with `make test requires PSF_BACKEND=cpu` instead of reusing a +stale CPU test binary. + +## Resolver-only benchmark + +One deterministic sparse impulse buffer, default PSF (FWHM 2.7, beta 4.5, +relative tail 1e-8), `OMP_NUM_THREADS=16`. `--spatial` additionally runs the +reference resolver on the identical buffer. Binary: +`build/Release/obj/minkowski/standard_sink1_cpu/benchmark_fast_psf_fftw`. + +### Default 320x240 N=2, FFTW_ESTIMATE + +```sh +OMP_NUM_THREADS=16 ./build/Release/obj/minkowski/standard_sink1_cpu/benchmark_fast_psf_fftw \ + --width 320 --height 240 --supersample 2 --repeats 7 --spatial +``` + +```text +impulse: 320x240 supersample=2 deposit=nearest repeats=7 spatial=1 +plan_mode=estimate +setup: init=0.039444 s fftw_plan=0.008085 s kernel_fft=0.005237 s +Fast FFTW: linear min=666x826, fft=672x840, workers=16, plan=estimate, plan=0.008085 s, kernel_fft=0.005237 s, scratch=30.2 MiB +impulse_hash=5334fa90cfbb1089 +fftw: min=0.005602 s median=0.006455 s max=0.008144 s +fftw samples: 0.006697 0.005784 0.005602 0.005898 0.006455 0.007675 0.008144 +fftw_hdr_hash=9ad60d92992c9bbd +spatial: 0.702470 s +spatial_hdr_hash=9327ac3e8e10d670 +difference: max_abs=7.49e-16 peak=0.359 speedup=108.831x +rss_max_kb=52908 +``` + +### Default 320x240 N=4, FFTW_ESTIMATE + +```sh +OMP_NUM_THREADS=16 ./build/Release/obj/minkowski/standard_sink1_cpu/benchmark_fast_psf_fftw \ + --width 320 --height 240 --supersample 4 --repeats 7 --spatial +``` + +```text +impulse: 320x240 supersample=4 deposit=nearest repeats=7 spatial=1 +plan_mode=estimate +setup: init=0.114719 s fftw_plan=0.007169 s kernel_fft=0.011942 s +Fast FFTW: linear min=1330x1650, fft=1344x1680, workers=16, plan=estimate, plan=0.007169 s, kernel_fft=0.011942 s, scratch=120.7 MiB +impulse_hash=2639864ea9d3a26a +fftw: min=0.034756 s median=0.039591 s max=0.042474 s +fftw samples: 0.041352 0.039591 0.034756 0.040866 0.042474 0.035364 0.036241 +fftw_hdr_hash=1ff070fdaac85f4a +spatial: 11.054385 s +spatial_hdr_hash=7682ee8550f88efe +difference: max_abs=5.38e-15 peak=0.533 speedup=279.214x +rss_max_kb=161032 +``` + +### Medium 960x540 N=2, FFTW_ESTIMATE (FFTW only) + +```sh +OMP_NUM_THREADS=16 ./build/Release/obj/minkowski/standard_sink1_cpu/benchmark_fast_psf_fftw \ + --width 960 --height 540 --supersample 2 --repeats 7 +``` + +```text +impulse: 960x540 supersample=2 deposit=nearest repeats=7 spatial=0 +plan_mode=estimate +setup: init=0.042270 s fftw_plan=0.006317 s kernel_fft=0.010302 s +Fast FFTW: linear min=1266x2106, fft=1280x2160, workers=16, plan=estimate, plan=0.006317 s, kernel_fft=0.010302 s, scratch=147.7 MiB +impulse_hash=6e2de17f6859cd33 +fftw: min=0.038361 s median=0.039872 s max=0.055834 s +fftw samples: 0.055834 0.038745 0.039872 0.038823 0.038361 0.048436 0.040141 +fftw_hdr_hash=c493dbb793d04706 +rss_max_kb=216272 +``` + +### 4K 3840x2160 N=2, FFTW_ESTIMATE (FFTW only) + +```sh +OMP_NUM_THREADS=16 ./build/Release/obj/minkowski/standard_sink1_cpu/benchmark_fast_psf_fftw \ + --width 3840 --height 2160 --supersample 2 --repeats 5 +``` + +```text +impulse: 3840x2160 supersample=2 deposit=nearest repeats=5 spatial=0 +plan_mode=estimate +setup: init=0.524840 s fftw_plan=0.008154 s kernel_fft=0.482906 s +Fast FFTW: linear min=4506x7866, fft=4536x7875, workers=16, plan=estimate, plan=0.008154 s, kernel_fft=0.482906 s, scratch=1907.8 MiB +impulse_hash=e3b6d7f3d621630d +fftw: min=2.368634 s median=2.481927 s max=2.553895 s +fftw samples: 2.552212 2.553895 2.481927 2.389895 2.368634 +fftw_hdr_hash=1c5720f72a040e6d +rss_max_kb=2879596 +``` + +### 4K 3840x2160 N=4, FFTW_ESTIMATE (FFTW only) + +```sh +OMP_NUM_THREADS=16 ./build/Release/obj/minkowski/standard_sink1_cpu/benchmark_fast_psf_fftw \ + --width 3840 --height 2160 --supersample 4 --repeats 5 +``` + +```text +impulse: 3840x2160 supersample=4 deposit=nearest repeats=5 spatial=0 +plan_mode=estimate +setup: init=0.738771 s fftw_plan=0.003636 s kernel_fft=0.599235 s +Fast FFTW: linear min=9010x15730, fft=9072x15750, workers=16, plan=estimate, plan=0.003636 s, kernel_fft=0.599235 s, scratch=7631.4 MiB +impulse_hash=e6a528db0fb16c33 +fftw: min=3.276609 s median=3.339285 s max=3.369411 s +fftw samples: 3.309935 3.369411 3.276609 3.345614 3.339285 +fftw_hdr_hash=e55bed24ef3f5506 +rss_max_kb=10919208 +``` + +### FFTW_ESTIMATE versus FFTW_MEASURE + +`FFTW_MEASURE` cuts steady-state time but costs a large one-time plan. Because +the plan is created once per process and reused across frames, it only pays off +for long movie renders; production therefore defaults to `FFTW_ESTIMATE`. +Wisdom import/export would make `FFTW_MEASURE` the better default and remains a +follow-up. Measured FFTW_MEASURE plan output is not bit-reproducible (plan +measurement is timing-dependent), so its `fftw_hdr_hash` may differ run to run; +the spatial comparison stays at roundoff. + +```sh +OMP_NUM_THREADS=16 ./build/Release/obj/minkowski/standard_sink1_cpu/benchmark_fast_psf_fftw \ + --width 320 --height 240 --supersample 2 --repeats 7 --measure --spatial +``` + +```text +impulse: 320x240 supersample=2 deposit=nearest repeats=7 spatial=1 +plan_mode=measure +setup: init=1.584388 s fftw_plan=1.091687 s kernel_fft=0.467405 s +Fast FFTW: linear min=666x826, fft=672x840, workers=16, plan=measure, plan=1.091687 s, kernel_fft=0.467405 s, scratch=30.2 MiB +impulse_hash=5334fa90cfbb1089 +fftw: min=0.003995 s median=0.004702 s max=0.005804 s +fftw samples: 0.003995 0.004702 0.004330 0.005804 0.004939 0.004986 0.004619 +fftw_hdr_hash=626e386cba92ffa5 +spatial: 0.706244 s +spatial_hdr_hash=9327ac3e8e10d670 +difference: max_abs=8.05e-16 peak=0.359 speedup=150.198x +rss_max_kb=51860 +``` + +```sh +OMP_NUM_THREADS=16 ./build/Release/obj/minkowski/standard_sink1_cpu/benchmark_fast_psf_fftw \ + --width 3840 --height 2160 --supersample 2 --repeats 5 --measure +``` + +```text +impulse: 3840x2160 supersample=2 deposit=nearest repeats=5 spatial=0 +plan_mode=measure +setup: init=41.408673 s fftw_plan=41.239987 s kernel_fft=0.134272 s +Fast FFTW: linear min=4506x7866, fft=4536x7875, workers=16, plan=measure, plan=41.239987 s, kernel_fft=0.134272 s, scratch=1907.8 MiB +impulse_hash=e3b6d7f3d621630d +fftw: min=0.518707 s median=0.522181 s max=0.539465 s +fftw samples: 0.523380 0.539465 0.518707 0.520357 0.522181 +fftw_hdr_hash=ab90c4b076c775e1 +rss_max_kb=2884132 +``` + +## End-to-end bounded comparison + +Same catalog, geometry, exposure, supersample factors, and worker environment as +`benchmarks/fast_mode_cpu_2026-09-18.md`. Raw terminal output including the +timing wrapper: + +```sh +OMP_NUM_THREADS=16 TIMEFMT='real=%*E user=%U sys=%S cpu=%P' \ + time ./build/Release/minkowski_sky --all-sky-catalog assets/2mass/processed/all_sky \ + --width 320 --height 240 --fov-deg 10 --look-ra-deg 266 --look-dec-deg -29 \ + --exposure 1e13 --coarse-cell-pixels 16 --refine-max-level 0 \ + --output /tmp/opencode/fastmode_fftw_bench/normal.png +``` + +```text +Blackbody LUT: 1024 CIE 1931 2-deg 1 nm-linear XYZ nodes, T=[670.146556, 101408.88] K, loaded assets/blackbody/cie1931_2deg_xyz_1024.grbblut in 0.000046 s; payload fnv1a64=0e39867b0e70a809 +Blackbody backend: lut +PSF cache ready: 64x64 phases, radius 47 px, relative tail 1e-08, tail abs 1e-06, boundary 1e-07, build 0.776 s +Catalog prefetch: 96 requested, 96 newly loaded (6618807 stars), 0 unavailable in 0.729 s; 4 loader workers +Rendered 5912840 images from 6618807 catalog stars to /tmp/opencode/fastmode_fftw_bench/normal.png (ok) +PSF splats: cached 5912840, cached wing-clipped 0, direct fallbacks 0, discarded below min-Y 0 +real=17.526 user=225.46s sys=2.44s cpu=1300% +``` + +```sh +OMP_NUM_THREADS=16 TIMEFMT='real=%*E user=%U sys=%S cpu=%P' \ + time ./build/Release/minkowski_sky --all-sky-catalog assets/2mass/processed/all_sky \ + --width 320 --height 240 --fov-deg 10 --look-ra-deg 266 --look-dec-deg -29 \ + --exposure 1e13 --coarse-cell-pixels 16 --refine-max-level 0 \ + --fast-mode --fast-supersample 2 \ + --output /tmp/opencode/fastmode_fftw_bench/fast_n2.png +``` + +```text +Fast mode is a preview approximation; --max-cache-psf-flux is ignored and one global kernel is used for every event. +Blackbody LUT: 1024 CIE 1931 2-deg 1 nm-linear XYZ nodes, T=[670.146556, 101408.88] K, loaded assets/blackbody/cie1931_2deg_xyz_1024.grbblut in 0.000049 s; payload fnv1a64=0e39867b0e70a809 +Blackbody backend: lut +Fast PSF: supersample 2x, deposit nearest, kernel radius 93 ss px (46.50 final px), retained flux 1.000000, cached 4-point quadrature +Fast FFTW: linear min=666x826, fft=672x840, workers=16, plan=estimate, plan=0.007770 s, kernel_fft=0.005020 s, scratch=30.2 MiB +Catalog prefetch: 96 requested, 96 newly loaded (6618807 stars), 0 unavailable in 0.735 s; 4 loader workers +Rendered 5912840 images from 6618807 catalog stars to /tmp/opencode/fastmode_fftw_bench/fast_n2.png (ok) +Fast PSF splats: deposited 5912840, wing-clipped 0, discarded below min-Y 0 +real=0.995 user=5.30s sys=0.18s cpu=550% +``` + +```sh +OMP_NUM_THREADS=16 TIMEFMT='real=%*E user=%U sys=%S cpu=%P' \ + time ./build/Release/minkowski_sky --all-sky-catalog assets/2mass/processed/all_sky \ + --width 320 --height 240 --fov-deg 10 --look-ra-deg 266 --look-dec-deg -29 \ + --exposure 1e13 --coarse-cell-pixels 16 --refine-max-level 0 \ + --fast-mode --fast-supersample 4 \ + --output /tmp/opencode/fastmode_fftw_bench/fast_n4.png +``` + +```text +Fast mode is a preview approximation; --max-cache-psf-flux is ignored and one global kernel is used for every event. +Blackbody LUT: 1024 CIE 1931 2-deg 1 nm-linear XYZ nodes, T=[670.146556, 101408.88] K, loaded assets/blackbody/cie1931_2deg_xyz_1024.grbblut in 0.000044 s; payload fnv1a64=0e39867b0e70a809 +Blackbody backend: lut +Fast PSF: supersample 4x, deposit nearest, kernel radius 185 ss px (46.25 final px), retained flux 1.000000, cached 4-point quadrature +Fast FFTW: linear min=1330x1650, fft=1344x1680, workers=16, plan=estimate, plan=0.007294 s, kernel_fft=0.012307 s, scratch=120.7 MiB +Catalog prefetch: 96 requested, 96 newly loaded (6618807 stars), 0 unavailable in 0.722 s; 4 loader workers +Rendered 5912840 images from 6618807 catalog stars to /tmp/opencode/fastmode_fftw_bench/fast_n4.png (ok) +Fast PSF splats: deposited 5912840, wing-clipped 0, discarded below min-Y 0 +real=1.033 user=5.56s sys=0.22s cpu=559% +``` + +Image counts are identical across normal and both fast runs (5,912,840), so the +FFTW resolver did not change catalog/deposit behavior. Output SHA-256: + +```text +normal.png e9f0a98ae7a3c3db546a3f276f400be98c96446fab7bfbc632f388702d4722d1 +fast_n2.png 1f144aefb4a5520d9d76a13ab69026a72580663e5e45f53128a87a0801176ccf +fast_n4.png 679f14a67a085f76c144ecbe36af51b3b2534c152e8ef59708a7884a3f7c8464 +``` + +## Reading + +- Steady-state FFTW is ~109x (N=2) and ~279x (N=4) faster than the spatial + reference at 320x240, and stays under ~0.06 s per frame at 960x540. +- 4K N=2 resolves in ~2.5 s (ESTIMATE) or ~0.52 s (MEASURE) on this host; 4K N=4 + needs ~7.6 GiB scratch and ~3.3 s (ESTIMATE). +- FFTW-versus-spatial differences are floating-point roundoff (`max_abs` at the + `1e-15` level against peaks near `0.5`), not crop, wrap, channel, or + normalization errors. diff --git a/build.md b/build.md index dc55da9..1ffc96e 100644 --- a/build.md +++ b/build.md @@ -10,9 +10,13 @@ The default CPU build needs: - A C11 compiler and OpenMP runtime (for example, GCC with libgomp). - GNU Make. - libpng headers and library for PNG output. +- `pkg-config`, used to locate and link FFTW. +- FFTW 3 built with OpenMP, discovered through `pkg-config fftw3_omp` + (`sci-libs/fftw[openmp]` on Gentoo). This is required for the CPU fast-mode + fast convolution; HIP and dummy PSF builds do not link FFTW. -Optional dependencies are CFITSIO and `pkg-config` for HDR/FITS output, and -HIP/ROCm with `hipcc` and a compatible GPU for HIP PSF accumulation. +Optional dependencies are CFITSIO for HDR/FITS output, and HIP/ROCm with +`hipcc` and a compatible GPU for HIP PSF accumulation. ## CPU build @@ -181,7 +185,7 @@ the test executable; run it separately: ```sh make PSF_BACKEND=hip SPACETIME=minkowski hip-psf-test -./build/Release/test_hip_psf +./build/Release/obj/minkowski/standard_sink1_hip/test_hip_psf ``` For rendering with real stellar data, continue with [README.md](README.md#prepare-the-stellar-catalog) diff --git a/fast_mode_fftw_optimization_plan.md b/fast_mode_fftw_optimization_plan.md new file mode 100644 index 0000000..18e63e2 --- /dev/null +++ b/fast_mode_fftw_optimization_plan.md @@ -0,0 +1,799 @@ +# Fast-mode FFTW convolution optimization plan + +## 1. Status and objective + +This plan starts from commit `3deebfb` (`Feat: Add --fast-mode supersampled +point-source accumulation`). The current CPU fast mode already: + +- maps catalog images to continuous image positions; +- deposits each image into one nearest supersampled cell by default, or four + bilinear cells when explicitly requested; +- stores a supersampled double-RGB impulse buffer; +- applies one global pixel-integrated Moffat kernel; +- averages each `N x N` supersampled block into the final HDR framebuffer; +- preserves additive HDR semantics, `--psf-min-y`, wing-clipping statistics, + and the documented preview-only boundary. + +The remaining bottleneck is `fast_psf_accumulator_resolve()`, which evaluates +the global convolution by nested spatial loops. The earlier HIP work already +established PSF accumulation as the dominant stage, and the target machine is +the local 128 GB host. FFTW is an accepted dependency. This work therefore +does **not** repeat a feasibility study or reconsider the selected nearest-cell +visual tradeoff. Its objective is: + +> Replace the spatial global convolution with a reusable, multithreaded FFTW +> linear-convolution path while reproducing the current discrete fast-mode +> result to floating-point tolerance and retaining the spatial implementation +> as a test/benchmark reference. + +The target is the CPU PSF backend. HIP and dummy PSF builds continue to reject +`--fast-mode` and must not acquire an unnecessary FFTW link dependency. + +## 2. Fixed requirements and non-goals + +### 2.1 Required behavior + +- Use double-precision FFTW (`fftw3`), not the float API. +- Use the locally available FFTW OpenMP backend (`fftw3_omp`); the current host + reports FFTW `3.3.10` and `pkg-config --libs fftw3_omp` returns + `-lfftw3_omp -lfftw3`. +- On Gentoo this requires `sci-libs/fftw[openmp]`, which is already enabled on + the target host. It does **not** require the separate `threads` USE flag: + `threads` selects FFTW's pthread backend, while this plan deliberately uses + the OpenMP backend so it shares the renderer's OpenMP runtime. +- Compute zero-padded **linear** convolution. Circular wraparound is a failure. +- Match the current circular kernel traversal, not the unused square corners of + the allocated `weights` array. +- Preserve the current impulse-buffer contents so the same buffer can be sent + through the spatial reference and FFTW implementations in one process. +- Normalize the unscaled FFTW inverse exactly once. +- Apply the existing `1/N^2` box-average normalization exactly once. +- Crop the full linear convolution at the mathematically correct offset before + downsampling. +- Add the result to the caller's HDR buffer; never overwrite existing HDR. +- Cache plans, the kernel spectrum, FFT buffers, and dimension metadata across + movie frames. +- Keep FFTW planning outside OpenMP producer regions and avoid nested + OpenMP/FFTW oversubscription. +- Report one-time planning/kernel-transform time separately from per-frame FFT + execution time. +- Propagate allocation or plan failures as explicit fast-mode initialization + failures. Do not silently fall back to a potentially multi-minute spatial + convolution. + +### 2.2 Preserved preview semantics + +This optimization must not change: + +- nearest as the default deposit mode; +- bilinear as the explicit position-over-shape alternative; +- supersample factor or coordinate convention; +- the Moffat `FWHM`, `beta`, quadrature, normalization, or support radius; +- flux, color, magnification, frequency-shift, exposure, or `min-Y` decisions; +- fixed global-kernel wing clipping and its statistics; +- zero-outside-image behavior, including the existing narrow bilinear edge + approximation; +- tone mapping, PNG/FITS output, lens-map reuse, or catalog traversal. + +### 2.3 Non-goals + +- Do not redesign the producer/catalog path. +- Do not introduce a GPU FFT backend in this change. +- Do not add polyphase convolution as a competing production algorithm. +- Do not change nearest/bilinear acceptance criteria. +- Do not add an artificial memory cap for the 128 GB target host. +- Do not make bitwise identity an acceptance criterion; FFT summation order is + different from direct spatial accumulation. +- Do not combine this work with impulse-buffer persistence or PSF/exposure sweep + orchestration. The internal state should permit those later, but they are a + separate feature. + +## 3. Exact discrete operation to preserve + +Let the supersampled impulse buffer be + +```text +B[c][y][x], c in {R,G,B}, 0 <= x < Wss, 0 <= y < Hss +``` + +with `Wss = N * width` and `Hss = N * height`. Let the stored kernel radius be +`R`, its side be `K = 2R + 1`, and its circularly retained weights be + +```text +A[dy + R][dx + R] = weights[(dy + R) * K + (dx + R)] +``` + +only where + +```text +-R <= dy <= R +|dx| <= floor(sqrt(R^2 - dy^2)). +``` + +All other elements of `A` are zero. The current fine-grid result is + +```text +C[c][sy][sx] = sum_dy sum_dx A[dy + R][dx + R] + * B[c][sy - dy][sx - dx], +``` + +where samples of `B` outside its image are zero. The final operation is + +```text +hdr[c][row][column] += + (1 / N^2) * sum_j=0..N-1 sum_i=0..N-1 + C[c][N*row + j][N*column + i]. +``` + +The FFTW implementation must preserve this discrete operation, including the +current float kernel weights promoted to double during multiplication. + +## 4. Linear-convolution layout and crop convention + +This is the most error-prone part of the implementation and must be fixed by +tests before performance work. + +### 4.1 Padding bounds + +Choose FFT dimensions satisfying + +```text +Pwidth >= Wss + K - 1 = Wss + 2R +Pheight >= Hss + K - 1 = Hss + 2R. +``` + +`B` is copied at padded origin `(0,0)`. The kernel is also stored at padded +origin with its array indices unchanged: + +```text +padded_kernel[dy + R][dx + R] = A[dy + R][dx + R]. +``` + +Do **not** apply an additional `fftshift`/`ifftshift` in this convention. The +ordinary full convolution produced by FFT multiplication has dimensions +`Wss + K - 1` by `Hss + K - 1`, and the desired same-size fine-grid result is + +```text +C[sy][sx] = full_convolution[sy + R][sx + R]. +``` + +Thus the fine-grid crop begins exactly at `(R,R)` and has extent +`Wss x Hss`. + +An alternative center-at-frequency-origin layout is mathematically possible, +but it changes the crop/wrap convention and is deliberately excluded from the +first implementation. + +### 4.2 FFT-friendly dimensions + +Add a checked helper such as `next_smooth_size()` that returns the first value +at least the required extent whose prime factors are limited to small FFTW- +friendly factors, initially `{2,3,5,7}`. The helper must: + +- accept and return `size_t` while checking overflow; +- reject a selected dimension above `INT_MAX` before passing it to FFTW's + `int n[2]` interface; +- have unit tests for exact, next-size, and near-overflow inputs; +- record both the minimum linear-convolution extent and selected FFT extent. + +Do not hard-code powers of two: they can waste substantially more memory and +work than a nearby smooth composite size. + +## 5. Source and module boundaries + +### 5.1 New private FFTW module + +Keep FFTW details out of the public optics API and out of unrelated backends. +Add: + +```text +src/fast_psf_fftw.c +src/fast_psf_fftw.h +``` + +`fast_psf_fftw.h` is a private renderer header. It declares an opaque state +and functions conceptually equivalent to: + +```c +typedef struct FastPsfFftwState FastPsfFftwState; + +FastPsfFftwState *fast_psf_fftw_create( + int ss_width, int ss_height, + int final_width, int final_height, + int supersample, int kernel_radius, + const float *weights, + int fft_workers, + FastPsfFftwTiming *timing); + +int fast_psf_fftw_resolve( + FastPsfFftwState *state, + const double *interleaved_ss_rgb, + double *interleaved_hdr_rgb, + FastPsfFftwTiming *timing); + +void fast_psf_fftw_destroy(FastPsfFftwState *state); +``` + +Exact names may follow project style, but ownership must stay explicit: + +- `FastPsfAccumulator` owns one FFTW state pointer; +- the FFTW state owns every FFTW allocation and plan; +- the input impulse buffer remains owned by `FastPsfAccumulator`; +- the output HDR remains caller-owned; +- the kernel spectrum is immutable after initialization; +- destruction is safe for a zero/partially initialized state. + +Avoid exposing `fftw_plan` or `fftw_complex` through `optics.h`. If the public +`FastPsfAccumulator` must carry the private pointer, use a forward-declared +opaque struct. + +### 5.2 Spatial reference path + +Extract the current nested-loop resolver into a clearly named reference helper, +for example: + +```c +fast_psf_accumulator_resolve_spatial_reference(...) +``` + +It must remain callable by a dedicated FFTW test/benchmark with the identical +impulse buffer. Production `fast_psf_accumulator_resolve()` dispatches to +FFTW in CPU builds after validation. The spatial path is not a silent runtime +fallback. + +The reference helper should reuse a cached circular `row_span` array rather +than allocating it on every call. This small cleanup also guarantees that the +spatial and FFTW kernel masks are built from one definition. + +## 6. FFTW state and memory layout + +### 6.1 Planar FFT scratch + +Keep the existing interleaved impulse buffer because the deposit hot path and +its atomic RGB updates already use that layout. At resolve time, pack it into +FFTW-aligned planar scratch: + +```text +real_rgb[channel][padded_y][padded_x] +frequency_rgb[channel][ky][kx], kx extent = Pwidth/2 + 1 +kernel_frequency[ky][kx] +``` + +Use: + +- `fftw_alloc_real()` for real scratch; +- `fftw_alloc_complex()` for channel spectra and the kernel spectrum; +- out-of-place R2C and C2R transforms in the first implementation; +- `fftw_plan_many_dft_r2c()` and `fftw_plan_many_dft_c2r()` for three planar + RGB transforms with unit stride and per-plane distance; +- one shared kernel spectrum multiplied into each of the three channel + spectra. + +Out-of-place planar buffers are intentionally preferred over an in-place or +strided-interleaved first version. The target host has ample memory, while the +simpler layout reduces alignment, padding, and plan-many mistakes. + +### 6.2 Checked allocation sizes + +Before allocation, check every product used for: + +```text +3 * Pheight * Pwidth * sizeof(double) +3 * Pheight * (Pwidth/2 + 1) * sizeof(fftw_complex) +Pheight * (Pwidth/2 + 1) * sizeof(fftw_complex) +``` + +Allocation failure is a hard fast-mode initialization error with a diagnostic +that includes requested dimensions and byte counts. These numbers are logged +for provenance, not used as a policy cap. + +### 6.3 Packing + +For each frame: + +1. zero the full padded planar real scratch; +2. copy the active `Hss x Wss` region from interleaved RGB into three planes; +3. leave all right and bottom padding zero; +4. keep the original impulse buffer unchanged. + +The zero and pack loops may use one ordinary OpenMP parallel-for region. They +must finish before FFTW execution starts. + +## 7. Kernel transform + +Build the padded real kernel once after plans and scratch buffers exist: + +1. clear a single padded real plane; +2. for `dy = -R..R`, compute the same circular `row_span` as the spatial + reference; +3. copy only `dx = -row_span..row_span` to index `(dy+R, dx+R)`; +4. leave square corners and all FFT padding zero; +5. execute one R2C kernel transform; +6. store its spectrum in immutable `kernel_frequency`; +7. compute the reported retained-flux sum from this same circular mask. + +Do not assume the kernel spectrum is purely real. The unshifted full- +convolution layout places the kernel center at `(R,R)`, so the spectrum +generally has a phase. Complex multiplication must use both real and imaginary +components: + +```text +(ar + i ai) * (br + i bi) + = (ar*br - ai*bi) + i(ar*bi + ai*br). +``` + +The one-time kernel transform and plan creation must not mutate the production +impulse buffer. + +## 8. Per-frame FFTW resolve + +The production resolver executes these stages in order: + +1. **Zero and pack** interleaved impulses into padded planar real scratch. +2. **Forward FFT** all three planes with the cached R2C plan. +3. **Frequency multiply** each RGB spectrum by the shared kernel spectrum. +4. **Inverse FFT** all three planes with the cached C2R plan. +5. **Crop, normalize, and box-downsample** into the caller HDR. + +FFTW's transforms are unnormalized. Let + +```text +fft_scale = 1 / (Pwidth * Pheight) +box_scale = 1 / (N * N). +``` + +For output pixel `(row,column)` and fine offsets `(j,i)`, read inverse scratch +at + +```text +y = R + N*row + j +x = R + N*column + i. +``` + +Then accumulate + +```text +hdr += inverse[y][x] * fft_scale * box_scale. +``` + +Fuse crop, inverse normalization, box summation, RGB interleaving, and HDR +addition in one OpenMP parallel-for over final output rows. Do not materialize +a separate cropped fine-grid image. + +The frequency multiply is embarrassingly parallel over frequency bins and may +use OpenMP if measurement shows it is not already hidden by FFT cost. It must +not run inside an active FFTW OpenMP region. + +## 9. FFTW planning, threads, and wisdom + +### 9.1 Process-level threading lifecycle + +The local host provides `libfftw3_omp`. Use the FFTW threading API exported by +that OpenMP library: + +```c +fftw_init_threads(); +fftw_plan_with_nthreads(fft_workers); +``` + +Requirements: + +- initialize FFTW threading before plan creation; +- create and destroy plans only from the serial control thread; +- never create a plan inside the catalog OpenMP region; +- do not call `fftw_cleanup_threads()` while any accumulator or plan exists; +- keep process-global FFTW runtime ownership/reference counting in one private + module rather than letting individual accumulators independently clean up + global state; +- choose `fft_workers` from the same OpenMP worker policy already used by the + renderer and report it; +- do not set a nested OpenMP region around `fftw_execute()`; +- after FFTW returns, use a separate OpenMP region for crop/downsample. + +`OMP_DYNAMIC` and the user's OpenMP environment remain authoritative. Do not +silently force a different global thread count. + +### 9.2 Planning flags + +Implement in two steps: + +1. correctness and unit tests use `FFTW_ESTIMATE` so tests are fast and do not + depend on host wisdom; +2. the release benchmark compares `FFTW_ESTIMATE` with `FFTW_MEASURE`, recording + plan time separately from steady-state execution. + +Select the release default only from that measured comparison. If +`FFTW_MEASURE` materially improves repeated movie-frame execution, add optional +wisdom import/export in a follow-up or in the same implementation phase: + +- wisdom is host/FFTW-version/dimension specific; +- a missing or incompatible wisdom file must be explicit in verbose output; +- wisdom files are local artifacts and are not committed as portable project + data; +- planning may destroy scratch contents, so plans are created before the + kernel and frame inputs are populated. + +Do not mix plan time into the per-frame convolution metric. + +## 10. Build-system integration + +### 10.1 Dependency detection + +For `PSF_BACKEND=cpu`, require: + +```sh +pkg-config --exists fftw3_omp +pkg-config --modversion fftw3_omp +pkg-config --cflags --libs fftw3_omp +``` + +The current host returns version `3.3.10` and libraries +`-lfftw3_omp -lfftw3`. + +On Gentoo the package requirement is `sci-libs/fftw[openmp]`; do not require +`sci-libs/fftw[threads]` and do not link `libfftw3_threads`. The similarly +named FFTW threading API is implemented by both backends, but exactly one +backend library should provide it in this executable. + +Add CPU-only build variables such as: + +```make +FFTW_CPPFLAGS := $(shell pkg-config --cflags fftw3_omp) +FFTW_LDLIBS := $(shell pkg-config --libs fftw3_omp) +``` + +and a CPU-only definition such as `-DFAST_PSF_FFTW`. Fail early with a clear +Make error if `fftw3_omp` is unavailable in a CPU build. + +Do not append FFTW flags to HIP or dummy backend links. CPU tests that compile +`src/optics.c` and the new FFTW module must receive the same FFTW compile/link +flags. + +The current Makefile obtains common sources through a `src/*.c` wildcard. Do +not let that wildcard pull the FFTW module into HIP and dummy binaries. Exclude +`src/fast_psf_fftw.c` from `COMMON_SOURCES` and append it only to the CPU +renderer and CPU-test source lists. This makes the dependency boundary visible. +Guard calls from `src/optics.c` consistently so HIP and dummy links have no +unresolved FFTW-module symbols. + +Audit at least: + +- renderer backend link; +- `test_frame` and the other CPU tests using `COMMON_SOURCES`; +- the new FFTW unit test; +- `make_psf_fixture` if it links the FFTW-enabled optics object; +- capture/replay tools, ensuring HIP-only tools do not accidentally gain a CPU + FFTW requirement. + +Keep dependency tracking (`-MMD -MP`) for the new module. + +### 10.2 Build identities + +The CPU build always includes FFTW after this change, so a separate executable +suffix is not required. If a temporary compile-time A/B switch is introduced, +include it in the object-directory tag to prevent stale-object reuse. Prefer a +single CPU binary with FFTW production resolution and a test-only spatial +reference over persistent compile-time variants. + +## 11. Diagnostics and timing + +Extend fast-mode reporting with one-time state: + +```text +Fast FFTW: linear min=x, fft=x, + workers=, plan=, plan=, + kernel_fft=, scratch= +``` + +Under `--verbose`, report per-frame stages: + +```text +Fast FFTW frame: zero_pack=, forward=, multiply=, + inverse=, crop_downsample=, total= +``` + +Timing rules: + +- use wall time, not sums of worker-local time; +- plan/kernel setup is one-time and separate; +- `total` spans zero/pack through completed HDR addition; +- keep catalog prefetch/deposit timing outside the FFTW resolve metric; +- include selected dimensions, supersample factor, FWHM/beta, kernel radius, + FFTW version, planning flag, worker environment, and binary revision in + benchmark records. + +These timings validate and attribute the implementation; they are not another +feasibility gate for fast mode. + +## 12. Correctness test plan + +### 12.1 Dedicated FFTW-versus-spatial unit test + +Add `tests/test_fast_psf_fftw.c`. It constructs one accumulator/impulse buffer, +runs both resolvers without re-depositing, and compares their final HDR arrays. +Cover: + +- non-power-of-two image sizes; +- `N = 1, 2, 3, 4`; +- nearest and bilinear deposit modes; +- one white center impulse; +- distinct R/G/B values to catch channel-layout mistakes; +- impulses near all four edges and four corners; +- an impulse exactly on a supersampled cell center; +- phases immediately on both sides of a nearest-cell boundary; +- multiple separated impulses; +- many impulses in one supersampled cell; +- a deterministic dense pseudo-random impulse field; +- an empty impulse buffer; +- a nonzero prefilled HDR background; +- a kernel radius that makes the selected FFT padding larger in both axes; +- dimensions where only one axis needs the next smooth FFT size. + +For each case compute: + +- maximum absolute RGB error; +- maximum relative error above a reference-magnitude floor; +- RMS error; +- total per-channel flux difference; +- worst sample coordinate and channel; +- NaN/Inf count. + +Use provisional double-precision acceptance limits of: + +```text +max_abs <= 1e-10 * max(1, reference_peak) +max_rel <= 1e-9 for |reference| above 1e-12 * reference_peak +relative total-flux error <= 1e-10 +``` + +Tighten these after observing stable results on the target host; do not loosen +them merely to hide crop, wrap, channel, or normalization errors. Any coherent +edge band, shifted peak, factor-of-area scaling, or opposite-edge ghost is an +algorithmic failure regardless of aggregate tolerance. + +### 12.2 Analytic impulse checks + +Before relying only on spatial comparison, add exact index tests with a single +unit impulse: + +- verify the full-convolution peak appears at impulse index plus `R`; +- verify the cropped fine-grid center returns to the original impulse index; +- verify no signal appears at the opposite edge; +- verify a unit impulse's cropped result equals the circular kernel samples; +- verify box-downsample indices for every phase `0..N-1` in both axes; +- verify RGB channels do not cross-contaminate. + +These tests diagnose layout errors that a broad image statistic can obscure. + +### 12.3 Existing regression suite + +After FFTW becomes the production resolver: + +```sh +make -j4 BUILD_TYPE=Debug PSF_BACKEND=cpu test +``` + +must pass. In particular: + +- the existing fast-mode snapped-kernel regression now exercises FFTW; +- the additive-HDR regression must still preserve its `0.25` background; +- ordinary non-fast FITS references remain unchanged; +- Minkowski and Schwarzschild camera/lens-map regressions remain unchanged. + +Add CLI smoke coverage for `--fast-mode` in a CPU build. HIP and dummy builds +must continue to reject it with the existing message. + +### 12.4 Scientific/visual metrics + +Using identical impulses, compare FFTW and spatial fast-mode outputs for: + +- centroid; +- fitted FWHM; +- fitted Moffat beta relative to the same direct-reference fitting procedure; +- total linear-HDR flux; +- center, edge, and corner crops. + +The FFTW change should add only floating-point roundoff to the already measured +nearest/bilinear approximation. It must not change the deposition-error table +or the rationale for nearest as default. + +## 13. Performance benchmark plan + +### 13.1 Resolver-only benchmark + +Add `tests/benchmark_fast_psf_fftw.c` or an equivalent bounded benchmark target +that: + +- initializes one accumulator; +- deposits or directly fills one deterministic impulse buffer once; +- runs spatial and FFTW resolution on the identical buffer; +- separates first-use plan/setup from repeated execution; +- performs warm-up before measured FFT executions; +- records at least five steady-state executions when runtime permits; +- reports median, minimum, and all raw samples; +- hashes the impulse buffer and both HDR outputs; +- records RSS and FFTW dimensions for provenance. + +Required cases: + +1. current default PSF, N=2, 320x240; +2. current default PSF, N=4, 320x240; +3. a medium frame that is large enough for FFT execution to dominate startup; +4. the intended 4K N=2 target; +5. 4K N=4 if the user requests that mode to be production-supported. + +The 4K cases are bounded synthetic/identical-buffer resolver tests, not another +catalog or geodesic run. + +### 13.2 End-to-end comparison + +Repeat the recorded 320x240 Galactic-center command from +`benchmarks/fast_mode_cpu_2026-09-18.md` with the same catalog, geometry, +exposure, supersample factor, deposit mode, worker environment, and revision +metadata. Preserve raw terminal output rather than only a timing table. + +Do not claim an end-to-end speedup unless inputs, event/image counts, PSF +parameters, and output mode match. Do not add worker-summed times to process +wall time. + +An expensive full-sky or production 4K catalog render is not automatically run +as part of implementation. It requires separate explicit authorization and +must use a frozen lens map/input provenance when used for final evidence. + +### 13.3 Performance acceptance + +The FFTW implementation is ready to become the default resolver when: + +- all correctness tests pass; +- the resolver-only benchmark shows a clear steady-state reduction for default + N=2 and the N=4 case that exposed spatial scaling; +- plan/setup time is reported separately and is acceptable for the intended + single-frame/movie usage, or wisdom support addresses it; +- the recorded end-to-end bounded command does not regress catalog/deposit + behavior or image count; +- FFTW output differences remain within the linear-HDR tolerances and do not + change measured centroid/FWHM/beta beyond roundoff. + +No new arbitrary global wall-time target is introduced by this plan. + +## 14. Failure handling and cleanup + +- Reject invalid/overflowed FFT dimensions before allocation or plan creation. +- Treat `fftw_alloc_* == NULL` or a null plan as initialization failure. +- Print which allocation/plan failed, dimensions, plan mode, and worker count. +- Never continue with a partially initialized kernel spectrum. +- Never silently switch to spatial convolution after an FFTW failure. +- Destroy forward/inverse/kernel plans before freeing their buffers. +- Release every FFTW allocation in the partial-initialization error path. +- Ensure movie failure exits and imported-lens-map early returns destroy the + FFTW state through the existing accumulator cleanup. +- Do not call global FFTW cleanup while another accumulator remains alive. +- Keep `fast_psf_accumulator_destroy()` idempotent on a zero state. + +## 15. Implementation sequence + +### Phase 1: Freeze the spatial reference + +- Extract the existing resolver without changing its arithmetic. +- Cache/reuse its circular row spans. +- Add a small digest/reference test for the current spatial output. +- Add checked smooth-size and product helpers with unit tests. + +Exit criterion: current fast-mode tests and the spatial output digest pass. + +### Phase 2: Add FFTW build and private state + +- Add CPU-only `fftw3_omp` detection and link flags. +- Add the private FFTW module and opaque ownership. +- Allocate planar real/frequency scratch. +- Initialize process-level FFTW threading. +- Create `FFTW_ESTIMATE` plans and build the circular kernel spectrum. +- Report plan dimensions and setup time. + +Exit criterion: all CPU/HIP/dummy build variants still compile with the correct +dependency boundary; state construction/destruction passes allocation tests. + +### Phase 3: Implement exact FFT resolution + +- Implement zero/pack, batched forward FFT, shared-kernel multiplication, + batched inverse FFT, `(R,R)` crop, normalization, box downsample, and additive + HDR write. +- Add analytic impulse/crop tests first. +- Add dense FFTW-versus-spatial comparisons. + +Exit criterion: dedicated correctness tests pass for all required dimensions, +phases, channels, and edge cases. + +### Phase 4: Integrate the production path + +- Make `fast_psf_accumulator_resolve()` use FFTW in CPU fast mode. +- Retain the spatial helper only for tests/benchmarks. +- Add CLI smoke coverage and verbose timing. +- Run the full CPU Debug suite. + +Exit criterion: existing behavior and references pass, and additive HDR +semantics remain covered. + +### Phase 5: Benchmark planning and execution + +- Run the resolver-only matrix. +- Compare `FFTW_ESTIMATE` and `FFTW_MEASURE`. +- Decide whether wisdom support is required. +- Repeat the bounded Galactic-center command with full provenance. +- Save raw logs and output hashes. + +Exit criterion: correctness remains within tolerance and FFTW is demonstrably +faster than the spatial resolver for the default and N=4 cases. + +### Phase 6: Documentation finalization + +- Update `usage.md` with FFTW dependency and planning/wisdom behavior. +- Update `nr_spacetime_movie_renderer_design.md` to state that the global fast + convolution is implemented as zero-padded FFTW linear convolution. +- Add a dated benchmark record with commands, revision, environment, raw logs, + hashes, plan/setup time, and steady-state execution samples. +- Keep the existing deposition-error record unchanged except for links to the + FFTW validation record; FFTW does not change nearest/bilinear errors. + +## 16. Review checklist + +Before considering the implementation complete, review the following +explicitly: + +- [ ] FFTW is linked only for the CPU PSF backend and CPU tests. +- [ ] Local dependency is `fftw3_omp`, version and libraries are reported. +- [ ] Padding satisfies `Wss + 2R` and `Hss + 2R` before smooth rounding. +- [ ] Kernel square corners outside the circular support are zero. +- [ ] Kernel is stored unshifted at indices `(dy+R, dx+R)`. +- [ ] Fine-grid crop begins at `(R,R)`. +- [ ] Inverse normalization uses `1/(Pwidth*Pheight)` exactly once. +- [ ] Box normalization uses `1/N^2` exactly once. +- [ ] Output uses `+=` into HDR. +- [ ] RGB batch layout and distances are covered by distinct-channel tests. +- [ ] Empty input produces no HDR change. +- [ ] No opposite-edge wraparound appears. +- [ ] Original impulse buffer is unchanged by resolve. +- [ ] Plans and kernel spectrum are reused across frames. +- [ ] Planning occurs outside OpenMP producer regions. +- [ ] FFTW and renderer OpenMP stages do not nest/oversubscribe. +- [ ] Partial initialization and all early exits free FFTW state. +- [ ] Spatial reference and FFTW consume the identical impulse buffer. +- [ ] Linear-HDR errors, flux, centroid, FWHM, and beta are recorded. +- [ ] Plan time and steady-state execution time are reported separately. +- [ ] Benchmark commands, revisions, inputs, event counts, and raw output are + preserved. +- [ ] No expensive full-sky/4K catalog render is run without explicit + authorization. + +## 17. Expected final architecture + +```text +catalog / inverse lens map / g / blackbody / flux + | + v + atomic nearest or bilinear SS delta deposit + | + v + interleaved double-RGB impulse buffer + | + v + zero + pack to FFTW planar padded real arrays + | + v + batched RGB R2C forward FFT + | + v + multiply by one cached Moffat kernel spectrum + | + v + batched RGB C2R inverse FFT + | + v + crop at (R,R) + normalize + N x N downsample + | + v + additive double-RGB HDR framebuffer +``` + +The physical/catalog decisions remain on the CPU producer side; FFTW replaces +only the global convolution implementation inside the explicit preview path. diff --git a/nr_spacetime_movie_renderer_design.md b/nr_spacetime_movie_renderer_design.md index c9ddb2f..46e76fd 100644 --- a/nr_spacetime_movie_renderer_design.md +++ b/nr_spacetime_movie_renderer_design.md @@ -851,7 +851,13 @@ pixel-area 积分，所以核本身不改变指定的 FWHM/beta。`nearest` 在 沉积误差测量、`N` 与取舍见 [`benchmarks/fast_mode_deposit_2026-09-18.md`](benchmarks/fast_mode_deposit_2026-09-18.md)。 该模式只支持 CPU PSF 后端；`--max-cache-psf-flux` 不适用（所有事件共用同一 -个全局核半径）。 +个全局核半径）。CPU 构建中该全局卷积实现为零填充的 double 精度 FFTW 线性 +卷积（`fftw3`/`fftw3_omp`）：核频谱、FFT plans 与 planar scratch 在初始化时 +一次性建立并跨帧复用，`resolve()` 只做 zero/pack、批量 R2C、频域乘、批量 +C2R，以及在 `(R,R)` 处的裁剪、`1/(Pwidth·Pheight)` 与 `1/N²` 归一化和 +`+=` 到 HDR。旧的嵌套空间循环作为测试/基准参考保留，不是运行时 fallback。 +依赖与 plan 模式取舍见 `build.md` 与 +[`benchmarks/fast_mode_fftw_2026-09-25.md`](benchmarks/fast_mode_fftw_2026-09-25.md)。 PSF 第一版可用 Gaussian； 以后可换成 Airy 或其他相机模型。 diff --git a/src/fast_psf_fftw.c b/src/fast_psf_fftw.c new file mode 100644 index 0000000..df6f592 --- /dev/null +++ b/src/fast_psf_fftw.c @@ -0,0 +1,393 @@ +#include "fast_psf_fftw.h" + +#include +#include +#include +#include +#include +#include +#include + +/* Process-global FFTW threading lifecycle. Plans are created only from the + * serial control thread, and global cleanup happens only after the last live + * state is gone. This keeps the individual states from racing on FFTW's + * process-wide planner state. */ +static int g_threading_initialized = 0; +static int g_live_states = 0; +static int g_plan_measure = 0; + +void fast_psf_fftw_set_plan_mode(int measure) +{ + g_plan_measure = measure ? 1 : 0; +} + +static unsigned plan_flags(void) +{ + return g_plan_measure ? FFTW_MEASURE : FFTW_ESTIMATE; +} + +struct FastPsfFftwState { + int fft_width, fft_height; + int ss_width, ss_height; + int final_width, final_height; + int supersample; + int kernel_radius; + int workers; + int registered; + int plan_measure; + size_t scratch_bytes; + double setup_seconds; + double kernel_seconds; + double fft_scale; + double box_scale; + + double *real_rgb; /* 3 planes of fft_height * fft_width */ + fftw_complex *freq_rgb; /* 3 planes of fft_height * (fft_width/2 + 1) */ + fftw_complex *kernel_freq; /* one plane of fft_height * (fft_width/2 + 1) */ + fftw_plan forward_plan; + fftw_plan inverse_plan; +}; + +static int checked_mul_size(size_t a, size_t b, size_t *out) +{ + if (a != 0 && b > SIZE_MAX / a) + return -1; + *out = a * b; + return 0; +} + +int fast_psf_fftw_next_smooth_size(size_t min_extent, size_t *out) +{ + if (out == NULL || min_extent == 0) + return -1; + if (min_extent > (size_t)INT_MAX) + return -1; + static const unsigned factors[] = {2, 3, 5, 7}; + for (size_t candidate = min_extent; candidate <= (size_t)INT_MAX; + ++candidate) { + size_t remaining = candidate; + for (size_t i = 0; i < sizeof factors / sizeof factors[0]; ++i) + while (remaining % factors[i] == 0) + remaining /= factors[i]; + if (remaining == 1) { + *out = candidate; + return 0; + } + } + return -1; +} + +void fast_psf_fftw_destroy(FastPsfFftwState *state) +{ + if (state == NULL) + return; + if (state->forward_plan != NULL) + fftw_destroy_plan(state->forward_plan); + if (state->inverse_plan != NULL) + fftw_destroy_plan(state->inverse_plan); + if (state->real_rgb != NULL) + fftw_free(state->real_rgb); + if (state->freq_rgb != NULL) + fftw_free(state->freq_rgb); + if (state->kernel_freq != NULL) + fftw_free(state->kernel_freq); + const int registered = state->registered; + free(state); + + if (!registered) + return; + g_live_states--; + if (g_live_states <= 0) { + g_live_states = 0; + if (g_threading_initialized) { + fftw_cleanup_threads(); + g_threading_initialized = 0; + } + } +} + +FastPsfFftwState *fast_psf_fftw_create(int ss_width, int ss_height, + int final_width, int final_height, + int supersample, int kernel_radius, + const float *weights, + const int *row_span, int fft_workers, + double *setup_seconds, + double *kernel_seconds) +{ + if (setup_seconds != NULL) + *setup_seconds = 0.0; + if (kernel_seconds != NULL) + *kernel_seconds = 0.0; + if (ss_width <= 0 || ss_height <= 0 || final_width <= 0 || + final_height <= 0 || supersample < 1 || kernel_radius < 0 || + weights == NULL || row_span == NULL || fft_workers < 1) { + fprintf(stderr, + "Fast FFTW init failed (arguments): ss=%dx%d final=%dx%d " + "supersample=%d kernel_radius=%d workers=%d\n", + ss_width, ss_height, final_width, final_height, supersample, + kernel_radius, fft_workers); + return NULL; + } + if ((size_t)final_width * (size_t)supersample != (size_t)ss_width || + (size_t)final_height * (size_t)supersample != (size_t)ss_height) { + fprintf(stderr, + "Fast FFTW init failed (dimensions): ss=%dx%d is not " + "final=%dx%d times supersample=%d\n", + ss_width, ss_height, final_width, final_height, supersample); + return NULL; + } + + const char *stage = "state allocation"; + FastPsfFftwState *state = calloc(1, sizeof *state); + if (state == NULL) { + fprintf(stderr, "Fast FFTW init failed (state allocation): %zu bytes\n", + sizeof *state); + return NULL; + } + state->ss_width = ss_width; + state->ss_height = ss_height; + state->final_width = final_width; + state->final_height = final_height; + state->supersample = supersample; + state->kernel_radius = kernel_radius; + state->workers = fft_workers; + state->plan_measure = g_plan_measure; + + size_t real_bytes = 0, freq_bytes = 0, kernel_bytes = 0; + size_t fft_width = 0, fft_height = 0; + const size_t radius = (size_t)kernel_radius; + size_t width_min, height_min; + stage = "kernel extent"; + if (checked_mul_size(radius, 2, &width_min) || + checked_mul_size(radius, 2, &height_min)) + goto fail; + if ((size_t)ss_width > SIZE_MAX - width_min || + (size_t)ss_height > SIZE_MAX - height_min) + goto fail; + width_min += (size_t)ss_width; + height_min += (size_t)ss_height; + stage = "smooth FFT size"; + if (fast_psf_fftw_next_smooth_size(width_min, &fft_width) || + fast_psf_fftw_next_smooth_size(height_min, &fft_height)) + goto fail; + state->fft_width = (int)fft_width; + state->fft_height = (int)fft_height; + + const size_t complex_width = fft_width / 2 + 1; + size_t real_plane, real_count; + size_t freq_plane, freq_count; + stage = "allocation size"; + if (checked_mul_size(fft_height, fft_width, &real_plane) || + checked_mul_size(real_plane, 3, &real_count) || + checked_mul_size(real_count, sizeof(double), &real_bytes) || + checked_mul_size(fft_height, complex_width, &freq_plane) || + checked_mul_size(freq_plane, 3, &freq_count) || + checked_mul_size(freq_count, sizeof(fftw_complex), &freq_bytes) || + checked_mul_size(freq_plane, sizeof(fftw_complex), &kernel_bytes) || + real_bytes > SIZE_MAX - freq_bytes - kernel_bytes) + goto fail; + /* FFTW's plan_many interface takes idist/odist as int, so the per-plane + * element counts must fit before the cast below even when each axis is a + * valid int. */ + stage = "plan_many stride"; + if (real_plane > (size_t)INT_MAX || freq_plane > (size_t)INT_MAX) + goto fail; + state->scratch_bytes = real_bytes + freq_bytes + kernel_bytes; + + stage = "FFTW threading initialization"; + if (g_live_states == 0) { + if (!fftw_init_threads()) + goto fail; + g_threading_initialized = 1; + } + g_live_states++; + state->registered = 1; + fftw_plan_with_nthreads(fft_workers); + + const double plan_start = omp_get_wtime(); + stage = "FFTW scratch allocation"; + state->real_rgb = fftw_alloc_real(real_count); + state->freq_rgb = fftw_alloc_complex(freq_count); + state->kernel_freq = fftw_alloc_complex(freq_plane); + if (state->real_rgb == NULL || state->freq_rgb == NULL || + state->kernel_freq == NULL) + goto fail; + + int n[2] = {state->fft_height, state->fft_width}; + int real_embed[2] = {state->fft_height, state->fft_width}; + int complex_embed[2] = {state->fft_height, (int)complex_width}; + const int real_dist = (int)real_plane; + const int complex_dist = (int)freq_plane; + stage = "FFTW plan creation"; + state->forward_plan = fftw_plan_many_dft_r2c( + 2, n, 3, state->real_rgb, real_embed, 1, real_dist, state->freq_rgb, + complex_embed, 1, complex_dist, plan_flags()); + state->inverse_plan = fftw_plan_many_dft_c2r( + 2, n, 3, state->freq_rgb, complex_embed, 1, complex_dist, + state->real_rgb, real_embed, 1, real_dist, plan_flags()); + if (state->forward_plan == NULL || state->inverse_plan == NULL) + goto fail; + state->setup_seconds = omp_get_wtime() - plan_start; + if (setup_seconds != NULL) + *setup_seconds = state->setup_seconds; + + /* One-time kernel transform. The kernel is stored unshifted at padded origin + * indices (dy + R, dx + R), retaining only the circular mask the spatial + * reference traverses; square corners and FFT padding stay zero. The plan is + * created before the kernel is written because FFTW_MEASURE overwrites its + * input during planning. */ + stage = "kernel scratch allocation"; + double *kernel_real = fftw_alloc_real(real_plane); + fftw_plan kernel_plan = NULL; + if (kernel_real == NULL) + goto fail; + const double kernel_start = omp_get_wtime(); + stage = "kernel plan creation"; + kernel_plan = fftw_plan_dft_r2c_2d(state->fft_height, state->fft_width, + kernel_real, state->kernel_freq, + plan_flags()); + if (kernel_plan == NULL) { + fftw_free(kernel_real); + goto fail; + } + memset(kernel_real, 0, real_plane * sizeof *kernel_real); + const size_t side = 2 * radius + 1; + for (int dy = -kernel_radius; dy <= kernel_radius; ++dy) { + const int span = row_span[dy + kernel_radius]; + for (int dx = -span; dx <= span; ++dx) + kernel_real[(size_t)(dy + kernel_radius) * fft_width + + (size_t)(dx + kernel_radius)] = + (double)weights[(size_t)(dy + kernel_radius) * side + + (size_t)(dx + kernel_radius)]; + } + fftw_execute(kernel_plan); + state->kernel_seconds = omp_get_wtime() - kernel_start; + if (kernel_seconds != NULL) + *kernel_seconds = state->kernel_seconds; + fftw_destroy_plan(kernel_plan); + fftw_free(kernel_real); + + state->fft_scale = 1.0 / ((double)fft_width * (double)fft_height); + state->box_scale = 1.0 / ((double)supersample * (double)supersample); + return state; + +fail: + fprintf(stderr, + "Fast FFTW init failed at '%s': ss=%dx%d, linear min=%zux%zu, " + "fft=%zux%zu, plan=%s, workers=%d, bytes real=%zu freq=%zu " + "kernel=%zu total=%zu\n", + stage, ss_width, ss_height, + (size_t)ss_width + 2 * radius, + (size_t)ss_height + 2 * radius, fft_width, fft_height, + state->plan_measure ? "measure" : "estimate", fft_workers, + real_bytes, freq_bytes, kernel_bytes, + real_bytes + freq_bytes + kernel_bytes); + fast_psf_fftw_destroy(state); + return NULL; +} + +int fast_psf_fftw_resolve(FastPsfFftwState *state, + const double *interleaved_ss_rgb, + double *interleaved_hdr_rgb, + FastPsfFftwFrameTiming *timing) +{ + if (state == NULL || interleaved_ss_rgb == NULL || + interleaved_hdr_rgb == NULL) + return -1; + FastPsfFftwFrameTiming local = {0}; + const double total_start = omp_get_wtime(); + + const size_t plane = (size_t)state->fft_height * state->fft_width; + const size_t complex_width = (size_t)state->fft_width / 2 + 1; + const size_t complex_plane = (size_t)state->fft_height * complex_width; + const size_t real_count = 3 * plane; + const size_t freq_count = 3 * complex_plane; + + double start = omp_get_wtime(); +#pragma omp parallel for schedule(static) num_threads(state->workers) \ + if (state->workers > 1) + for (size_t i = 0; i < real_count; ++i) + state->real_rgb[i] = 0.0; +#pragma omp parallel for schedule(static) num_threads(state->workers) \ + if (state->workers > 1) + for (int y = 0; y < state->ss_height; ++y) { + for (int x = 0; x < state->ss_width; ++x) { + const size_t source = + 3 * ((size_t)y * state->ss_width + (size_t)x); + const size_t target = (size_t)y * state->fft_width + (size_t)x; + state->real_rgb[target] = interleaved_ss_rgb[source]; + state->real_rgb[plane + target] = interleaved_ss_rgb[source + 1]; + state->real_rgb[2 * plane + target] = interleaved_ss_rgb[source + 2]; + } + } + local.zero_pack_seconds = omp_get_wtime() - start; + + start = omp_get_wtime(); + fftw_execute(state->forward_plan); + local.forward_seconds = omp_get_wtime() - start; + + start = omp_get_wtime(); +#pragma omp parallel for schedule(static) num_threads(state->workers) \ + if (state->workers > 1) + for (size_t i = 0; i < freq_count; ++i) { + const double ar = state->freq_rgb[i][0]; + const double ai = state->freq_rgb[i][1]; + const double br = state->kernel_freq[i % complex_plane][0]; + const double bi = state->kernel_freq[i % complex_plane][1]; + state->freq_rgb[i][0] = ar * br - ai * bi; + state->freq_rgb[i][1] = ar * bi + ai * br; + } + local.multiply_seconds = omp_get_wtime() - start; + + start = omp_get_wtime(); + fftw_execute(state->inverse_plan); + local.inverse_seconds = omp_get_wtime() - start; + + start = omp_get_wtime(); + const int radius = state->kernel_radius; + const int supersample = state->supersample; + const double scale = state->fft_scale * state->box_scale; +#pragma omp parallel for schedule(static) num_threads(state->workers) \ + if (state->workers > 1) + for (int row = 0; row < state->final_height; ++row) { + for (int column = 0; column < state->final_width; ++column) { + double sum[3] = {0.0, 0.0, 0.0}; + const int base_y = radius + supersample * row; + const int base_x = radius + supersample * column; + for (int j = 0; j < supersample; ++j) + for (int i = 0; i < supersample; ++i) { + const size_t offset = + (size_t)(base_y + j) * state->fft_width + + (size_t)(base_x + i); + sum[0] += state->real_rgb[offset]; + sum[1] += state->real_rgb[plane + offset]; + sum[2] += state->real_rgb[2 * plane + offset]; + } + const size_t out = + 3 * ((size_t)row * state->final_width + (size_t)column); + interleaved_hdr_rgb[out] += sum[0] * scale; + interleaved_hdr_rgb[out + 1] += sum[1] * scale; + interleaved_hdr_rgb[out + 2] += sum[2] * scale; + } + } + local.crop_downsample_seconds = omp_get_wtime() - start; + local.total_seconds = omp_get_wtime() - total_start; + if (timing != NULL) + *timing = local; + return 0; +} + +void fast_psf_fftw_report(const FastPsfFftwState *state, FILE *stream) +{ + if (state == NULL || stream == NULL) + return; + fprintf(stream, + "Fast FFTW: linear min=%dx%d, fft=%dx%d, workers=%d, " + "plan=%s, plan=%.6f s, kernel_fft=%.6f s, scratch=%.1f MiB\n", + state->ss_height + 2 * state->kernel_radius, + state->ss_width + 2 * state->kernel_radius, state->fft_height, + state->fft_width, state->workers, + state->plan_measure ? "measure" : "estimate", state->setup_seconds, + state->kernel_seconds, + (double)state->scratch_bytes / (1024.0 * 1024.0)); +} diff --git a/src/fast_psf_fftw.h b/src/fast_psf_fftw.h new file mode 100644 index 0000000..7e6c818 --- /dev/null +++ b/src/fast_psf_fftw.h @@ -0,0 +1,67 @@ +#ifndef FAST_PSF_FFTW_H +#define FAST_PSF_FFTW_H + +#include +#include + +/* Private FFTW-backed global convolution for fast-mode PSF accumulation. + * + * This module owns every FFTW allocation and plan. It replaces only the + * spatial convolution inside the explicit fast-mode preview path; the impulse + * buffer stays owned by the caller, and the HDR framebuffer is caller-owned and + * updated with +=. It is compiled only for the CPU PSF backend (FAST_PSF_FFTW) + * so the HIP and dummy backends do not link FFTW. */ + +typedef struct FastPsfFftwState FastPsfFftwState; + +typedef struct { + double zero_pack_seconds; + double forward_seconds; + double multiply_seconds; + double inverse_seconds; + double crop_downsample_seconds; + double total_seconds; +} FastPsfFftwFrameTiming; + +/* Builds one reusable convolution state from the immutable kernel `weights` + * (side 2*kernel_radius+1, row-major) and the cached circular `row_span` + * (length 2*kernel_radius+1, the number of retained dx at each dy). The state + * computes a zero-padded linear convolution over Wss x Hss and downsamples by + * `supersample`. `setup_seconds` receives plan creation time and + * `kernel_seconds` receives the one-time kernel transform time. Returns NULL + * on any allocation or plan failure; nothing is leaked. */ +FastPsfFftwState *fast_psf_fftw_create(int ss_width, int ss_height, + int final_width, int final_height, + int supersample, int kernel_radius, + const float *weights, + const int *row_span, int fft_workers, + double *setup_seconds, + double *kernel_seconds); + +/* Convolves `interleaved_ss_rgb` (3 doubles per supersampled pixel, row-major) + * into `interleaved_hdr_rgb` using +=. The input buffer is not modified. + * Returns 0 on success and -1 on invalid state. */ +int fast_psf_fftw_resolve(FastPsfFftwState *state, + const double *interleaved_ss_rgb, + double *interleaved_hdr_rgb, + FastPsfFftwFrameTiming *timing); + +/* Safe for NULL and partially initialized states; releases plans before the + * buffers they reference. Does not touch process-global FFTW state except on + * the final live state. */ +void fast_psf_fftw_destroy(FastPsfFftwState *state); + +/* Prints one diagnostic line when `state` is non-NULL. */ +void fast_psf_fftw_report(const FastPsfFftwState *state, FILE *stream); + +/* First value >= min_extent whose prime factors are limited to {2,3,5,7}. + * Returns 0 and stores the result in *out on success; returns -1 when the + * result would exceed INT_MAX or on overflow. */ +int fast_psf_fftw_next_smooth_size(size_t min_extent, size_t *out); + +/* Selects FFTW_MEASURE (measure != 0) instead of FFTW_ESTIMATE for plans + * created after this call. Intended for the benchmark's plan-mode comparison; + * production uses the default FFTW_ESTIMATE. */ +void fast_psf_fftw_set_plan_mode(int measure); + +#endif diff --git a/src/main.c b/src/main.c index 8b5be48..3571aa3 100644 --- a/src/main.c +++ b/src/main.c @@ -977,6 +977,7 @@ static int render_lens_map(const Settings *s, StarCatalog *catalog) { if (s->fast_mode) { fast = &local_fast; fast_psf_accumulator_report(&local_fast, stderr); + fast_psf_accumulator_set_verbose(&local_fast, s->verbose); } for (size_t i = 0; i < map.frame_count; ++i) { const char *output_path = s->output_path; @@ -1218,6 +1219,7 @@ int main(int argc, char **argv) { } settings.fast_psf = &fast_accumulator; fast_psf_accumulator_report(&fast_accumulator, stderr); + fast_psf_accumulator_set_verbose(&fast_accumulator, settings.verbose); } else if (settings.fast_mode) { fputs("Fast mode on an imported lens map uses the map's own dimensions.\n", stderr); diff --git a/src/optics.c b/src/optics.c index 89cdb42..a84608f 100644 --- a/src/optics.c +++ b/src/optics.c @@ -1,5 +1,9 @@ #include "optics.h" +#ifdef FAST_PSF_FFTW +#include "fast_psf_fftw.h" +#endif + #include #include #include @@ -500,12 +504,19 @@ int fast_psf_accumulator_init(FastPsfAccumulator *accumulator, int width, side * side > SIZE_MAX / sizeof(float)) return -1; float *weights = malloc(side * side * sizeof *weights); + int *row_span = malloc(side * sizeof *row_span); double *buffer = calloc(ss_width * ss_height * 3, sizeof *buffer); - if (weights == NULL || buffer == NULL) { + if (weights == NULL || row_span == NULL || buffer == NULL) { free(weights); + free(row_span); free(buffer); return -1; } + for (int dy = -radius; dy <= radius; ++dy) { + const double remaining = + (double)radius * radius - (double)dy * dy; + row_span[dy + radius] = remaining > 0.0 ? (int)sqrt(remaining) : 0; + } *accumulator = (FastPsfAccumulator){ .fwhm_pixels = psf->fwhm_pixels, .moffat_beta = psf->moffat_beta, @@ -522,6 +533,7 @@ int fast_psf_accumulator_init(FastPsfAccumulator *accumulator, int width, .radius_pixels = radius, .use_reference = use_reference, .weights = weights, + .row_span = row_span, .buffer = buffer}; /* k[m] is the pixel-area integral of I(z/N; alpha, beta) over ss cell m, * i.e. the final-normalized Moffat evaluated at the supersampled scale. @@ -539,6 +551,23 @@ int fast_psf_accumulator_init(FastPsfAccumulator *accumulator, int width, weights[(size_t)(dy + radius) * side + (dx + radius)] = (float)(weight * normalization_scale); } +#ifdef FAST_PSF_FFTW + int fft_workers = omp_get_max_threads(); + if (fft_workers < 1) + fft_workers = 1; + accumulator->fftw = fast_psf_fftw_create( + (int)ss_width, (int)ss_height, width, height, supersample, radius, + weights, row_span, fft_workers, &accumulator->fftw_setup_seconds, + &accumulator->fftw_kernel_seconds); + if (accumulator->fftw == NULL) { + free(weights); + free(row_span); + free(buffer); + *accumulator = (FastPsfAccumulator){0}; + return -1; + } + accumulator->fftw_enabled = 1; +#endif return 0; } @@ -590,24 +619,20 @@ int fast_psf_accumulator_deposit(FastPsfAccumulator *accumulator, double x, return wing_clipped; } -int fast_psf_accumulator_resolve(const FastPsfAccumulator *accumulator, - double *hdr, int worker_count) +/* Spatial reference global convolution. Kept as a test/benchmark reference + * with the identical arithmetic the FFTW path must reproduce; it is never a + * silent runtime fallback. */ +int fast_psf_accumulator_resolve_spatial_reference( + FastPsfAccumulator *accumulator, double *hdr, int worker_count) { const int supersample = accumulator == NULL ? 0 : accumulator->supersample; const int radius = accumulator == NULL ? 0 : accumulator->radius_pixels; const size_t side = (size_t)2 * radius + 1; if (accumulator == NULL || accumulator->weights == NULL || - accumulator->buffer == NULL || hdr == NULL || supersample <= 0) + accumulator->row_span == NULL || accumulator->buffer == NULL || + hdr == NULL || supersample <= 0) return -1; - int *row_span = malloc(side * sizeof *row_span); - if (row_span == NULL) - return -1; - for (int dy = -radius; dy <= radius; ++dy) { - const double remaining = - (double)radius * radius - (double)dy * dy; - row_span[dy + radius] = - remaining > 0.0 ? (int)sqrt(remaining) : 0; - } + const int *row_span = accumulator->row_span; const double inverse_block = 1.0 / ((double)supersample * supersample); const int threads = worker_count > 0 ? worker_count : 1; #pragma omp parallel for schedule(static) num_threads(threads) if (threads > 1) @@ -652,15 +677,54 @@ int fast_psf_accumulator_resolve(const FastPsfAccumulator *accumulator, hdr[out + 2] += sum[2] * inverse_block; } } - free(row_span); return 0; } +int fast_psf_accumulator_resolve(FastPsfAccumulator *accumulator, + double *hdr, int worker_count) +{ + if (accumulator == NULL || accumulator->buffer == NULL || hdr == NULL) + return -1; +#ifdef FAST_PSF_FFTW + if (accumulator->fftw_enabled && accumulator->fftw != NULL) { + FastPsfFftwFrameTiming timing = {0}; + const double start = omp_get_wtime(); + if (fast_psf_fftw_resolve(accumulator->fftw, accumulator->buffer, hdr, + &timing)) + return -1; + accumulator->fftw_frame_seconds = omp_get_wtime() - start; + if (accumulator->verbose) + fprintf(stderr, + "Fast FFTW frame: zero_pack=%.6f forward=%.6f " + "multiply=%.6f inverse=%.6f crop_downsample=%.6f " + "total=%.6f\n", + timing.zero_pack_seconds, timing.forward_seconds, + timing.multiply_seconds, timing.inverse_seconds, + timing.crop_downsample_seconds, timing.total_seconds); + return 0; + } +#endif + return fast_psf_accumulator_resolve_spatial_reference(accumulator, hdr, + worker_count); +} + +void fast_psf_accumulator_set_verbose(FastPsfAccumulator *accumulator, + int verbose) +{ + if (accumulator != NULL) + accumulator->verbose = verbose ? 1 : 0; +} + void fast_psf_accumulator_destroy(FastPsfAccumulator *accumulator) { if (accumulator == NULL) return; +#ifdef FAST_PSF_FFTW + fast_psf_fftw_destroy(accumulator->fftw); + accumulator->fftw = NULL; +#endif free(accumulator->weights); + free(accumulator->row_span); free(accumulator->buffer); *accumulator = (FastPsfAccumulator){0}; } @@ -687,6 +751,9 @@ void fast_psf_accumulator_report(const FastPsfAccumulator *accumulator, sum / ((double)accumulator->supersample * accumulator->supersample), accumulator->use_reference ? "reference 8-point" : "cached 4-point"); +#ifdef FAST_PSF_FFTW + fast_psf_fftw_report(accumulator->fftw, stream); +#endif } diff --git a/src/optics.h b/src/optics.h index 2306a7d..ffdb25e 100644 --- a/src/optics.h +++ b/src/optics.h @@ -26,6 +26,10 @@ typedef enum { FAST_PSF_DEPOSIT_BILINEAR = 1, } FastPsfDeposit; +/* Private FFTW convolution state. Defined in fast_psf_fftw.c; only the CPU + * PSF backend may reference the implementation. */ +typedef struct FastPsfFftwState FastPsfFftwState; + /* Fast point-source accumulation: every image event is deposited as a delta * (one nearest supersampled pixel, or 4 bilinear pixels) into one shared * supersampled HDR buffer. A single immutable global kernel is convolved once @@ -45,7 +49,16 @@ typedef struct { int radius_pixels; /* kernel radius in supersampled pixels */ int use_reference; /* 8-point instead of 4-point kernel quadrature */ float *weights; /* (2R+1)^2 pixel-area kernel, row-major */ + int *row_span; /* cached circular half-width per dy, length 2R+1 */ double *buffer; /* supersampled HDR, 3 channels per pixel */ +#ifdef FAST_PSF_FFTW + FastPsfFftwState *fftw; /* owned private convolution state */ + int fftw_enabled; + double fftw_setup_seconds; + double fftw_kernel_seconds; + double fftw_frame_seconds; /* most recent resolve */ +#endif + int verbose; } FastPsfAccumulator; typedef struct { @@ -110,8 +123,14 @@ void fast_psf_accumulator_clear(FastPsfAccumulator *accumulator); * kernel radius (wing clipped), and 3 when discarded by --psf-min-y. */ int fast_psf_accumulator_deposit(FastPsfAccumulator *accumulator, double x, double y, LinearRgb color, double flux); -int fast_psf_accumulator_resolve(const FastPsfAccumulator *accumulator, +int fast_psf_accumulator_resolve(FastPsfAccumulator *accumulator, double *hdr, int worker_count); +/* Test/benchmark reference only: the original nested-loop global convolution. + * Production resolve never falls back to it. */ +int fast_psf_accumulator_resolve_spatial_reference( + FastPsfAccumulator *accumulator, double *hdr, int worker_count); +void fast_psf_accumulator_set_verbose(FastPsfAccumulator *accumulator, + int verbose); void fast_psf_accumulator_destroy(FastPsfAccumulator *accumulator); void fast_psf_accumulator_report(const FastPsfAccumulator *accumulator, FILE *stream); diff --git a/tests/benchmark_fast_psf_fftw.c b/tests/benchmark_fast_psf_fftw.c new file mode 100644 index 0000000..1dcbba8 --- /dev/null +++ b/tests/benchmark_fast_psf_fftw.c @@ -0,0 +1,223 @@ +/* Resolver-only benchmark for the fast-mode FFTW global convolution. + * + * Fills one deterministic impulse buffer, then times the spatial reference and + * the FFTW production resolver on that identical buffer. First-use plan/setup + * is reported separately from steady-state execution. No catalog or geodesic + * work is involved. */ +#ifndef FAST_PSF_FFTW +#error "benchmark_fast_psf_fftw requires -DFAST_PSF_FFTW (PSF_BACKEND=cpu)" +#endif + +#include "fast_psf_fftw.h" +#include "optics.h" + +#include +#include +#include +#include +#include +#include +#include + +typedef struct { + int width, height, supersample, repeats; + int run_spatial; + int measure; + FastPsfDeposit deposit; +} Options; + +static int parse_int(const char *text, int *out) +{ + char *end = NULL; + const long value = strtol(text, &end, 10); + if (end == text || *end != '\0' || value <= 0 || value > 1000000) + return -1; + *out = (int)value; + return 0; +} + +static int parse_options(int argc, char **argv, Options *options) +{ + *options = (Options){.width = 320, + .height = 240, + .supersample = 2, + .repeats = 5, + .run_spatial = 0, + .measure = 0, + .deposit = FAST_PSF_DEPOSIT_NEAREST}; + for (int i = 1; i < argc; ++i) { + if (!strcmp(argv[i], "--spatial")) { + options->run_spatial = 1; + } else if (!strcmp(argv[i], "--measure")) { + options->measure = 1; + } else if (!strcmp(argv[i], "--bilinear")) { + options->deposit = FAST_PSF_DEPOSIT_BILINEAR; + } else if (!strcmp(argv[i], "--width") && i + 1 < argc) { + if (parse_int(argv[++i], &options->width)) + return -1; + } else if (!strcmp(argv[i], "--height") && i + 1 < argc) { + if (parse_int(argv[++i], &options->height)) + return -1; + } else if (!strcmp(argv[i], "--supersample") && i + 1 < argc) { + if (parse_int(argv[++i], &options->supersample)) + return -1; + } else if (!strcmp(argv[i], "--repeats") && i + 1 < argc) { + if (parse_int(argv[++i], &options->repeats)) + return -1; + } else { + fprintf(stderr, "unknown argument '%s'\n", argv[i]); + return -1; + } + } + return 0; +} + +static uint64_t hash_bytes(const double *values, size_t count) +{ + const unsigned char *bytes = (const unsigned char *)values; + const size_t nbytes = count * sizeof *values; + uint64_t hash = 1469598103934665603ULL; + for (size_t i = 0; i < nbytes; ++i) { + hash ^= bytes[i]; + hash *= 1099511628211ULL; + } + return hash; +} + +static double wall_seconds(void) +{ + return omp_get_wtime(); +} + +static int compare_double(const void *a, const void *b) +{ + const double x = *(const double *)a; + const double y = *(const double *)b; + return x < y ? -1 : x > y ? 1 : 0; +} + +int main(int argc, char **argv) +{ + Options options; + if (parse_options(argc, argv, &options)) { + fprintf(stderr, + "usage: %s [--width W] [--height H] [--supersample N] " + "[--repeats R] [--bilinear] [--spatial] [--measure]\n", + argv[0]); + return 2; + } + const PointSpreadFunction psf = {.fwhm_pixels = 2.7, .moffat_beta = 4.5}; + const double relative_tail = 1e-8; + if (options.repeats > 128) + options.repeats = 128; + + printf("impulse: %dx%d supersample=%d deposit=%s repeats=%d spatial=%d\n", + options.width, options.height, options.supersample, + options.deposit == FAST_PSF_DEPOSIT_NEAREST ? "nearest" : "bilinear", + options.repeats, options.run_spatial); + printf("plan_mode=%s\n", options.measure ? "measure" : "estimate"); + fast_psf_fftw_set_plan_mode(options.measure); + + FastPsfAccumulator acc = {0}; + const double init_start = wall_seconds(); + if (fast_psf_accumulator_init(&acc, options.width, options.height, + options.supersample, options.deposit, &psf, + relative_tail, 0.0, 1)) { + fputs("accumulator initialization failed\n", stderr); + return 1; + } + const double init_seconds = wall_seconds() - init_start; + printf("setup: init=%.6f s fftw_plan=%.6f s kernel_fft=%.6f s\n", + init_seconds, acc.fftw_setup_seconds, acc.fftw_kernel_seconds); + fast_psf_fftw_report(acc.fftw, stdout); + + const size_t ss_count = (size_t)acc.supersampled_width * + acc.supersampled_height; + uint64_t state = 0x9e3779b97f4a7c15ULL; + for (size_t i = 0; i < ss_count; ++i) { + state = state * 6364136223846793005ULL + 1442695040888963407ULL; + if ((state >> 58) != 0) + continue; + state = state * 6364136223846793005ULL + 1442695040888963407ULL; + acc.buffer[3 * i] = (double)(state >> 40) / (double)(1ULL << 24); + state = state * 6364136223846793005ULL + 1442695040888963407ULL; + acc.buffer[3 * i + 1] = (double)(state >> 40) / (double)(1ULL << 24); + state = state * 6364136223846793005ULL + 1442695040888963407ULL; + acc.buffer[3 * i + 2] = (double)(state >> 40) / (double)(1ULL << 24); + } + const size_t hdr_count = (size_t)options.width * options.height * 3; + const size_t buffer_count = ss_count * 3; + const uint64_t impulse_hash = hash_bytes(acc.buffer, buffer_count); + printf("impulse_hash=%016llx\n", (unsigned long long)impulse_hash); + + double *fftw_hdr = calloc(hdr_count, sizeof *fftw_hdr); + double *spatial_hdr = calloc(hdr_count, sizeof *spatial_hdr); + if (fftw_hdr == NULL || spatial_hdr == NULL) { + fputs("HDR allocation failed\n", stderr); + free(fftw_hdr); + free(spatial_hdr); + fast_psf_accumulator_destroy(&acc); + return 1; + } + + /* Warm-up: first execution includes any lazy per-frame allocation. */ + if (fast_psf_accumulator_resolve(&acc, fftw_hdr, omp_get_max_threads())) { + fputs("FFTW warm-up resolve failed\n", stderr); + return 1; + } + double *samples = calloc((size_t)options.repeats, sizeof *samples); + if (samples == NULL) { + fputs("sample allocation failed\n", stderr); + return 1; + } + for (int i = 0; i < options.repeats; ++i) { + memset(fftw_hdr, 0, hdr_count * sizeof *fftw_hdr); + const double start = wall_seconds(); + if (fast_psf_accumulator_resolve(&acc, fftw_hdr, omp_get_max_threads())) { + fputs("FFTW resolve failed\n", stderr); + return 1; + } + samples[i] = wall_seconds() - start; + } + const uint64_t fftw_hash = hash_bytes(fftw_hdr, hdr_count); + double sorted[128]; + memcpy(sorted, samples, (size_t)options.repeats * sizeof *sorted); + qsort(sorted, (size_t)options.repeats, sizeof *sorted, compare_double); + printf("fftw: min=%.6f s median=%.6f s max=%.6f s\n", sorted[0], + sorted[options.repeats / 2], sorted[options.repeats - 1]); + printf("fftw samples:"); + for (int i = 0; i < options.repeats; ++i) + printf(" %.6f", samples[i]); + printf("\nfftw_hdr_hash=%016llx\n", (unsigned long long)fftw_hash); + + if (options.run_spatial) { + memset(spatial_hdr, 0, hdr_count * sizeof *spatial_hdr); + const double start = wall_seconds(); + if (fast_psf_accumulator_resolve_spatial_reference( + &acc, spatial_hdr, omp_get_max_threads())) { + fputs("spatial resolve failed\n", stderr); + return 1; + } + const double spatial_seconds = wall_seconds() - start; + printf("spatial: %.6f s\n", spatial_seconds); + printf("spatial_hdr_hash=%016llx\n", + (unsigned long long)hash_bytes(spatial_hdr, hdr_count)); + double max_abs = 0.0, peak = 0.0; + for (size_t i = 0; i < hdr_count; ++i) { + peak = fmax(peak, fabs(spatial_hdr[i])); + max_abs = fmax(max_abs, fabs(fftw_hdr[i] - spatial_hdr[i])); + } + printf("difference: max_abs=%.3g peak=%.3g speedup=%.3fx\n", max_abs, peak, + spatial_seconds / sorted[options.repeats / 2]); + } + + struct rusage usage; + if (getrusage(RUSAGE_SELF, &usage) == 0) + printf("rss_max_kb=%ld\n", usage.ru_maxrss); + + free(samples); + free(fftw_hdr); + free(spatial_hdr); + fast_psf_accumulator_destroy(&acc); + return 0; +} diff --git a/tests/test_camera_cli.py b/tests/test_camera_cli.py index 5e6572b..30931c3 100644 --- a/tests/test_camera_cli.py +++ b/tests/test_camera_cli.py @@ -9,6 +9,7 @@ import tempfile import zlib BUILD = Path(sys.argv[1] if len(sys.argv) > 1 else 'build/Release').resolve() +TESTDIR = Path(sys.argv[2]).resolve() if len(sys.argv) > 2 else BUILD ENV = dict(os.environ, OMP_NUM_THREADS='4') @@ -67,6 +68,15 @@ with tempfile.TemporaryDirectory(prefix='gr-camera-cli-') as directory: run(binary, *common, '--output', path, *options) return image_payload(path) + # CPU fast-mode CLI smoke test: the FFTW resolve must run and report its + # one-time setup line. + if backend == 'minkowski': + fast_path = tmp / 'minkowski_fast.png' + fast = run(binary, *common, '--fast-mode', '--fast-supersample', 2, + '--output', fast_path) + assert 'Fast FFTW:' in fast.stderr, fast.stderr + assert image_payload(fast_path) + # Equivalent independently specified and inferred camera geometry. inferred = render('position', '--observer-position', -30, 0, 0) explicit = render('explicit', '--observer-position', -30, 0, 0, @@ -117,7 +127,7 @@ with tempfile.TemporaryDirectory(prefix='gr-camera-cli-') as directory: assert not missing_catalog.exists(), result.stderr assert 'PSF cache ready' not in result.stderr track = tmp / f'{backend}.csv' - run(BUILD / f'test_observer_{backend}', track) + run(TESTDIR / f'test_observer_{backend}', track) single_map, movie_map = tmp / 'single.grlens', tmp / 'movie.grlens' single = render('moving', '--observer-position', 3, -4, 5, '--observer-velocity', 0.2, -0.1, 0.3, diff --git a/tests/test_fast_psf_fftw.c b/tests/test_fast_psf_fftw.c new file mode 100644 index 0000000..4c68aae --- /dev/null +++ b/tests/test_fast_psf_fftw.c @@ -0,0 +1,378 @@ +/* Dedicated FFTW-versus-spatial fast-mode convolution regression. + * + * Both resolvers run on the identical impulse buffer and must agree to double + * rounding. The spatial resolver is the production reference, not a fallback. + * Compile only for the CPU PSF backend. */ +#ifndef FAST_PSF_FFTW +#error "test_fast_psf_fftw requires -DFAST_PSF_FFTW (PSF_BACKEND=cpu)" +#endif + +#include "optics.h" +#include "fast_psf_fftw.h" + +#include +#include +#include +#include +#include +#include + +#ifndef INT_MAX +#define INT_MAX 2147483647 +#endif + +static int g_failures = 0; + +static void fail(const char *what) +{ + fprintf(stderr, "FAIL: %s\n", what); + ++g_failures; +} + +static PointSpreadFunction default_psf(void) +{ + PointSpreadFunction psf = {.fwhm_pixels = 2.7, .moffat_beta = 4.5}; + return psf; +} + +static uint64_t hash_bytes(const double *values, size_t count) +{ + const unsigned char *bytes = (const unsigned char *)values; + const size_t nbytes = count * sizeof *values; + uint64_t hash = 1469598103934665603ULL; + for (size_t i = 0; i < nbytes; ++i) { + hash ^= bytes[i]; + hash *= 1099511628211ULL; + } + return hash; +} + +typedef enum { + SCENE_CENTER_WHITE, + SCENE_DISTINCT_RGB, + SCENE_EDGES, + SCENE_BOUNDARY_PHASES, + SCENE_MULTIPLE, + SCENE_DENSE_CELL, + SCENE_RANDOM, + SCENE_EMPTY, +} Scene; + +static void fill_scene(FastPsfAccumulator *acc, Scene scene, int width, + int height) +{ + const LinearRgb white = {1.0, 1.0, 1.0}; + switch (scene) { + case SCENE_CENTER_WHITE: + fast_psf_accumulator_deposit(acc, 0.5 * width, 0.5 * height, white, 1.0); + break; + case SCENE_DISTINCT_RGB: + fast_psf_accumulator_deposit(acc, 0.5 * width, 0.5 * height, + (LinearRgb){1.0, 0.25, 0.05}, 0.75); + break; + case SCENE_EDGES: + fast_psf_accumulator_deposit(acc, 0.3, 0.3, white, 1.0); + fast_psf_accumulator_deposit(acc, width - 0.7, 0.4, white, 1.0); + fast_psf_accumulator_deposit(acc, 0.5, height - 0.6, white, 1.0); + fast_psf_accumulator_deposit(acc, width - 0.4, height - 0.3, white, 1.0); + break; + case SCENE_BOUNDARY_PHASES: + fast_psf_accumulator_deposit(acc, 10.499, 12.499, white, 1.0); + fast_psf_accumulator_deposit(acc, 10.501, 12.501, white, 1.0); + break; + case SCENE_MULTIPLE: + fast_psf_accumulator_deposit(acc, 5.5, 5.5, (LinearRgb){1.0, 0.0, 0.0}, 0.5); + fast_psf_accumulator_deposit(acc, 12.25, 7.75, (LinearRgb){0.0, 1.0, 0.0}, 0.25); + fast_psf_accumulator_deposit(acc, 8.1, 14.9, (LinearRgb){0.0, 0.0, 1.0}, 1.5); + break; + case SCENE_DENSE_CELL: + for (int i = 0; i < 32; ++i) + fast_psf_accumulator_deposit(acc, 9.0 + 0.01 * i, 9.0 + 0.013 * i, + (LinearRgb){1.0, 0.5, 0.25}, 0.05); + break; + case SCENE_RANDOM: { + uint64_t state = 0x9e3779b97f4a7c15ULL; + const size_t count = (size_t)acc->supersampled_width * + acc->supersampled_height * 3; + for (size_t i = 0; i < count; ++i) { + state = state * 6364136223846793005ULL + 1442695040888963407ULL; + acc->buffer[i] += (double)(state >> 40) / (double)(1ULL << 24) - 0.5; + } + break; + } + case SCENE_EMPTY: + break; + } +} + +static int run_compare(const char *name, int width, int height, + int supersample, FastPsfDeposit deposit, Scene scene, + int prefill) +{ + FastPsfAccumulator acc = {0}; + const PointSpreadFunction psf = default_psf(); + const double relative_tail = 1e-8; + if (fast_psf_accumulator_init(&acc, width, height, supersample, deposit, &psf, + relative_tail, 0.0, 1)) { + fprintf(stderr, "FAIL: %s: accumulator init failed\n", name); + ++g_failures; + return -1; + } + if (!acc.fftw_enabled) { + fprintf(stderr, "FAIL: %s: FFTW path not enabled\n", name); + ++g_failures; + fast_psf_accumulator_destroy(&acc); + return -1; + } + const size_t hdr_count = (size_t)width * height * 3; + const size_t buffer_count = (size_t)acc.supersampled_width * + acc.supersampled_height * 3; + double *fftw_hdr = malloc(hdr_count * sizeof *fftw_hdr); + double *spatial_hdr = malloc(hdr_count * sizeof *spatial_hdr); + if (fftw_hdr == NULL || spatial_hdr == NULL) { + fprintf(stderr, "FAIL: %s: HDR allocation failed\n", name); + ++g_failures; + free(fftw_hdr); + free(spatial_hdr); + fast_psf_accumulator_destroy(&acc); + return -1; + } + for (size_t i = 0; i < hdr_count; ++i) + fftw_hdr[i] = spatial_hdr[i] = prefill ? 0.25 : 0.0; + fill_scene(&acc, scene, width, height); + const uint64_t buffer_before = hash_bytes(acc.buffer, buffer_count); + + if (fast_psf_accumulator_resolve(&acc, fftw_hdr, 4)) { + fprintf(stderr, "FAIL: %s: FFTW resolve failed\n", name); + ++g_failures; + goto cleanup; + } + if (fast_psf_accumulator_resolve_spatial_reference(&acc, spatial_hdr, 4)) { + fprintf(stderr, "FAIL: %s: spatial reference resolve failed\n", name); + ++g_failures; + goto cleanup; + } + if (hash_bytes(acc.buffer, buffer_count) != buffer_before) { + fprintf(stderr, "FAIL: %s: resolve mutated the impulse buffer\n", name); + ++g_failures; + goto cleanup; + } + + double peak = 0.0, max_abs = 0.0, max_rel = 0.0, sum_sq = 0.0; + double flux_fftw[3] = {0.0, 0.0, 0.0}; + double flux_spatial[3] = {0.0, 0.0, 0.0}; + size_t nan_count = 0; + for (size_t i = 0; i < hdr_count; ++i) { + const double reference = spatial_hdr[i]; + const double error = fabs(fftw_hdr[i] - reference); + peak = fmax(peak, fabs(reference)); + max_abs = fmax(max_abs, error); + sum_sq += error * error; + if (!isfinite(fftw_hdr[i]) || !isfinite(reference)) + ++nan_count; + flux_fftw[i % 3] += fftw_hdr[i]; + flux_spatial[i % 3] += reference; + } + /* Mixed absolute/relative acceptance: the absolute floor covers the tiny + * Moffat far-wing samples where the FFT and the direct sum disagree only by + * roundoff, while the relative term checks significant samples. A crop, + * wrap, channel, or normalization bug produces O(1) errors well above both. */ + const double rms = sqrt(sum_sq / (double)hdr_count); + const double abs_tol = 1e-10 * fmax(1.0, peak); + const double rel_tol = 1e-9; + const double rel_floor = 1e-6 * fmax(peak, 1.0); + double max_violation = 0.0; + for (size_t i = 0; i < hdr_count; ++i) { + const double reference = spatial_hdr[i]; + const double error = fabs(fftw_hdr[i] - reference); + max_violation = + fmax(max_violation, error - (abs_tol + rel_tol * fabs(reference))); + if (fabs(reference) > rel_floor) + max_rel = fmax(max_rel, error / fabs(reference)); + } + int worst = -1; + double worst_error = 0.0; + for (size_t i = 0; i < hdr_count; ++i) { + const double error = fabs(fftw_hdr[i] - spatial_hdr[i]); + if (error > worst_error) { + worst_error = error; + worst = (int)(i / 3); + } + } + double flux_rel = 0.0; + for (int c = 0; c < 3; ++c) { + const double denom = fmax(fabs(flux_spatial[c]), 1e-30); + flux_rel = fmax(flux_rel, fabs(flux_fftw[c] - flux_spatial[c]) / denom); + } + + if (nan_count != 0 || max_violation > 0.0 || max_rel > rel_tol || + flux_rel > 1e-10) { + fprintf(stderr, + "FAIL: %s: peak=%.6g max_abs=%.3g (tol %.3g) max_rel=%.3g " + "rms=%.3g flux_rel=%.3g nan=%zu worst_pixel=%d\n", + name, peak, max_abs, abs_tol, max_rel, rms, flux_rel, nan_count, + worst); + ++g_failures; + } else { + printf("ok %-28s peak=%.4g max_abs=%.3g max_rel=%.3g rms=%.3g\n", name, + peak, max_abs, max_rel, rms); + } + +cleanup: + free(fftw_hdr); + free(spatial_hdr); + fast_psf_accumulator_destroy(&acc); + return g_failures == 0 ? 0 : -1; +} + +static void test_next_smooth_size(void) +{ + struct { + size_t input; + long expected; /* -1 = must fail */ + } cases[] = { + {1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, + {7, 7}, {8, 8}, {9, 9}, {11, 12}, {13, 14}, {15, 15}, + {121, 125}, {127, 128}, {200, 200}, {241, 243}, + {(size_t)INT_MAX + 1, -1}, + {0, -1}, + }; + for (size_t i = 0; i < sizeof cases / sizeof cases[0]; ++i) { + size_t out = 0; + const int rc = fast_psf_fftw_next_smooth_size(cases[i].input, &out); + if (cases[i].expected < 0) { + if (rc == 0) + fprintf(stderr, "FAIL: next_smooth_size(%zu) unexpectedly returned %zu\n", + cases[i].input, out), ++g_failures; + } else if (rc != 0 || out != (size_t)cases[i].expected) { + fprintf(stderr, "FAIL: next_smooth_size(%zu) = rc %d, %zu (want %ld)\n", + cases[i].input, rc, out, cases[i].expected); + ++g_failures; + } + } + size_t out = 0; + if (fast_psf_fftw_next_smooth_size((size_t)1 << 30, &out) != 0 || + out > (size_t)INT_MAX) + fail("next_smooth_size near INT_MAX did not return a valid FFT size"); +} + +/* A single impulse must not wrap to the opposite edge, and an empty buffer must + * leave the framebuffer untouched. */ +static void test_no_wraparound_and_empty(void) +{ + PointSpreadFunction psf = default_psf(); + /* A compact kernel keeps the image much wider than the support so a genuine + * circular-wrap bug would show up as an opposite-edge ghost. */ + psf.fwhm_pixels = 0.02; + const double relative_tail = 1e-8; + const int width = 12, height = 12, supersample = 2; + FastPsfAccumulator acc = {0}; + if (fast_psf_accumulator_init(&acc, width, height, supersample, + FAST_PSF_DEPOSIT_NEAREST, &psf, relative_tail, + 0.0, 1)) { + fail("wraparound accumulator init"); + return; + } + const size_t count = (size_t)width * height * 3; + double *hdr = calloc(count, sizeof *hdr); + if (hdr == NULL) { + fail("wraparound allocation"); + fast_psf_accumulator_destroy(&acc); + return; + } + fast_psf_accumulator_deposit(&acc, 0.6, 0.6, (LinearRgb){1.0, 1.0, 1.0}, 1.0); + if (fast_psf_accumulator_resolve(&acc, hdr, 2)) { + fail("wraparound resolve"); + } else { + const int far_x = width - 1, far_y = height - 1; + const double ghost = hdr[3 * (far_y * width + far_x)]; + if (fabs(ghost) > 1e-15) + fprintf(stderr, "FAIL: opposite-edge ghost value %.3g\n", ghost), + ++g_failures; + if (hdr[0] <= 0.0) + fail("impulse peak missing near the deposited corner"); + } + /* Empty buffer: the same accumulator, after clearing, must not change HDR. */ + fast_psf_accumulator_clear(&acc); + double background[3] = {0.25, 0.5, 0.75}; + for (size_t i = 0; i < count; ++i) + hdr[i] = background[i % 3]; + if (fast_psf_accumulator_resolve(&acc, hdr, 2)) + fail("empty resolve"); + for (size_t i = 0; i < count; ++i) { + if (hdr[i] != background[i % 3]) { + fail("empty input changed the HDR framebuffer"); + break; + } + } + free(hdr); + fast_psf_accumulator_destroy(&acc); +} + +/* Channel separation: a pure-red impulse must leave green and blue at zero. */ +static void test_channel_isolation(void) +{ + const PointSpreadFunction psf = default_psf(); + const int width = 20, height = 18; + FastPsfAccumulator acc = {0}; + if (fast_psf_accumulator_init(&acc, width, height, 2, + FAST_PSF_DEPOSIT_NEAREST, &psf, 1e-8, 0.0, 1)) { + fail("channel isolation init"); + return; + } + const size_t count = (size_t)width * height * 3; + double *hdr = calloc(count, sizeof *hdr); + fast_psf_accumulator_deposit(&acc, 10.0, 9.0, (LinearRgb){1.0, 0.0, 0.0}, 1.0); + if (fast_psf_accumulator_resolve(&acc, hdr, 2)) + fail("channel isolation resolve"); + for (size_t i = 0; i < count; i += 3) { + if (hdr[i + 1] != 0.0 || hdr[i + 2] != 0.0) { + fail("green/blue channel leaked into a pure-red impulse"); + break; + } + } + free(hdr); + fast_psf_accumulator_destroy(&acc); +} + +int main(void) +{ + test_next_smooth_size(); + test_no_wraparound_and_empty(); + test_channel_isolation(); + + const int sizes[][2] = {{17, 13}, {16, 16}, {23, 31}, {33, 17}}; + const int supersamples[] = {1, 2, 3, 4}; + const FastPsfDeposit deposits[] = {FAST_PSF_DEPOSIT_NEAREST, + FAST_PSF_DEPOSIT_BILINEAR}; + for (size_t s = 0; s < sizeof sizes / sizeof sizes[0]; ++s) { + for (size_t n = 0; n < sizeof supersamples / sizeof supersamples[0]; ++n) { + for (size_t d = 0; d < sizeof deposits / sizeof deposits[0]; ++d) { + for (Scene scene = SCENE_CENTER_WHITE; scene <= SCENE_EMPTY; + ++scene) { + char name[128]; + snprintf(name, sizeof name, "%dx%d N%d %s scene%d", sizes[s][0], + sizes[s][1], supersamples[n], + deposits[d] == FAST_PSF_DEPOSIT_NEAREST ? "near" : "bilin", + (int)scene); + if (run_compare(name, sizes[s][0], sizes[s][1], supersamples[n], + deposits[d], scene, 0) != 0) { + /* Keep going: collect all failures before summarizing. */ + } + } + } + } + } + run_compare("prefilled background", 24, 20, 2, FAST_PSF_DEPOSIT_NEAREST, + SCENE_MULTIPLE, 1); + /* One axis (height) stays at the exact FFT size while the other grows. */ + run_compare("single-axis padding", 32, 8, 2, FAST_PSF_DEPOSIT_NEAREST, + SCENE_CENTER_WHITE, 0); + + if (g_failures != 0) { + fprintf(stderr, "%d fast-PSF FFTW regression failure(s)\n", g_failures); + return 1; + } + puts("fast PSF FFTW regressions passed"); + return 0; +} diff --git a/usage.md b/usage.md index 2710529..70f54d1 100644 --- a/usage.md +++ b/usage.md @@ -291,6 +291,20 @@ bright event whose requested support exceeds that radius is wing-clipped and counted in the report. `--psf-min-y` is still applied per event. Fast mode is a preview approximation, not the physically exact per-event PSF path. +In CPU builds the global convolution is a zero-padded FFTW linear convolution +(`fftw3`/`fftw3_omp`, double precision). The immutable kernel spectrum, FFTW +plans, and scratch buffers are built once when the accumulator is initialized +and reused for every frame; the disposable spatial reference remains only for +tests and benchmarks. Startup prints one `Fast FFTW:` line with the minimum +linear-convolution extent, the selected smooth FFT dimensions, worker count, +plan mode, one-time plan and kernel-transform time, and scratch bytes. +`--verbose` additionally prints per-frame stage timings (zero/pack, forward FFT, +frequency multiply, inverse FFT, crop/downsample, total). The default FFTW +planning mode is `FFTW_ESTIMATE`; `FFTW_MEASURE` can cut steady-state execution +substantially but costs a much longer one-time plan for large frames, as +recorded in +[benchmarks/fast_mode_fftw_2026-09-25.md](benchmarks/fast_mode_fftw_2026-09-25.md). + ## Progress and diagnostics The PSF-cache completion line is printed before tracing and catalog splatting