Merge pull request #3458 from andrew-platt/b/fftpack_oneapi

FFTPACK 5.1: suppress default-real promotion for Intel and Flang builds

Author: andrew-platt

modules/nwtc-library/CMakeLists.txt

@@ -99,6 +99,7 @@ set(NWTCLIBS_SOURCES
 
   # NetLib sources
   src/NetLib/fftpack/fftpack5.1.f
+  src/NetLib/fftpack/fftpack_kind.f
   src/NetLib/scalapack/dlasrt2.f
   src/NetLib/scalapack/slasrt2.f
   src/NetLib/fftpack/NWTC_FFTPACK.f90
@@ -126,11 +127,49 @@ get_filename_component(FCNAME ${CMAKE_Fortran_COMPILER} NAME)
 
 # FFTPACK 5.1 uses legacy Fortran 77 type punning (COMPLEX<->REAL) and must
 # NOT be promoted to double precision (its callers pass explicit SiKi/R4Ki arrays).
-# -Wno-pedantic is required because -pedantic (added for Debug builds) promotes the
-# argument mismatches back to errors, overriding -fallow-argument-mismatch.
+#
+# fftpack5.1.f declares all of its arrays as bare REAL/COMPLEX, but NWTC_FFTPACK.f90
+# hands it explicitly kinded REAL(SiKi)/COMPLEX(SiKi) buffers and passes their element
+# counts as LENSAV/LENWRK.  If DOUBLE_PRECISION promotes the default REAL to 8 bytes in
+# these files, FFTPACK writes twice as many bytes as those buffers hold, silently
+# corrupting the heap and stack (observed as SIGBUS/SIGSEGV far from the FFT call).
+# fftpack_kind.f reports the kind actually used so NWTC_FFTPACK can verify it at run
+# time; it must therefore be compiled with exactly the same flags as fftpack5.1.f.
+#
+# The Visual Studio build does the same thing via RealKIND="realKIND4" on both files
+# in vs-build/modules/NWTC-Library.vfproj.
+set(FFTPACK_SOURCES
+   src/NetLib/fftpack/fftpack5.1.f
+   src/NetLib/fftpack/fftpack_kind.f
+)
 if (${CMAKE_Fortran_COMPILER_ID} STREQUAL "GNU")
-   set_source_files_properties(src/NetLib/fftpack/fftpack5.1.f PROPERTIES COMPILE_FLAGS
+   # -Wno-pedantic is required because -pedantic (added for Debug builds) promotes the
+   # argument mismatches back to errors, overriding -fallow-argument-mismatch.
+   set_source_files_properties(${FFTPACK_SOURCES} PROPERTIES COMPILE_FLAGS
       "-fallow-argument-mismatch -Wno-pedantic -fno-default-real-8 -fno-default-double-8")
+elseif (${CMAKE_Fortran_COMPILER_ID} MATCHES "^Intel")
+   # Intel's -real-size 64 promotes REAL->REAL(8) *and* COMPLEX->COMPLEX(8), so it must be
+   # undone here.  -double-size 64 restores the Intel default so the genuine DOUBLE
+   # PRECISION accumulators in fftpack5.1.f (DSUM, TPI, ARGH, ...) stay 8 bytes, matching
+   # gfortran's -fno-default-double-8.  Intel needs no -fallow-argument-mismatch analogue:
+   # these are F77 externals with no explicit interfaces, so it never diagnoses them.
+   if("${CMAKE_Fortran_COMPILER_VERSION}" VERSION_GREATER "19")
+      if (WIN32)
+         set(FFTPACK_REAL4_FLAGS "/real-size:32 /double-size:64")
+      else()
+         set(FFTPACK_REAL4_FLAGS "-real-size 32 -double-size 64")
+      endif()
+   else()   # the hyphenated spelling above is only supported from version 19 onwards
+      if (WIN32)
+         set(FFTPACK_REAL4_FLAGS "/real_size:32 /double_size:64")
+      else()
+         set(FFTPACK_REAL4_FLAGS "-real_size 32 -double_size 64")
+      endif()
+   endif()
+   set_source_files_properties(${FFTPACK_SOURCES} PROPERTIES COMPILE_FLAGS "${FFTPACK_REAL4_FLAGS}")
+elseif (${CMAKE_Fortran_COMPILER_ID} STREQUAL "Flang")
+   set_source_files_properties(${FFTPACK_SOURCES} PROPERTIES COMPILE_FLAGS
+      "-fno-default-real-8")
 endif()
 
 # Recursive use of routine in qk61/dqk61 will trigger errors in debug

modules/nwtc-library/src/NetLib/fftpack/NWTC_FFTPACK.f90

@@ -7,6 +7,7 @@
    !                  Also updated to check that transform has been initialized for the
    !                    correct type (to avoid having wSave too small)
    ! ADP: 07/28/2014: Added in the complex FFT routines from fftpack v. 4.1
+   ! ADP: 08/15/2026: upgraded from fftpack v4.1 to v5.1 and added interfaces for 2D fft
 !=======================================================================
 MODULE NWTC_FFTPACK
 !-----------------------------------------------------------------------
@@ -674,6 +675,38 @@ SUBROUTINE ExitSINT(FFT_Data, ErrStat)
 
    END SUBROUTINE ExitSINT
   !------------------------------------------------------------------------
+   SUBROUTINE CheckFFTPACKRealKind( ErrStat )
+
+        ! This subroutine verifies that fftpack5.1.f was compiled with a default REAL
+        ! kind of SiKi.  FFTPACK 5.1 declares its arrays as bare REAL/COMPLEX, but this
+        ! wrapper hands it explicitly kinded REAL(SiKi)/COMPLEX(SiKi) buffers and passes
+        ! their element counts as LENSAV/LENWRK.  If the build promotes the default REAL
+        ! to 8 bytes in fftpack5.1.f (DOUBLE_PRECISION does this by default, via
+        ! -fdefault-real-8 for GNU or -real-size 64 for Intel), FFTPACK writes twice as
+        ! many bytes as those buffers hold and silently corrupts the heap and stack.
+        ! See the FFTPACK_SOURCES block in modules/nwtc-library/CMakeLists.txt.
+
+      IMPLICIT                         NONE
+
+      INTEGER, INTENT(OUT),OPTIONAL :: ErrStat        ! returns non-zero if an error occurred
+
+      INTEGER, EXTERNAL             :: FFTPACK_REALKIND  ! from src/NetLib/fftpack/fftpack_kind.f
+
+
+      IF ( PRESENT(ErrStat) ) ErrStat = ErrID_None
+
+      IF ( FFTPACK_REALKIND() /= SiKi ) THEN
+         CALL ProgAbort ( 'FFTPACK 5.1 was compiled with a default REAL kind of '// &
+                          TRIM(Num2LStr(FFTPACK_REALKIND()))//', but NWTC_FFTPACK requires '// &
+                          TRIM(Num2LStr(SiKi))//'.  The build must suppress default-real promotion '// &
+                          'for fftpack5.1.f and fftpack_kind.f.', PRESENT(ErrStat) )
+         IF ( PRESENT(ErrStat) ) ErrStat = ErrID_Fatal
+      ENDIF
+
+
+   END SUBROUTINE CheckFFTPACKRealKind
+  !------------------------------------------------------------------------
+
    SUBROUTINE InitCOST( NumSteps, FFT_Data, NormalizeIn, ErrStat )
 
         ! This subroutine initializes the cosine transform working space
@@ -694,6 +727,14 @@ SUBROUTINE InitCOST( NumSteps, FFT_Data, NormalizeIn, ErrStat )
 
       IF ( PRESENT(ErrStat) ) ErrStat = ErrID_None
 
+        ! Verify FFTPACK's default REAL kind matches this wrapper's (SiKi)
+
+      CALL CheckFFTPACKRealKind( ErrStat )
+      IF ( PRESENT(ErrStat) ) THEN
+         IF ( ErrStat >= AbortErrLev ) RETURN
+      ENDIF
+
+
         ! Number of timesteps in the time series returned from the cosine transform
         ! N should be odd:
 
@@ -763,6 +804,14 @@ SUBROUTINE InitCFFT( NumSteps, FFT_Data, NormalizeIn, ErrStat )
 
       IF ( PRESENT(ErrStat) ) ErrStat = ErrID_None
 
+        ! Verify FFTPACK's default REAL kind matches this wrapper's (SiKi)
+
+      CALL CheckFFTPACKRealKind( ErrStat )
+      IF ( PRESENT(ErrStat) ) THEN
+         IF ( ErrStat >= AbortErrLev ) RETURN
+      ENDIF
+
+
         ! Number of timesteps in the time series returned from the backward FFT
         ! N should be even:
 
@@ -832,6 +881,14 @@ SUBROUTINE InitFFT( NumSteps, FFT_Data, NormalizeIn, ErrStat )
 
       IF ( PRESENT(ErrStat) ) ErrStat = ErrID_None
 
+        ! Verify FFTPACK's default REAL kind matches this wrapper's (SiKi)
+
+      CALL CheckFFTPACKRealKind( ErrStat )
+      IF ( PRESENT(ErrStat) ) THEN
+         IF ( ErrStat >= AbortErrLev ) RETURN
+      ENDIF
+
+
         ! Number of timesteps in the time series returned from the backward FFT
         ! N should be even:
 
@@ -902,6 +959,14 @@ SUBROUTINE InitSINT( NumSteps, FFT_Data, NormalizeIn, ErrStat )
 
       IF ( PRESENT(ErrStat) ) ErrStat = ErrID_None
 
+        ! Verify FFTPACK's default REAL kind matches this wrapper's (SiKi)
+
+      CALL CheckFFTPACKRealKind( ErrStat )
+      IF ( PRESENT(ErrStat) ) THEN
+         IF ( ErrStat >= AbortErrLev ) RETURN
+      ENDIF
+
+
         ! Number of timesteps in the time series returned from the sine transform
         ! N should be odd:
 
@@ -994,6 +1059,14 @@ SUBROUTINE InitFFT2D( L, M, FFT_Data, NormalizeIn, ErrStat )
 
       IF ( PRESENT(ErrStat) ) ErrStat = ErrID_None
 
+        ! Verify FFTPACK's default REAL kind matches this wrapper's (SiKi)
+
+      CALL CheckFFTPACKRealKind( ErrStat )
+      IF ( PRESENT(ErrStat) ) THEN
+         IF ( ErrStat >= AbortErrLev ) RETURN
+      ENDIF
+
+
       FFT_Data%L = L
       FFT_Data%M = M
 
@@ -1155,6 +1228,14 @@ SUBROUTINE InitCFFT2D( L, M, FFT_Data, NormalizeIn, ErrStat )
 
       IF ( PRESENT(ErrStat) ) ErrStat = ErrID_None
 
+        ! Verify FFTPACK's default REAL kind matches this wrapper's (SiKi)
+
+      CALL CheckFFTPACKRealKind( ErrStat )
+      IF ( PRESENT(ErrStat) ) THEN
+         IF ( ErrStat >= AbortErrLev ) RETURN
+      ENDIF
+
+
       FFT_Data%L = L
       FFT_Data%M = M
 

modules/nwtc-library/src/NetLib/fftpack/fftpack_kind.f

@@ -0,0 +1,29 @@
+C=======================================================================
+C  FFTPACK_REALKIND
+C
+C  Returns the default REAL kind that this file -- and therefore
+C  fftpack5.1.f -- was compiled with.
+C
+C  This file MUST be compiled with exactly the same real-size flags as
+C  fftpack5.1.f.  See the FFTPACK_SOURCES block in
+C  modules/nwtc-library/CMakeLists.txt and the RealKIND="realKIND4"
+C  file configurations in vs-build/modules/NWTC-Library.vfproj.
+C
+C  FFTPACK 5.1 declares its arrays as bare REAL/COMPLEX, while the
+C  NWTC_FFTPACK wrapper passes it explicitly kinded REAL(SiKi) and
+C  COMPLEX(SiKi) buffers together with their element counts (LENSAV,
+C  LENWRK).  A build that promotes the default REAL to 8 bytes in
+C  fftpack5.1.f -- which DOUBLE_PRECISION does by default via
+C  -fdefault-real-8 (GNU) or -real-size 64 (Intel) -- makes FFTPACK
+C  write twice as many bytes as those buffers hold.  Nothing diagnoses
+C  that at compile time, and the resulting heap/stack corruption
+C  surfaces as a SIGBUS or SIGSEGV far away from the FFT call.
+C
+C  NWTC_FFTPACK calls this at initialization and aborts with a clear
+C  message if the answer is not SiKi, so a build system that forgets
+C  the flag fails loudly instead of corrupting memory.
+C=======================================================================
+      INTEGER FUNCTION FFTPACK_REALKIND ()
+      FFTPACK_REALKIND = KIND(1.0)
+      RETURN
+      END

modules/nwtc-library/tests/test_NWTC_FFTPACK.F90

@@ -17,6 +17,7 @@ module test_NWTC_FFTPACK
 subroutine test_NWTC_FFTPACK_suite(testsuite)
    type(unittest_type), allocatable, intent(out) :: testsuite(:)
    testsuite = [ &
+      new_unittest("FFTPACK_real_kind", test_fftpack_real_kind), &
       new_unittest("FFT_roundtrip", test_fft_roundtrip), &
       new_unittest("FFT_forward_known", test_fft_forward_known), &
       new_unittest("FFT_forward_sign", test_fft_forward_sign), &
@@ -32,6 +33,17 @@ subroutine test_NWTC_FFTPACK_suite(testsuite)
    ]
 end subroutine
 
+! FFTPACK must be compiled with 4-byte default REAL, even in a DOUBLE_PRECISION
+! build: the NWTC_FFTPACK wrapper hands it REAL(SiKi) wSave/wWork buffers sized in
+! 4-byte elements.  A build that forgot to suppress default-real promotion for
+! fftpack5.1.f would have FFTPACK write 8-byte elements into them and corrupt memory.
+subroutine test_fftpack_real_kind(error)
+   type(error_type), allocatable, intent(out) :: error
+   integer, external :: FFTPACK_REALKIND
+
+   call check(error, FFTPACK_REALKIND(), SiKi)
+end subroutine
+
 ! Forward then backward (no normalization) gives x*N
 subroutine test_fft_roundtrip(error)
    type(error_type), allocatable, intent(out) :: error

vs-build/modules/NWTC-Library.vfproj

@@ -122,6 +122,32 @@
 				<File RelativePath="..\..\modules\nwtc-library\src\NetLib\slatec\dqk61.f"/>
 				<File RelativePath="..\..\modules\nwtc-library\src\NetLib\slatec\fdump.f"/>
 				<File RelativePath="..\..\modules\nwtc-library\src\NetLib\fftpack\fftpack5.1.f">
+					<FileConfiguration Name="Double_Debug|x64">
+						<Tool Name="VFFortranCompilerTool" WarnInterfaces="false" RealKIND="realKIND4"/>
+					</FileConfiguration>
+					<FileConfiguration Name="Double_Release|x64">
+						<Tool Name="VFFortranCompilerTool" WarnInterfaces="false" RealKIND="realKIND4"/>
+					</FileConfiguration>
+					<FileConfiguration Name="Double_OpenMP_Release|x64">
+						<Tool Name="VFFortranCompilerTool" WarnInterfaces="false" RealKIND="realKIND4"/>
+					</FileConfiguration>
+					<FileConfiguration Name="Matlab_Release|x64">
+						<Tool Name="VFFortranCompilerTool" WarnInterfaces="false" RealKIND="realKIND4"/>
+					</FileConfiguration>
+					<FileConfiguration Name="Matlab_Debug|x64">
+						<Tool Name="VFFortranCompilerTool" WarnInterfaces="false" RealKIND="realKIND4"/>
+					</FileConfiguration>
+					<FileConfiguration Name="Release|x64">
+						<Tool Name="VFFortranCompilerTool" WarnInterfaces="false" RealKIND="realKIND4"/>
+					</FileConfiguration>
+					<FileConfiguration Name="Debug|x64">
+						<Tool Name="VFFortranCompilerTool" WarnInterfaces="false" RealKIND="realKIND4"/>
+					</FileConfiguration>
+					<FileConfiguration Name="OpenMP_Release|x64">
+						<Tool Name="VFFortranCompilerTool" WarnInterfaces="false" RealKIND="realKIND4"/>
+					</FileConfiguration>
+				</File>
+				<File RelativePath="..\..\modules\nwtc-library\src\NetLib\fftpack\fftpack_kind.f">
 					<FileConfiguration Name="Double_Debug|x64">
 						<Tool Name="VFFortranCompilerTool" RealKIND="realKIND4"/>
 					</FileConfiguration>
@@ -218,6 +244,8 @@
 			<Filter Name="ranlux">
 				<File RelativePath="..\..\modules\nwtc-library\src\ranlux\RANLUX.f90"/>
 			</Filter>
+			<File RelativePath="..\..\modules\nwtc-library\src\GridInterp.f90"/>
+			<File RelativePath="..\..\modules\nwtc-library\src\GridInterp_Types.f90"/>
 			<File RelativePath="..\..\modules\nwtc-library\src\JSON.f90"/>
 			<File RelativePath="..\..\modules\nwtc-library\src\KdTree.f90"/>
 			<File RelativePath="..\..\modules\nwtc-library\src\ModMesh.f90"/>
@@ -257,8 +285,6 @@
 			</File>
 			<File RelativePath="..\..\modules\nwtc-library\src\VTK.f90"/>
 			<File RelativePath="..\..\modules\nwtc-library\src\YAML.f90"/>
-			<File RelativePath="..\..\modules\nwtc-library\src\GridInterp.f90"/>
-			<File RelativePath="..\..\modules\nwtc-library\src\GridInterp_Types.f90"/>
 		</Filter>
 	</Files>
 	<Globals/>