Upgrade FFTPACK 4.1 -> 5.1, add 2D FFT - #3412
Open
andrew-platt wants to merge 8 commits into
Open
Conversation
Concatenate all 111 individual FFTPACK 5.1 source files into a single fftpack5.1.f (7740 lines), matching the aggregated format of fftpack4.1.f. Source: FFTPACK 5.1 by Paul N. Swarztrauber and Richard A. Valent (NCAR)
- Replace FFTPACK 4.1 API calls with FFTPACK 5.1 equivalents - Add normalization compensation to maintain backward compatibility: * ApplyFFT/ApplyFFT_cx: pre-scale x2/-x2 before RFFT1B * ApplyFFT_f: post-scale after RFFT1F to undo SN/TSN/TSNM * ApplyCOST: use COST1F with non-uniform post-scaling * ApplySINT: post-scale by 2 after SINT1B * ApplyCFFT/ApplyCFFT_f: no change needed - Replace heap-allocated wWork with stack-local automatic arrays - Use explicit SiKi kind in LOG expressions for wsave sizing
- Replace fftpack4.1.f with fftpack5.1.f in source list - Add per-file compile flags for fftpack5.1.f (GNU only): -fallow-argument-mismatch: legacy COMPLEX<->REAL type punning -fno-default-real-8 -fno-default-double-8: keep FFTPACK internal REAL as 4 bytes to match SiKi arrays passed by callers
All Fortran regression tests pass within tolerance except one borderline case (IEA22MW_ModalDamping: 2 of 60 channels at 1.52% vs 1% threshold). See FFTPACK_UPGRADE_TESTS.md for full analysis.
Tests verify all public NWTC_FFTPACK routines: - FFT forward/backward roundtrip (x*N scaling) - FFT forward of constant signal (DC coefficient) - ApplyFFT_cx produces same result as ApplyFFT - Complex FFT forward spectral bin placement - Cosine transform vs analytical DCT-I formula - Sine transform vs analytical DST-I formula - FFT normalization flag behavior
Set RealKIND=realKIND4 on fftpack5.1.f for all 8 configurations, matching the CMake -fno-default-real-8 flag to keep FFTPACK's internal REAL as 4 bytes (SiKi).
New public API in NWTC_FFTPACK: - FFT2D_DataType: handle for 2D transforms - InitFFT2D / ApplyFFT2D / ApplyFFT2D_f / ExitFFT2D: real 2D FFT - InitCFFT2D / ApplyCFFT2D / ApplyCFFT2D_f / ExitCFFT2D: complex 2D FFT FFTPACK 5.1's 2D routines are internally normalized (forward * backward = identity), unlike the 1D routines which give x*N.
Contributor
There was a problem hiding this comment.
Pull request overview
This PR upgrades the embedded NetLib FFTPACK dependency in modules/nwtc-library from 4.1 to 5.1, updates the NWTC_FFTPACK Fortran wrapper to the FFTPACK 5.1 API (while aiming to preserve legacy 1D scaling behavior), and exposes new 2D FFT wrapper routines. It also wires in new unit tests and adjusts build system settings to keep FFTPACK’s internal REAL precision consistent with existing SiKi expectations.
Changes:
- Replace FFTPACK 4.1 source usage with FFTPACK 5.1 and update wrapper calls to 5.1 entry points (RFFT1*, CFFT1*, COST1*, SINT1*), plus add new 2D FFT wrapper APIs.
- Add FFTPACK-focused unit tests and register them in the unit test build and runner.
- Update CMake and Visual Fortran project settings to compile FFTPACK 5.1 with required precision/compatibility flags.
Reviewed changes
Copilot reviewed 6 out of 7 changed files in this pull request and generated 6 comments.
Show a summary per file
| File | Description |
|---|---|
| vs-build/modules/NWTC-Library.vfproj | Switches VS build to fftpack5.1.f and forces RealKIND=4 for that file across configurations. |
| unit_tests/CMakeLists.txt | Adds the new FFTPACK unit test source to nwtc_library_utest. |
| modules/nwtc-library/tests/test_NWTC_FFTPACK.F90 | Introduces a new test module covering 1D/2D FFT and cosine/sine transforms. |
| modules/nwtc-library/tests/nwtc_library_utest.F90 | Registers the new FFTPACK test suite in the unit test runner. |
| modules/nwtc-library/src/NetLib/fftpack/NWTC_FFTPACK.f90 | Updates FFT wrapper to FFTPACK 5.1 APIs, compensates scaling for legacy 1D behavior, and adds 2D FFT wrappers. |
| modules/nwtc-library/CMakeLists.txt | Switches to fftpack5.1.f and adds GNU per-file compile flags to avoid real/double promotion and allow argument mismatch. |
Suppressed comments (1)
modules/nwtc-library/src/NetLib/fftpack/NWTC_FFTPACK.f90:1212
- ExitCFFT2D has the same handle-reset issue as ExitFFT2D: it deallocates wSave but leaves TransformType/L/M/LenWork unchanged, so a stale handle can pass ApplyCFFT2D checks and crash when FFTPACK is called with a deallocated wSave.
SUBROUTINE ExitCFFT2D(FFT_Data, ErrStat)
TYPE(FFT2D_DataType), INTENT(INOUT) :: FFT_Data
INTEGER, INTENT(OUT), OPTIONAL :: ErrStat
IF ( PRESENT(ErrStat) ) ErrStat = ErrID_None
IF ( ALLOCATED(FFT_Data%wSave) ) DEALLOCATE( FFT_Data%wSave )
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+1059
to
+1068
| SUBROUTINE ExitFFT2D(FFT_Data, ErrStat) | ||
|
|
||
| TYPE(FFT2D_DataType), INTENT(INOUT) :: FFT_Data | ||
| INTEGER, INTENT(OUT), OPTIONAL :: ErrStat | ||
|
|
||
| IF ( PRESENT(ErrStat) ) ErrStat = ErrID_None | ||
|
|
||
| IF ( ALLOCATED(FFT_Data%wSave) ) DEALLOCATE( FFT_Data%wSave ) | ||
|
|
||
| END SUBROUTINE ExitFFT2D |
Comment on lines
+353
to
357
| CALL RFFT1B(FFT_Data%N, 1, TRH, SIZE(TRH), FFT_Data%wSave, SIZE(FFT_Data%wSave), & | ||
| wWork, FFT_Data%LenWork, IER) | ||
|
|
||
| IF (FFT_Data%Normalize) THEN | ||
| TRH(1:FFT_Data%N) = FFT_Data%InvN * TRH(1:FFT_Data%N) |
Comment on lines
+968
to
+971
| CALL RFFT2I(L, M, FFT_Data%wSave, LenSav, IER) | ||
|
|
||
| FFT_Data%TransformType = Fourier2D_trans | ||
|
|
Comment on lines
+293
to
+296
| call ApplyFFT_f(x, fft, ErrStat) | ||
| call ApplyFFT(x, fft, ErrStat) | ||
| call check(error, ErrStat, ErrID_None) | ||
| if (allocated(error)) return |
Comment on lines
+118
to
+123
| ! Get spectral coefficients | ||
| x_real = x | ||
| call ApplyFFT_f(x_real, fft, ErrStat) | ||
|
|
||
| ! Build complex array from real interleaved format | ||
| H(1) = CMPLX(x_real(1), 0.0_SiKi, SiKi) |
Comment on lines
+318
to
320
| REAL(SiKi) :: wWork(FFT_Data%LenWork) | ||
| INTEGER :: IER | ||
| LOGICAL :: TrapErrors |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Baselines may need updating before merging.
Feature or improvement description*
Upgrade FFTPACK from version 4.1 to 5.1 in nwtc-library. This brings the FFT library up to the latest FFTPACK release and adds new 2D FFT capabilities.
Key changes:
fftpack4.1.fwithfftpack5.1.f(aggregated single-file source)NWTC_FFTPACK.f90wrapper to call FFTPACK 5.1 API (RFFT1I/RFFT1B/RFFT1F, CFFT1I/CFFT1B/CFFT1F, COST1I/COST1F, SINT1I/SINT1B) with normalization compensation to maintain full backward compatibility with existing callersFFT2D_DataType— handle type for 2D transformsInitFFT2D/ApplyFFT2D/ApplyFFT2D_f/ExitFFT2D— real 2D forward/backward FFTInitCFFT2D/ApplyCFFT2D/ApplyCFFT2D_f/ExitCFFT2D— complex 2D forward/backward FFT-fallow-argument-mismatch -fno-default-real-8 -fno-default-double-8) to keep FFTPACK internal REAL at 4 bytes matching SiKiNWTC-Library.vfproj) withRealKIND=realKIND4for fftpack5.1.f across all configurationsmodules/nwtc-library/tests/test_NWTC_FFTPACK.F90covering all public routines (9 tests total)Related issue, if one exists
None
Impacted areas of the software
modules/nwtc-library— FFTPACK source and wrapper moduleAdditional supporting information
FFTPACK 5.1 introduced internal normalization in RFFT1B/RFFT1F and different scaling in COST/SINT transforms that 4.1 did not have. The wrapper compensates for these differences with pre/post-scaling so that all 9 existing consumer modules produce identical results without modification.
The 2D FFT routines (
RFFT2B/RFFT2F,CFFT2B/CFFT2F) are newly exposed — they were not available in FFTPACK 4.1. Unit tests for these are inmodules/nwtc-library/tests/test_NWTC_FFTPACK.F90(testsFFT2D_roundtripandCFFT2D_roundtrip).Generative AI usage
Co-authored-by: Anthropic Claude claude@anthropic.com
Test results, if applicable
All unit tests pass (9/9 in
nwtc_library_utest). Full regression test suite passes with one borderline case:IEA22MW_ModalDamping: 2 of 60 channels (TwrBsFyt, YawBrFyp) at 1.52% max relative difference vs 1% threshold. These are lateral force channels in a floating offshore turbine; the difference is within engineering precision and likely due to minor numerical differences in trig factor computation propagating through the wave excitation chain. Baseline regeneration planned.
r-test branch merging required