Why do both of these subroutines compile and build with no error?

brianreinhold · ‎09-17-2025

I am using the FFTW3 wrappers in the MKL for real to real 2d transforms. My needs are, every time step in a dynamic model, to go from spectral space to grid point space, perform some nonlinear computations, transform back to spectral space without using complex arrays (replacing the old F77 FFT991 FFT routines used in many meteorological models).

So step 1 is create a subroutine to generate the needed plans based upon the current zonal and meridional dimensions.

I wrote two such routines. BOTH COMPILE and BUILD (though I don't use the plans yet because I have no figured out how the arrays are actually used in the transform yet).

The module with the two subroutines are shown below:

module Fftw3HelpersModule
    use, intrinsic :: iso_c_binding
    implicit none

    ! FFTW MKL plan handles for different transforms
    integer*8 :: plan_spectral_to_grid_sin = 0
    integer*8 :: plan_grid_to_spectral_sin = 0
    integer*8 :: plan_spectral_to_grid_cos = 0
    integer*8 :: plan_grid_to_spectral_cos = 0
    integer*8 :: plan_spectral_to_grid_zonal_1d = 0
    integer*8 :: plan_grid_to_spectral_zonal_1d = 0
    
    ! FFTW plan handles for different transforms
    type(c_ptr) :: plan_spectral_to_grid_sin2, plan_grid_to_spectral_sin2
    type(c_ptr) :: plan_spectral_to_grid_cos2, plan_grid_to_spectral_cos2
    type(c_ptr) :: plan_spectral_to_grid_zonal_1d2, plan_grid_to_spectral_zonal_1d2

    contains

    subroutine setup_fftw_plans(Nx, Ny)
        use, intrinsic :: iso_c_binding
        implicit none

        include 'fftw3.f'
        include 'fftw3_mkl.f'

        integer(c_int), intent(in) :: Nx, Ny
        real(C_FLOAT), allocatable :: dummy(:, :)
        real(C_FLOAT), allocatable :: dummy_1d(:)

        ! Cleanup existing plans
        if (plan_spectral_to_grid_sin .NE. 0) call sfftw_destroy_plan(plan_spectral_to_grid_sin)
        if (plan_grid_to_spectral_sin .NE. 0) call sfftw_destroy_plan(plan_grid_to_spectral_sin)
        if (plan_spectral_to_grid_cos .NE. 0) call sfftw_destroy_plan(plan_spectral_to_grid_cos)
        if (plan_grid_to_spectral_cos .NE. 0) call sfftw_destroy_plan(plan_grid_to_spectral_cos)
        if (plan_spectral_to_grid_zonal_1d .NE. 0) call sfftw_destroy_plan(plan_spectral_to_grid_zonal_1d)
        if (plan_grid_to_spectral_zonal_1d .NE. 0) call sfftw_destroy_plan(plan_grid_to_spectral_zonal_1d)

        allocate(dummy(Ny, Nx))
        allocate(dummy_1d(Ny))

        ! Plans for sin(kx)sin(ly) and cos(kx)sin(ly) components
        call sfftw_plan_r2r_2d(plan_grid_to_spectral_sin, Ny, Nx, dummy, dummy, FFTW_RODFT00, FFTW_R2HC, FFTW_MEASURE) ! MKL_RODFT00
        call sfftw_plan_r2r_2d(plan_spectral_to_grid_sin, Ny, Nx, dummy, dummy, FFTW_RODFT00, FFTW_HC2R, FFTW_MEASURE)
        ! Plans for sin(kx)cos(ly) and cos(kx)cos(ly) components
        call sfftw_plan_r2r_2d(plan_grid_to_spectral_cos, Ny, Nx, dummy, dummy, FFTW_REDFT00, FFTW_R2HC, FFTW_MEASURE)
        call sfftw_plan_r2r_2d(plan_spectral_to_grid_cos, Ny, Nx, dummy, dummy, FFTW_REDFT00, FFTW_HC2R, FFTW_MEASURE)
        ! Plans for zonal flow components cos(ly).
        call sfftw_plan_r2r_1d(plan_grid_to_spectral_zonal_1d, Ny, dummy_1d, dummy_1d, FFTW_REDFT00, FFTW_MEASURE)
        call sfftw_plan_r2r_1d(plan_spectral_to_grid_zonal_1d, Ny, dummy_1d, dummy_1d, FFTW_REDFT00, FFTW_MEASURE)

        deallocate(dummy)
        deallocate(dummy_1d)
    end subroutine setup_fftw_plans
    
    ! WHy does this version compile and build???
    subroutine setup_fftw_plans_2(Nx, Ny)
        use, intrinsic :: iso_c_binding
        implicit none
        include 'fftw3.f03'

        integer(c_int), intent(in) :: Nx, Ny
        real(C_FLOAT), allocatable :: dummy(:, :)
        real(C_FLOAT), allocatable :: dummy_1d(:)

        ! Cleanup existing plans
        if (c_associated(plan_spectral_to_grid_sin2)) call fftw_destroy_plan(plan_spectral_to_grid_sin2)
        if (c_associated(plan_grid_to_spectral_sin2)) call fftw_destroy_plan(plan_grid_to_spectral_sin2)
        if (c_associated(plan_spectral_to_grid_cos2)) call fftw_destroy_plan(plan_spectral_to_grid_cos2)
        if (c_associated(plan_grid_to_spectral_cos2)) call fftw_destroy_plan(plan_grid_to_spectral_cos2)
        if (c_associated(plan_spectral_to_grid_zonal_1d2)) call fftw_destroy_plan(plan_spectral_to_grid_zonal_1d2)
        if (c_associated(plan_grid_to_spectral_zonal_1d2)) call fftw_destroy_plan(plan_grid_to_spectral_zonal_1d2)

        allocate(dummy(Ny, Nx))
        allocate(dummy_1d(Ny))

        ! Plans for sin(kx)sin(ly) and cos(kx)sin(ly) components
        plan_grid_to_spectral_sin2 = fftwf_plan_r2r_2d(Ny, Nx, dummy, dummy, FFTW_RODFT00, FFTW_R2HC, FFTW_MEASURE)
        plan_spectral_to_grid_sin2 = fftwf_plan_r2r_2d(Ny, Nx, dummy, dummy, FFTW_RODFT00, FFTW_HC2R, FFTW_MEASURE)
        ! Plans for sin(kx)cos(ly) and cos(kx)cos(ly) components
        plan_grid_to_spectral_cos2 = fftwf_plan_r2r_2d(Ny, Nx, dummy, dummy, FFTW_REDFT00, FFTW_R2HC, FFTW_MEASURE)
        plan_spectral_to_grid_cos2 = fftwf_plan_r2r_2d(Ny, Nx, dummy, dummy, FFTW_REDFT00, FFTW_HC2R, FFTW_MEASURE)
        ! Plans for zonal flow components cos(ly).
        plan_grid_to_spectral_zonal_1d2 = fftwf_plan_r2r_1d(Ny, dummy_1d, dummy_1d, FFTW_REDFT00, FFTW_MEASURE)
        plan_spectral_to_grid_zonal_1d2 = fftwf_plan_r2r_1d(Ny, dummy_1d, dummy_1d, FFTW_REDFT00, FFTW_MEASURE)

        deallocate(dummy)
        deallocate(dummy_1d)
    end subroutine setup_fftw_plans_2

end module Fftw3HelpersModule

Why do they both build? I would expect that subroutine setup_fftw_plans(Nx, Ny) would compile as it is based upon the provided examples, and I would expect subroutine setup_fftw_plans2(Nx, Ny) to fail as that is based upon using the FFTW3 library. But both succeed.

I am using visual studio with the Intel compiler and the only item I have configured in the project properties is to use the MKL library (parallel). Additional library and include directories are empty as is the additional dependencies.

How and why do both build? Which one should I use?

Thanks for any help/clarification!

brianreinhold · ‎09-18-2025

Apparently if you include 'fftw3.f' and 'fftw3_mkl.f' you are using an old F77 compatible version of the FFTW3 wrappers. Unfortunately, the examples provided are based on these outdated interfaces. If, instead you include 'fftw3.f03' you are using the more modern F2003 interfaces. There are no provided examples for this more modern version, but in theory, I am assuming that is what one should use unless you are using the far less safe F77 fortran.

Fengrui · ‎09-18-2025

Both are supported in oneMKL fftw3 wrappers. setup_fftw_plans (fftw3_mkl.f, fftw3_mkl_f77.h) are defined in the F77 style while setup_fftw_plans2 (fftw3.f03) are defined in a modern Fortran style.

Fengrui · ‎09-18-2025

It looks you already figured that out. Didn't refresh the page for a while:)