From 15587f60d607a1dd311af5fb49050b70aa7b233d Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 15 Jun 2026 16:42:09 +0200 Subject: [PATCH] Replace GSL adaptive quadrature in KiLCA and QL-Balance with fortnum --- KiLCA/CMakeLists.txt | 5 ++++- .../conductivity/tests/calc_I_array_drift.cpp | 20 +++---------------- KiLCA/math/fourier/four_transf.cpp | 15 +++----------- QL-Balance/CMakeLists.txt | 4 ++++ QL-Balance/src/base/vel_integral.cpp | 18 +++++------------ 5 files changed, 19 insertions(+), 43 deletions(-) diff --git a/KiLCA/CMakeLists.txt b/KiLCA/CMakeLists.txt index 29c0307f..d8e63fa8 100644 --- a/KiLCA/CMakeLists.txt +++ b/KiLCA/CMakeLists.txt @@ -111,7 +111,7 @@ set (kilca_eigparam_main progs/main_eig_param.cpp) set (libs_root_path ${CMAKE_BINARY_DIR}/external-install) set (gsl_include_path ${libs_root_path}/gsl/include) -set(EXTERNAL_LIBS slatec lapack blas SUNDIALS::cvode SUNDIALS::nvecserial gsl gslcblas) +set(EXTERNAL_LIBS slatec lapack blas SUNDIALS::cvode SUNDIALS::nvecserial fortnum gsl gslcblas) # configure a header file to pass some of the CMake settings to the source code configure_file( @@ -122,6 +122,8 @@ configure_file( include_directories("${CMAKE_CURRENT_BINARY_DIR}") include_directories ("${gsl_include_path}") +# fortnum C ABI header (fortnum.h); the Fortran fortnum target is linked above. +include_directories ("${fortnum_SOURCE_DIR}/include") include_directories ("${PROJECT_BINARY_DIR}") include_directories ("${CMAKE_CURRENT_SOURCE_DIR}/antenna") include_directories ("${CMAKE_CURRENT_SOURCE_DIR}/background") @@ -283,6 +285,7 @@ target_include_directories(kilca_lib PUBLIC # External dependencies ${sundials_include_path} ${gsl_include_path} + ${fortnum_SOURCE_DIR}/include ) add_executable (kilca_normal_exe ${kilca_normal_main}) diff --git a/KiLCA/flre/conductivity/tests/calc_I_array_drift.cpp b/KiLCA/flre/conductivity/tests/calc_I_array_drift.cpp index 851bc0f4..f783c31f 100644 --- a/KiLCA/flre/conductivity/tests/calc_I_array_drift.cpp +++ b/KiLCA/flre/conductivity/tests/calc_I_array_drift.cpp @@ -5,8 +5,7 @@ #include #include -#include -#include +#include "fortnum.h" #include "constants.h" @@ -95,8 +94,6 @@ int calc_imnl_quad_ (int * m_max, int * n_max, int * l_max, double * x4_re, double * x4_im, double * Imnl) { -gsl_set_error_handler_off(); - complex x1(*x1_re, *x1_im); complex x2(*x2_re, *x2_im); complex x3(*x3_re, *x3_im); @@ -117,15 +114,8 @@ x4 = x4 / ep2a; quad_params qp = {x1, x2, x3, x4, 0, gamma, 0, 0, 0}; -size_t limit = 1024; double epsabs = 1.0e-14, epsrel = 1.0e-14, err; -gsl_integration_workspace * w = gsl_integration_workspace_alloc(limit); - -gsl_function F; -F.function = &subint; -F.params = &qp; - double x = 0.0; double y = 0.0; @@ -154,13 +144,11 @@ for(int l = 0; l <= *l_max; ++l) { double I_re; qp.part = 0; - //gsl_integration_qagiu(&F, 0.0e0, epsabs, epsrel, limit, w, &I_re, &err); - gsl_integration_qag(&F, t1, t2, epsabs, epsrel, limit, GSL_INTEG_GAUSS21, w, &I_re, &err); + fortnum_integrate_qag(&subint, t1, t2, epsabs, epsrel, 21, &I_re, &err, &qp); double I_im; qp.part = 1; - //gsl_integration_qagiu(&F, 0.0e0, epsabs, epsrel, limit, w, &I_im, &err); - gsl_integration_qag(&F, t1, t2, epsabs, epsrel, limit, GSL_INTEG_GAUSS21, w, &I_im, &err); + fortnum_integrate_qag(&subint, t1, t2, epsabs, epsrel, 21, &I_im, &err, &qp); InT = (I_re + I_im * I) / pow(ep2a, 0.5e0*(3 + m + n)); @@ -172,8 +160,6 @@ for(int l = 0; l <= *l_max; ++l) } } -gsl_integration_workspace_free(w); - return 0; } diff --git a/KiLCA/math/fourier/four_transf.cpp b/KiLCA/math/fourier/four_transf.cpp index f2d24eac..41b9e004 100644 --- a/KiLCA/math/fourier/four_transf.cpp +++ b/KiLCA/math/fourier/four_transf.cpp @@ -9,9 +9,7 @@ #include "spline.h" #include "constants.h" -#include -#include -#include +#include "fortnum.h" /*-----------------------------------------------------------------*/ @@ -72,14 +70,8 @@ spline_alloc_ (N, 1, dimx, x, C, &sid); spline_calc_ (sid, yt, 0, 1, NULL, &ierr); -gsl_integration_workspace *W = gsl_integration_workspace_alloc (1000); - struct quad_func_params P = {0, 0, sid, pi/delta, x0}; -gsl_function F; -F.function = &func; -F.params = &P; - double re, im, err_re, err_im; for (k=-M; k<=M; k++) @@ -87,10 +79,10 @@ for (k=-M; k<=M; k++) P.k = k; P.part = 0; - gsl_integration_qag (&F, x0-delta, x0+delta, 1.0e-8, 1.0e-8, 1000, 3, W, &re, &err_re); + fortnum_integrate_qag (&func, x0-delta, x0+delta, 1.0e-8, 1.0e-8, 31, &re, &err_re, &P); P.part = 1; - gsl_integration_qag (&F, x0-delta, x0+delta, 1.0e-8, 1.0e-8, 1000, 3, W, &im, &err_im); + fortnum_integrate_qag (&func, x0-delta, x0+delta, 1.0e-8, 1.0e-8, 31, &im, &err_im, &P); //fprintf (stdout, "\nk = %d:\tresult = %lg%+lgi\terror = %lg%+lgi", k, re, im, err_re, err_im); @@ -102,7 +94,6 @@ delete [] C; delete [] yt; spline_free_ (sid); -gsl_integration_workspace_free (W); return 0; } diff --git a/QL-Balance/CMakeLists.txt b/QL-Balance/CMakeLists.txt index 82690ceb..70a5d3d1 100644 --- a/QL-Balance/CMakeLists.txt +++ b/QL-Balance/CMakeLists.txt @@ -124,6 +124,9 @@ set_source_files_properties(${base_balance_src}/kim_wave_code_adapter.f90 # Include central slatec math sources via cache variable defined by common/math target_include_directories(ql-balance_lib PUBLIC ${SLATEC_INCLUDE_DIR}) +# fortnum C ABI header (fortnum.h) for vel_integral.cpp. +target_include_directories(ql-balance_lib PUBLIC ${fortnum_SOURCE_DIR}/include) + find_package(HDF5 REQUIRED COMPONENTS C HL Fortran Fortran_HL) message(" ====== HDF5 was found! ======") @@ -150,6 +153,7 @@ target_link_libraries(ql-balance_lib PUBLIC SUNDIALS::cvode SUNDIALS::nvecserial blas + fortnum gsl gslcblas hdf5::hdf5 diff --git a/QL-Balance/src/base/vel_integral.cpp b/QL-Balance/src/base/vel_integral.cpp index 1e354fb2..0e830c82 100644 --- a/QL-Balance/src/base/vel_integral.cpp +++ b/QL-Balance/src/base/vel_integral.cpp @@ -7,9 +7,7 @@ #include "constants.h" -#include -#include -#include +#include "fortnum.h" /*-----------------------------------------------------------------*/ @@ -65,14 +63,8 @@ return exp(-x*x)/(pow((P->a)*x + P->omE, 2)+(P->nu)*(P->nu))*real(field_fac)*ind void calc_velocity_integral_ (int ind, double vT, double ks, double kp, double omE, double nu, double *Es, double *Ep, double *res) { -gsl_integration_workspace *W = gsl_integration_workspace_alloc (10000); - struct quad_func_params P; -gsl_function F; -F.function = &func; -F.params = &P; - //set structure: P.ind = ind; P.a = sqrt(2.0)*vT*kp; @@ -100,10 +92,10 @@ for (i=0; i<10000; i++) double err; -gsl_integration_qagi (&F, 1.0e-6, 1.0e-6, 10000, W, res, &err); -//gsl_integration_qag (&F, -10.0, 10.0, 1.0e-6, 1.0e-6, 10000, 3, W, res, &err); - -gsl_integration_workspace_free (W); +// gsl_integration_qagi over (-inf, +inf): fortnum QAGIU with inf=2 splits at +// the bound and sums the two semi-infinite transforms (bound ignored except +// as the split point, matching QAGI's split at 0). +fortnum_integrate_qagiu (&func, 0.0, 2, 1.0e-6, 1.0e-6, res, &err, &P); } /*-----------------------------------------------------------------*/