From dae6879948f2f8e431b866a335395df9154cf33f Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sun, 14 Jun 2026 09:42:12 +0200 Subject: [PATCH 1/7] Replace GSL Bessel I_n with fortnum bessel_in_array Drop the gsl_sf_bessel_In_array C binding in mephit_pert and compute the modified Bessel functions through fortnum's bessel_in_array. fortnum fills I_0..I_nmax in one pass, so the orders |m|-1, |m|, |m|+1 are recovered with the symmetry I_{-n}(x) = I_n(x); for pol_mode = 0 the order -1 maps to I_1. Wire fortnum via CMake FetchContent at the pinned tag, guarded so a build that also pulls fortnum transitively through libneo does not declare the target twice. Add the fortnum include directory for the C/C++ sources and link the fortnum target. GSL stays linked until the remaining call sites move over. --- CMakeLists.txt | 14 +++++++++++++- src/mephit_pert.f90 | 28 ++++++++-------------------- 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f178408..cb7835f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -100,6 +100,17 @@ endif() include(Util) find_or_fetch(libneo) +### fortnum: numerical core replacing the former GSL routines. libneo may pull +### fortnum in transitively, so guard against declaring the target twice. +include(FetchContent) +if(NOT TARGET fortnum) + FetchContent_Declare(fortnum + GIT_REPOSITORY git@github.com:lazy-fortran/fortnum.git + GIT_TAG a7faa3c + ) + FetchContent_MakeAvailable(fortnum) +endif() + ### Project source files # SuiteSparse set(SUITESPARSE_SRC_FILES @@ -152,7 +163,7 @@ if(WITH_MFEM) ) endif() set_source_files_properties(src/mephit_run.c ${MEPHIT_C_SRC_FILES} ${MEPHIT_CPP_SRC_FILES} - PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} -DREAL=double -I${TRIANGLE_INCLUDE_DIR} -I${SUITESPARSE_INCLUDE_DIRS} -L${TRIANGLE_LIB_DIR}") + PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} -DREAL=double -I${TRIANGLE_INCLUDE_DIR} -I${SUITESPARSE_INCLUDE_DIRS} -I${fortnum_SOURCE_DIR}/include -L${TRIANGLE_LIB_DIR}") ### Define library @@ -176,6 +187,7 @@ set(MEPHIT_LIBS hdf5::hdf5_fortran hdf5::hdf5_hl_fortran GSL::gsl + fortnum BLAS::BLAS LAPACK::LAPACK ${FFTW_LIBRARIES} diff --git a/src/mephit_pert.f90 b/src/mephit_pert.f90 index d50f372..4605bbb 100644 --- a/src/mephit_pert.f90 +++ b/src/mephit_pert.f90 @@ -1,7 +1,6 @@ module mephit_pert use iso_fortran_env, only: dp => real64 - use iso_c_binding, only: c_int, c_double implicit none @@ -37,17 +36,6 @@ function vector_element_projection(ktri, weight, R, Z, n, f) end function vector_element_projection end interface - interface - function gsl_sf_bessel_icn_array(nmin, nmax, x, result_array) & - bind(C, name='gsl_sf_bessel_In_array') - import :: c_int, c_double - integer(c_int), intent(in), value :: nmin, nmax - real(c_double), intent(in), value :: x - real(c_double), intent(inout), dimension(nmax - nmin + 1) :: result_array - integer(c_int) :: gsl_sf_bessel_icn_array - end function gsl_sf_bessel_icn_array - end interface - type :: L1_t !> Number of points on which the L1 are defined integer :: npoint @@ -1427,21 +1415,21 @@ end subroutine compute_kilca_vac_coeff !> field subroutine kilca_vacuum_fourier(tor_mode, pol_mode, R_0, r, vac_coeff, & B_rad, B_pol, B_tor) - use mephit_conf, only: logger use mephit_util, only: imun + use fortnum_special, only: bessel_in_array integer, intent(in) :: tor_mode, pol_mode real(dp), intent(in) :: R_0, r complex(dp), intent(in) :: vac_coeff complex(dp), intent(out) :: B_rad, B_pol, B_tor - real(c_double) :: I_m(-1:1), k_z_r - integer(c_int) :: status + real(dp) :: I_n(0:abs(pol_mode) + 1), I_m(-1:1), k_z_r k_z_r = tor_mode / R_0 * r - status = gsl_sf_bessel_icn_array(abs(pol_mode)-1, abs(pol_mode)+1, k_z_r, I_m) - if (status /= 0 .and. logger%err) then - write (logger%msg, '("gsl_sf_bessel_In_array returned error ", i0)') status - call logger%write_msg - end if + ! fortnum fills I_n(0:nmax) with I_0..I_nmax; recover the order |pol_mode|-1 + ! (which is -1 for pol_mode = 0) via the symmetry I_{-n}(x) = I_n(x). + call bessel_in_array(abs(pol_mode) + 1, k_z_r, I_n) + I_m(-1) = I_n(abs(abs(pol_mode) - 1)) + I_m(0) = I_n(abs(pol_mode)) + I_m(1) = I_n(abs(pol_mode) + 1) B_rad = 0.5d0 * (I_m(-1) + I_m(1)) * vac_coeff B_pol = imun * pol_mode / k_z_r * I_m(0) * vac_coeff B_tor = imun * I_m(0) * vac_coeff From 8f3e113d195a699b7751a7d986692f2eed0be124 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sun, 14 Jun 2026 11:32:42 +0200 Subject: [PATCH 2/7] Fetch fortnum over https so CI can clone the public repo --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index cb7835f..976ed51 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -105,7 +105,7 @@ find_or_fetch(libneo) include(FetchContent) if(NOT TARGET fortnum) FetchContent_Declare(fortnum - GIT_REPOSITORY git@github.com:lazy-fortran/fortnum.git + GIT_REPOSITORY https://github.com/lazy-fortran/fortnum.git GIT_TAG a7faa3c ) FetchContent_MakeAvailable(fortnum) From 54232cd762a67a4ef7bd25c749f16c5e76034302 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sun, 14 Jun 2026 16:29:48 +0200 Subject: [PATCH 3/7] Pin fortnum to the v0.1.0 release tag --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 976ed51..f334069 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -106,7 +106,7 @@ include(FetchContent) if(NOT TARGET fortnum) FetchContent_Declare(fortnum GIT_REPOSITORY https://github.com/lazy-fortran/fortnum.git - GIT_TAG a7faa3c + GIT_TAG v0.1.0 ) FetchContent_MakeAvailable(fortnum) endif() From 592554461c912af0d8a6f43b712eab8aa976a72a Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sun, 14 Jun 2026 09:43:44 +0200 Subject: [PATCH 4/7] Replace GSL quadrature and Levin sum in hyper1F1 with fortnum Route the a = 1 confluent hypergeometric evaluation through the fortnum C ABI: the boundary integral uses fortnum_integrate_qag with the Gauss-Kronrod 21 rule, and the Kummer-series acceleration uses fortnum_levin_u_accel. Both fortnum routines take the integrand plus an opaque context, so the existing quad_params struct is forwarded unchanged and no workspace is allocated. Add float.h for DBL_EPSILON, which was previously pulled in transitively through a GSL header. --- src/hyper1F1.c | 37 ++++++++++--------------------------- 1 file changed, 10 insertions(+), 27 deletions(-) diff --git a/src/hyper1F1.c b/src/hyper1F1.c index a4d2791..f59e901 100644 --- a/src/hyper1F1.c +++ b/src/hyper1F1.c @@ -3,13 +3,12 @@ Both Kummer series and continued fractions are used. */ +#include #include #include #include -#include -#include -#include +#include "fortnum.h" #include "hyper1F1.h" @@ -51,30 +50,18 @@ void hypergeometric1f1_quad(double *b_re, double *b_im, double *f_re, double *f_im) { // computes function 1F1(a,b,z) for a = 1 and complex b & z by quadrature - // must be optimized: avoid memory allocation! - - gsl_set_error_handler_off(); complex_double b = CMPLX(*b_re, *b_im), z = CMPLX(*z_re, *z_im); - size_t limit = 100; double epsabs = 1.0e-12, epsrel = 1.0e-12, err; - gsl_integration_workspace *w = gsl_integration_workspace_alloc(limit); - struct quad_params qp = {b, z, 0}; - gsl_function F; - F.function = &exp1mt; - F.params = &qp; - qp.part = 0; - gsl_integration_qag(&F, 0.0, 1.0, epsabs, epsrel, limit, GSL_INTEG_GAUSS21, w, f_re, &err); + fortnum_integrate_qag(&exp1mt, 0.0, 1.0, epsabs, epsrel, 21, f_re, &err, &qp); qp.part = 1; - gsl_integration_qag(&F, 0.0, 1.0, epsabs, epsrel, limit, GSL_INTEG_GAUSS21, w, f_im, &err); - - gsl_integration_workspace_free(w); + fortnum_integrate_qag(&exp1mt, 0.0, 1.0, epsabs, epsrel, 21, f_im, &err, &qp); } /*******************************************************************/ @@ -305,23 +292,19 @@ void hypergeometric1f1_kummer_modified_0_accel(double *b_re, double *b_im, t_im[n] = cimag(term[n]); } - gsl_sum_levin_u_workspace *w = gsl_sum_levin_u_alloc((size_t) N); - double err; - gsl_sum_levin_u_accel(t_re, (size_t) N, w, f_re, &err); + fortnum_levin_u_accel(t_re, N, f_re, &err); if (err > 1.0e-16) { - fprintf(stdout, "\nerr = %.16le sum_re = %.16le using %ld terms", - err, *f_re, w->terms_used); + fprintf(stdout, "\nerr = %.16le sum_re = %.16le using %d terms", + err, *f_re, N); } - gsl_sum_levin_u_accel(t_im, (size_t) N, w, f_im, &err); + fortnum_levin_u_accel(t_im, N, f_im, &err); if (err > 1.0e-16) { - fprintf(stdout, "\nerr = %.16le sum_im = %.16le using %ld terms", - err, *f_im, w->terms_used); + fprintf(stdout, "\nerr = %.16le sum_im = %.16le using %d terms", + err, *f_im, N); } - - gsl_sum_levin_u_free(w); } /*******************************************************************/ From b9da62768df3cc370e5afdd4574648bcb724264a Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sun, 14 Jun 2026 09:45:10 +0200 Subject: [PATCH 5/7] Replace GSL fixed-rule quadrature and error handler with fortnum Move the remaining C-side GSL call sites onto the fortnum C ABI. The Gauss-Legendre nodes and weights on the unit interval now come from fortnum_gauss_legendre_ab, which folds the table allocation and the [a, b] mapping into one call. Drop gsl_errno_msg and the gsl_set_error_handler registration: fortnum routines report failures through return codes rather than a global GSL error handler, so the dedicated GSL message hook and its forward declaration are no longer reachable. After this change libmephit.so resolves no GSL symbols. --- src/mephit_fem.c | 13 ++----------- src/mephit_run.c | 2 -- src/mephit_util.c | 11 ----------- src/mephit_util.h | 1 - 4 files changed, 2 insertions(+), 25 deletions(-) diff --git a/src/mephit_fem.c b/src/mephit_fem.c index e6829db..ff14866 100644 --- a/src/mephit_fem.c +++ b/src/mephit_fem.c @@ -6,8 +6,7 @@ #include #include #include -#include -#include +#include "fortnum.h" #include "triangle.h" #include "mephit_util.h" #include "mephit_fem.h" @@ -232,15 +231,7 @@ void FEM_deinit(void) void gauss_legendre_unit_interval(int order, double *points, double *weights) { - gsl_integration_glfixed_table *table; - size_t i, n; - - n = (size_t) order; - table = gsl_integration_glfixed_table_alloc(n); - for (i = 0; i < n; ++i) { - gsl_integration_glfixed_point(0.0, 1.0, i, &points[i], &weights[i], table); - } - gsl_integration_glfixed_table_free(table); + fortnum_gauss_legendre_ab(order, 0.0, 1.0, points, weights); } void FEM_triangulate_external(const int npt_inner, diff --git a/src/mephit_run.c b/src/mephit_run.c index 327becf..e2e7255 100644 --- a/src/mephit_run.c +++ b/src/mephit_run.c @@ -8,7 +8,6 @@ #include #include #include -#include #include "mephit_util.h" static volatile sig_atomic_t caught_signal = 0; @@ -131,7 +130,6 @@ int main(int argc, char *argv[]) } mephit_fork: mephit.pid = fork(); if (mephit.pid == (pid_t) 0) { - gsl_set_error_handler(gsl_errno_msg); mephit_run(runmode, config, suffix); exit(0); } else if (mephit.pid == (pid_t) -1) { diff --git a/src/mephit_util.c b/src/mephit_util.c index e7eb9aa..c4235a1 100644 --- a/src/mephit_util.c +++ b/src/mephit_util.c @@ -3,7 +3,6 @@ #include #include #include -#include #include "mephit_util.h" void timestamp(char *buffer) { @@ -45,13 +44,3 @@ void errno_msg(void (*exit_func)(int), const char *file, int line, int errnum, c exit_func(1); } } - -void gsl_errno_msg(const char *reason, const char *file, int line, int gsl_errno) -{ - char now[72] = "1990-06-11 20:47:00"; - - timestamp(now); - fprintf(stderr, "[%s] %s:%i:%s: %s.\n", - now, file, line, gsl_strerror(gsl_errno), reason); - _exit(1); -} diff --git a/src/mephit_util.h b/src/mephit_util.h index 518b08e..7aa2559 100644 --- a/src/mephit_util.h +++ b/src/mephit_util.h @@ -15,7 +15,6 @@ extern "C" { void timestamp(char *buffer); void errno_msg(void (*exit_func)(int), const char *file, int line, int errnum, const char *msg_fmt, ...); -void gsl_errno_msg(const char *reason, const char *file, int line, int gsl_errno); #ifdef __cplusplus } From 83ac593982bd821db16a43f121c81b529a956ed5 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Sun, 14 Jun 2026 09:52:26 +0200 Subject: [PATCH 6/7] Replace FFTW with fortnum fft_c2c and drop the FFTW dependency Move the fft_t one-dimensional forward transform onto fortnum's fft_c2c, which matches the FFTW sign convention (sign = -1, unnormalized). The type now owns an allocatable samples buffer instead of FFTW-allocated pointers and an opaque plan; apply transforms a copy in place so the caller's samples survive, then extracts and normalizes the requested mode range exactly as before. Remove src/fftw3.f90, the cmake/FindFFTW.cmake module, and the empty FFTW_LIBRARIES link entry. MEPHIT no longer references any FFTW symbol; libmephit.so and the executables link without FFTW. --- CMakeLists.txt | 2 - cmake/FindFFTW.cmake | 419 ------------------------------------------- src/fftw3.f90 | 4 - src/mephit_util.f90 | 54 ++---- 4 files changed, 14 insertions(+), 465 deletions(-) delete mode 100644 cmake/FindFFTW.cmake delete mode 100644 src/fftw3.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index f334069..ad63636 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -128,7 +128,6 @@ set_source_files_properties(src/external/umf4_f77zwrapper.c PROPERTIES COMPILE_F set(COMMON_FORTRAN_SRC_FILES src/external/netlib_mod.f90 src/external/d1mach.f - src/fftw3.f90 src/field_line_integration_for_SYNCH.f90 src/magdata_in_symfluxcoord.f90 src/points_2d.f90 @@ -190,7 +189,6 @@ set(MEPHIT_LIBS fortnum BLAS::BLAS LAPACK::LAPACK - ${FFTW_LIBRARIES} ${triangle_lib} LIBNEO::magfie LIBNEO::neo diff --git a/cmake/FindFFTW.cmake b/cmake/FindFFTW.cmake deleted file mode 100644 index ec1981e..0000000 --- a/cmake/FindFFTW.cmake +++ /dev/null @@ -1,419 +0,0 @@ -# - Find the FFTW library -# -# Original version of this file: -# Copyright (c) 2015, Wenzel Jakob -# https://github.com/wjakob/layerlab/blob/master/cmake/FindFFTW.cmake, commit 4d58bfdc28891b4f9373dfe46239dda5a0b561c6 -# Modifications: -# Copyright (c) 2017, Patrick Bos -# -# Usage: -# find_package(FFTW [REQUIRED] [QUIET] [COMPONENTS component1 ... componentX] ) -# -# It sets the following variables: -# FFTW_FOUND ... true if fftw is found on the system -# FFTW_[component]_LIB_FOUND ... true if the component is found on the system (see components below) -# FFTW_LIBRARIES ... full paths to all found fftw libraries -# FFTW_[component]_LIB ... full path to one of the components (see below) -# FFTW_INCLUDE_DIRS ... fftw include directory paths -# -# The following variables will be checked by the function -# FFTW_USE_STATIC_LIBS ... if true, only static libraries are found, otherwise both static and shared. -# FFTW_ROOT ... if set, the libraries are exclusively searched -# under this path -# -# This package supports the following components: -# FLOAT_LIB -# DOUBLE_LIB -# LONGDOUBLE_LIB -# FLOAT_THREADS_LIB -# DOUBLE_THREADS_LIB -# LONGDOUBLE_THREADS_LIB -# FLOAT_OPENMP_LIB -# DOUBLE_OPENMP_LIB -# LONGDOUBLE_OPENMP_LIB -# - -# TODO (maybe): extend with ExternalProject download + build option -# TODO: put on conda-forge - - -if( NOT FFTW_ROOT AND DEFINED ENV{FFTWDIR} ) - set( FFTW_ROOT $ENV{FFTWDIR} ) -endif() - -# Check if we can use PkgConfig -find_package(PkgConfig) - -#Determine from PKG -if( PKG_CONFIG_FOUND AND NOT FFTW_ROOT ) - pkg_check_modules( PKG_FFTW QUIET "fftw3" ) -endif() - -#Check whether to search static or dynamic libs -set( CMAKE_FIND_LIBRARY_SUFFIXES_SAV ${CMAKE_FIND_LIBRARY_SUFFIXES} ) - -if( ${FFTW_USE_STATIC_LIBS} ) - set( CMAKE_FIND_LIBRARY_SUFFIXES ${CMAKE_STATIC_LIBRARY_SUFFIX} ) -else() - set( CMAKE_FIND_LIBRARY_SUFFIXES ${CMAKE_FIND_LIBRARY_SUFFIXES_SAV} ) -endif() - -if( FFTW_ROOT ) - # find libs - - find_library( - FFTW_DOUBLE_LIB - NAMES "fftw3" libfftw3-3 - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTW_DOUBLE_THREADS_LIB - NAMES "fftw3_threads" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTW_DOUBLE_OPENMP_LIB - NAMES "fftw3_omp" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTW_DOUBLE_MPI_LIB - NAMES "fftw3_mpi" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTW_FLOAT_LIB - NAMES "fftw3f" libfftw3f-3 - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTW_FLOAT_THREADS_LIB - NAMES "fftw3f_threads" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTW_FLOAT_OPENMP_LIB - NAMES "fftw3f_omp" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTW_FLOAT_MPI_LIB - NAMES "fftw3f_mpi" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTW_LONGDOUBLE_LIB - NAMES "fftw3l" libfftw3l-3 - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTW_LONGDOUBLE_THREADS_LIB - NAMES "fftw3l_threads" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTW_LONGDOUBLE_OPENMP_LIB - NAMES "fftw3l_omp" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - find_library( - FFTW_LONGDOUBLE_MPI_LIB - NAMES "fftw3l_mpi" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "lib" "lib64" - NO_DEFAULT_PATH - ) - - #find includes - find_path(FFTW_INCLUDE_DIRS - NAMES "fftw3.h" - PATHS ${FFTW_ROOT} - PATH_SUFFIXES "include" - NO_DEFAULT_PATH - ) - -else() - - find_library( - FFTW_DOUBLE_LIB - NAMES "fftw3" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_library( - FFTW_DOUBLE_THREADS_LIB - NAMES "fftw3_threads" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_library( - FFTW_DOUBLE_OPENMP_LIB - NAMES "fftw3_omp" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_library( - FFTW_DOUBLE_MPI_LIB - NAMES "fftw3_mpi" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_library( - FFTW_FLOAT_LIB - NAMES "fftw3f" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_library( - FFTW_FLOAT_THREADS_LIB - NAMES "fftw3f_threads" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_library( - FFTW_FLOAT_OPENMP_LIB - NAMES "fftw3f_omp" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_library( - FFTW_FLOAT_MPI_LIB - NAMES "fftw3f_mpi" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_library( - FFTW_LONGDOUBLE_LIB - NAMES "fftw3l" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_library( - FFTW_LONGDOUBLE_THREADS_LIB - NAMES "fftw3l_threads" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_library(FFTW_LONGDOUBLE_OPENMP_LIB - NAMES "fftw3l_omp" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_library(FFTW_LONGDOUBLE_MPI_LIB - NAMES "fftw3l_mpi" - PATHS ${PKG_FFTW_LIBRARY_DIRS} ${LIB_INSTALL_DIR} - ) - - find_path(FFTW_INCLUDE_DIRS - NAMES "fftw3.h" - PATHS ${PKG_FFTW_INCLUDE_DIRS} ${INCLUDE_INSTALL_DIR} - ) - -endif( FFTW_ROOT ) - -#--------------------------------------- components - -if (FFTW_DOUBLE_LIB) - set(FFTW_DOUBLE_LIB_FOUND TRUE) - set(FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTW_DOUBLE_LIB}) - add_library(FFTW::Double INTERFACE IMPORTED) - set_target_properties(FFTW::Double - PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${FFTW_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${FFTW_DOUBLE_LIB}" - ) -else() - set(FFTW_DOUBLE_LIB_FOUND FALSE) -endif() - -if (FFTW_FLOAT_LIB) - set(FFTW_FLOAT_LIB_FOUND TRUE) - set(FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTW_FLOAT_LIB}) - add_library(FFTW::Float INTERFACE IMPORTED) - set_target_properties(FFTW::Float - PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${FFTW_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${FFTW_FLOAT_LIB}" - ) -else() - set(FFTW_FLOAT_LIB_FOUND FALSE) -endif() - -if (FFTW_LONGDOUBLE_LIB) - set(FFTW_LONGDOUBLE_LIB_FOUND TRUE) - set(FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTW_LONGDOUBLE_LIB}) - add_library(FFTW::LongDouble INTERFACE IMPORTED) - set_target_properties(FFTW::LongDouble - PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${FFTW_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${FFTW_LONGDOUBLE_LIB}" - ) -else() - set(FFTW_LONGDOUBLE_LIB_FOUND FALSE) -endif() - -if (FFTW_DOUBLE_THREADS_LIB) - set(FFTW_DOUBLE_THREADS_LIB_FOUND TRUE) - set(FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTW_DOUBLE_THREADS_LIB}) - add_library(FFTW::DoubleThreads INTERFACE IMPORTED) - set_target_properties(FFTW::DoubleThreads - PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${FFTW_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${FFTW_DOUBLE_THREADS_LIB}" - ) -else() - set(FFTW_DOUBLE_THREADS_LIB_FOUND FALSE) -endif() - -if (FFTW_FLOAT_THREADS_LIB) - set(FFTW_FLOAT_THREADS_LIB_FOUND TRUE) - set(FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTW_FLOAT_THREADS_LIB}) - add_library(FFTW::FloatThreads INTERFACE IMPORTED) - set_target_properties(FFTW::FloatThreads - PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${FFTW_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${FFTW_FLOAT_THREADS_LIB}" - ) -else() - set(FFTW_FLOAT_THREADS_LIB_FOUND FALSE) -endif() - -if (FFTW_LONGDOUBLE_THREADS_LIB) - set(FFTW_LONGDOUBLE_THREADS_LIB_FOUND TRUE) - set(FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTW_LONGDOUBLE_THREADS_LIB}) - add_library(FFTW::LongDoubleThreads INTERFACE IMPORTED) - set_target_properties(FFTW::LongDoubleThreads - PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${FFTW_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${FFTW_LONGDOUBLE_THREADS_LIB}" - ) -else() - set(FFTW_LONGDOUBLE_THREADS_LIB_FOUND FALSE) -endif() - -if (FFTW_DOUBLE_OPENMP_LIB) - set(FFTW_DOUBLE_OPENMP_LIB_FOUND TRUE) - set(FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTW_DOUBLE_OPENMP_LIB}) - add_library(FFTW::DoubleOpenMP INTERFACE IMPORTED) - set_target_properties(FFTW::DoubleOpenMP - PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${FFTW_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${FFTW_DOUBLE_OPENMP_LIB}" - ) -else() - set(FFTW_DOUBLE_OPENMP_LIB_FOUND FALSE) -endif() - -if (FFTW_FLOAT_OPENMP_LIB) - set(FFTW_FLOAT_OPENMP_LIB_FOUND TRUE) - set(FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTW_FLOAT_OPENMP_LIB}) - add_library(FFTW::FloatOpenMP INTERFACE IMPORTED) - set_target_properties(FFTW::FloatOpenMP - PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${FFTW_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${FFTW_FLOAT_OPENMP_LIB}" - ) -else() - set(FFTW_FLOAT_OPENMP_LIB_FOUND FALSE) -endif() - -if (FFTW_LONGDOUBLE_OPENMP_LIB) - set(FFTW_LONGDOUBLE_OPENMP_LIB_FOUND TRUE) - set(FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTW_LONGDOUBLE_OPENMP_LIB}) - add_library(FFTW::LongDoubleOpenMP INTERFACE IMPORTED) - set_target_properties(FFTW::LongDoubleOpenMP - PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${FFTW_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${FFTW_LONGDOUBLE_OPENMP_LIB}" - ) -else() - set(FFTW_LONGDOUBLE_OPENMP_LIB_FOUND FALSE) -endif() - -if (FFTW_DOUBLE_MPI_LIB) - set(FFTW_DOUBLE_MPI_LIB_FOUND TRUE) - set(FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTW_DOUBLE_MPI_LIB}) - add_library(FFTW::DoubleMPI INTERFACE IMPORTED) - set_target_properties(FFTW::DoubleMPI - PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${FFTW_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${FFTW_DOUBLE_MPI_LIB}" - ) -else() - set(FFTW_DOUBLE_MPI_LIB_FOUND FALSE) -endif() - -if (FFTW_FLOAT_MPI_LIB) - set(FFTW_FLOAT_MPI_LIB_FOUND TRUE) - set(FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTW_FLOAT_MPI_LIB}) - add_library(FFTW::FloatMPI INTERFACE IMPORTED) - set_target_properties(FFTW::FloatMPI - PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${FFTW_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${FFTW_FLOAT_MPI_LIB}" - ) -else() - set(FFTW_FLOAT_MPI_LIB_FOUND FALSE) -endif() - -if (FFTW_LONGDOUBLE_MPI_LIB) - set(FFTW_LONGDOUBLE_MPI_LIB_FOUND TRUE) - set(FFTW_LIBRARIES ${FFTW_LIBRARIES} ${FFTW_LONGDOUBLE_MPI_LIB}) - add_library(FFTW::LongDoubleMPI INTERFACE IMPORTED) - set_target_properties(FFTW::LongDoubleMPI - PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${FFTW_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${FFTW_LONGDOUBLE_MPI_LIB}" - ) -else() - set(FFTW_LONGDOUBLE_MPI_LIB_FOUND FALSE) -endif() - -#--------------------------------------- end components - -set( CMAKE_FIND_LIBRARY_SUFFIXES ${CMAKE_FIND_LIBRARY_SUFFIXES_SAV} ) - -include(FindPackageHandleStandardArgs) - -find_package_handle_standard_args(FFTW - REQUIRED_VARS FFTW_INCLUDE_DIRS - HANDLE_COMPONENTS - ) - -mark_as_advanced( - FFTW_INCLUDE_DIRS - FFTW_LIBRARIES - FFTW_FLOAT_LIB - FFTW_DOUBLE_LIB - FFTW_LONGDOUBLE_LIB - FFTW_FLOAT_THREADS_LIB - FFTW_DOUBLE_THREADS_LIB - FFTW_LONGDOUBLE_THREADS_LIB - FFTW_FLOAT_OPENMP_LIB - FFTW_DOUBLE_OPENMP_LIB - FFTW_LONGDOUBLE_OPENMP_LIB - FFTW_FLOAT_MPI_LIB - FFTW_DOUBLE_MPI_LIB - FFTW_LONGDOUBLE_MPI_LIB - ) diff --git a/src/fftw3.f90 b/src/fftw3.f90 deleted file mode 100644 index f9e5f78..0000000 --- a/src/fftw3.f90 +++ /dev/null @@ -1,4 +0,0 @@ -module FFTW3 - use, intrinsic :: iso_c_binding - include 'fftw3.f03' -end module FFTW3 diff --git a/src/mephit_util.f90 b/src/mephit_util.f90 index 0bc1e8a..fefecca 100644 --- a/src/mephit_util.f90 +++ b/src/mephit_util.f90 @@ -1,7 +1,6 @@ module mephit_util use iso_fortran_env, only: dp => real64 - use iso_c_binding, only: c_ptr, c_null_ptr, c_double_complex implicit none @@ -36,9 +35,7 @@ module mephit_util type :: fft_t integer :: N = 0 - complex(c_double_complex), dimension(:), pointer :: samples => null() - complex(c_double_complex), dimension(:), pointer, private :: modes => null() - type(c_ptr), private :: plan_p = c_null_ptr, samples_p = c_null_ptr, modes_p = c_null_ptr + complex(dp), dimension(:), allocatable :: samples contains procedure :: init => fft_init procedure :: apply => fft_apply @@ -375,29 +372,21 @@ subroutine binsearch(x, lb, xi, k) end subroutine binsearch subroutine fft_init(fft, N) - use iso_c_binding, only: c_size_t, c_f_pointer - use fftw3, only: fftw_alloc_complex, fftw_plan_dft_1d, & - FFTW_FORWARD, FFTW_PATIENT, FFTW_DESTROY_INPUT class(fft_t), intent(inout) :: fft integer, intent(in) :: N call fft_deinit(fft) fft%N = N - fft%samples_p = fftw_alloc_complex(int(N, c_size_t)) - call c_f_pointer(fft%samples_p, fft%samples, [N]) - fft%modes_p = fftw_alloc_complex(int(N, c_size_t)) - call c_f_pointer(fft%modes_p, fft%modes, [N]) - fft%plan_p = fftw_plan_dft_1d(N, fft%samples, fft%modes, & - FFTW_FORWARD, ior(FFTW_PATIENT, FFTW_DESTROY_INPUT)) + allocate(fft%samples(N)) end subroutine fft_init subroutine fft_apply(fft, m_min, m_max, modes) - use iso_c_binding, only: c_associated - use fftw3, only: fftw_execute_dft use mephit_conf, only: logger + use fortnum_fft, only: fft_c2c class(fft_t), intent(inout) :: fft integer, intent(in) :: m_min, m_max complex(dp), dimension(m_min:), intent(inout) :: modes + complex(dp), dimension(:), allocatable :: spectrum if (ubound(modes, 1) /= m_max) then call logger%msg_arg_size('fft_apply', & @@ -406,44 +395,29 @@ subroutine fft_apply(fft, m_min, m_max, modes) if (logger%err) call logger%write_msg error stop end if - if (.not. (c_associated(fft%plan_p) .and. & - c_associated(fft%modes_p) .and. & - c_associated(fft%samples_p))) then - logger%msg = 'Attempt to dereference null pointer in fft_apply' + if (.not. allocated(fft%samples)) then + logger%msg = 'Attempt to use uninitialized samples in fft_apply' if (logger%err) call logger%write_msg error stop end if - call fftw_execute_dft(fft%plan_p, fft%samples, fft%modes) + ! fft_c2c is in-place; transform a copy so the caller's samples survive. + spectrum = fft%samples + call fft_c2c(spectrum, -1) if (m_min >= 0 .and. m_max >= 0) then - modes(:) = fft%modes(m_min+1:m_max+1) / dble(fft%N) + modes(:) = spectrum(m_min+1:m_max+1) / dble(fft%N) else if (m_min < 0 .and. m_max < 0) then - modes(:) = fft%modes(fft%N+m_min+1:fft%N+m_max+1) / dble(fft%N) + modes(:) = spectrum(fft%N+m_min+1:fft%N+m_max+1) / dble(fft%N) else - modes(m_min:-1) = fft%modes(fft%N+m_min+1:fft%N) / dble(fft%N) - modes(0:m_max) = fft%modes(:m_max+1) / dble(fft%N) + modes(m_min:-1) = spectrum(fft%N+m_min+1:fft%N) / dble(fft%N) + modes(0:m_max) = spectrum(:m_max+1) / dble(fft%N) end if end subroutine fft_apply subroutine fft_deinit(fft) - use iso_c_binding, only: c_associated, c_null_ptr - use fftw3, only: fftw_destroy_plan, fftw_free class(fft_t), intent(inout) :: fft fft%N = 0 - fft%samples => null() - fft%modes => null() - if (c_associated(fft%plan_p)) then - call fftw_destroy_plan(fft%plan_p) - fft%plan_p = c_null_ptr - end if - if (c_associated(fft%modes_p)) then - call fftw_free(fft%modes_p) - fft%modes_p = c_null_ptr - end if - if (c_associated(fft%samples_p)) then - call fftw_free(fft%samples_p) - fft%samples_p = c_null_ptr - end if + if (allocated(fft%samples)) deallocate(fft%samples) end subroutine fft_deinit !> Transform components of a vector \f$ \vec{v} \f$ from straight cylinder coordinates From 317d31292c5eedcf5b74e4768e3892f905d7c7ac Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 15 Jun 2026 20:11:47 +0200 Subject: [PATCH 7/7] Re-pin fortnum to current main after upstream history rewrite --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ad63636..2d94601 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -106,7 +106,7 @@ include(FetchContent) if(NOT TARGET fortnum) FetchContent_Declare(fortnum GIT_REPOSITORY https://github.com/lazy-fortran/fortnum.git - GIT_TAG v0.1.0 + GIT_TAG 92de6e949a772cfffc73bb5295fe5e2b056b9c18 ) FetchContent_MakeAvailable(fortnum) endif()