Migrate off GSL onto fortnum

CMakeLists.txt

@@ -59,7 +59,6 @@ endif()
 cmake_policy(SET CMP0074 NEW)
 find_package(HDF5 COMPONENTS C Fortran HL REQUIRED)
 find_package(Boost 1.74 REQUIRED)
-find_package(GSL REQUIRED)
 
 ### NetCDF
 find_program(NF_CONFIG "nf-config")
@@ -100,6 +99,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 main
+  )
+  FetchContent_MakeAvailable(fortnum)
+endif()
+
 ### Project source files
 # SuiteSparse
 set(SUITESPARSE_SRC_FILES
@@ -152,7 +162,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
@@ -175,7 +185,7 @@ set(MEPHIT_LIBS
   hdf5::hdf5_hl
   hdf5::hdf5_fortran
   hdf5::hdf5_hl_fortran
-  GSL::gsl
+  fortnum
   BLAS::BLAS
   LAPACK::LAPACK
   ${FFTW_LIBRARIES}

cmake/SetupCODE.cmake

@@ -3,9 +3,6 @@ if (DEFINED ENV{CODE})
     if(EXISTS $ENV{CODE}/libneo/build/)
         set(LIBNEO_DIR $ENV{CODE}/libneo/build/ CACHE STRING "libneo directory")
     endif()
-    if(EXISTS $ENV{CODE}/external/fgsl-1.6.0/)
-        set(FGSL_DIR $ENV{CODE}/external/fgsl-1.6.0/ CACHE STRING "FGSL directory")
-    endif()
     if(EXISTS $ENV{CODE}/external/triangle/)
         set(TRIANGLE_DIR $ENV{CODE}/external/triangle/ CACHE STRING "TRIANGLE directory")
     endif()

src/hyper1F1.c

@@ -3,13 +3,12 @@
   Both Kummer series and continued fractions are used.
 */
 
+#include <float.h>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_integration.h>
-#include <gsl/gsl_sum.h>
+#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);
 }
 
 /*******************************************************************/

src/mephit_fem.c

@@ -6,8 +6,7 @@
 #include <fcntl.h>
 #include <sys/types.h>
 #include <sys/stat.h>
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_integration.h>
+#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,

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

src/mephit_run.c

@@ -8,7 +8,6 @@
 #include <sys/stat.h>
 #include <sys/wait.h>
 #include <signal.h>
-#include <gsl/gsl_errno.h>
 #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) {

src/mephit_util.c

@@ -3,7 +3,6 @@
 #include <string.h>
 #include <time.h>
 #include <unistd.h>
-#include <gsl/gsl_errno.h>
 #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);
-}

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
 }