Hello!
I am working with a fortran project which uses fgsl (v1.2.0, together
with gsl v2.3.0).
The code reports an arithmetic exception in an integration routine, when
trying to integrate exp(-x^2) from x=19.3 to 20.0.
I would expect the code to work without problems (maybe underflow or
denormal).
Trying to create a minimal working example i came up with the attached
file (the cmake input i used for compilation is also attached, sorry I
could not boil this further down).
As a corresponding version for c++ works (also attached), this might
relate to a bug in fgsl, but for me (a person without knowledge of the
internals), that seems questionable according to what causes the problem
and what not.
The output is as follows:
Program received signal SIGFPE: Floating-point exception - erroneous
arithmetic operation.
Backtrace for this error:
#0 0x7fa39acc5d1d in ???
#1 0x7fa39acc4f7d in ???
#2 0x7fa39a1e405f in ???
#3 0x560ac54c9eff in ???
#4 0x560ac54860cf in neo2
at
/afs/itp.tugraz.at/user/buchholz/Programs/neo2-minimal/NEO-2-QL/neo2.f90:31
#5 0x560ac54861d8 in main
at
/afs/itp.tugraz.at/user/buchholz/Programs/neo2-minimal/NEO-2-QL/neo2.f90:4
It does not occur with the alternative integration routine commented in
the code.
This does not occur, if one reduces the lower bound to 19.2 or smaller
(which is a reason why I suspect that this might be a gsl bug - a bug in
the interface should have no influence on this).
With kind regards
Rico Buchholz
PS: I did also send the bug report to fgsl.
### Initialize CMake (some internal stuff)
cmake_minimum_required(VERSION 3.0)
set(CMAKE_DISABLE_IN_SOURCE_BUILD ON)
# Enable project
project(NEO-2-NTV)
enable_language(Fortran)
set(EXE_NAME "neo_2.x")
### Define paths to external libraries (load external file)
### Installation directory
set (CMAKE_INSTALL_PREFIX ${CMAKE_SOURCE_DIR})
### Path to the external libraries
set(PROJLIBS /proj/plasma/Libs/ CACHE STRING "Common library path")
set(GSL_PATH ${PROJLIBS}GSL_gcc-6.3/gsl-2.3/LIB)
set(GSL_LIB ${GSL_PATH}/lib/ CACHE STRING "GSL path")
set(FGSL_PATH /proj/plasma/Libs/FGSL_gcc-6.3/fgsl-1.2.0/LIB CACHE STRING "FGSL
path")
set(FGSL_LIB ${FGSL_PATH}/lib/)
set(FGSL_INC ${FGSL_PATH}/include/fgsl/)
### Path to the ATLAS libraries
set(ATLASLIBS /usr/lib/atlas-base/)
set(LAPACK_LIB ${ATLASLIBS}/atlas/)
set(F77BLAS_LIB ${ATLASLIBS})
set(ATLASBLAS_LIB ${ATLASLIBS})
### Path to OpenBLAS
set(OPEN_BLAS_PATH /usr/lib/openblas-base/)
set(OPEN_BLAS_LAPACK_LIB ${OPEN_BLAS_PATH})
set(OPEN_BLAS_LIB ${OPEN_BLAS_PATH})
# Print the current compiler
message(STATUS "The Compiler ID is ${CMAKE_Fortran_COMPILER_ID}")
# Include directories for NEO-2
include_directories(${NEO2_Inc})
set(CMAKE_Fortran_FLAGS_DEBUG "-O0 -cpp -g -ffpe-trap=invalid,zero,overflow
-fcheck=all -fbacktrace -Wall -Wno-unused-dummy-argument")
set(CFLAGS "-O2 -DDOUBLE_APPEND_FORTRAN")
### Define container for *.o and *.mod files
set(CMAKE_Fortran_MODULE_DIRECTORY ${PROJECT_BINARY_DIR}/OBJS)
### Find libraries
# GSL
find_library(gsl_lib gsl ${GSL_LIB})
# FGSL
find_library(fgsl_lib fgsl ${FGSL_LIB} NO_DEFAULT_PATH)
include_directories(${FGSL_INC})
# LAPACK and BLAS
if (CMAKE_Fortran_COMPILER_ID MATCHES "GNU")
find_library(lapack_lib lapack ${LAPACK_LIB} NO_DEFAULT_PATH)
find_library(f77blas_lib f77blas ${F77BLAS_LIB} NO_DEFAULT_PATH)
find_library(atlas_lib atlas ${ATLASBLAS_LIB} NO_DEFAULT_PATH)
find_library(open_blas_lapack_lib lapack ${OPEN_BLAS_LAPACK_LIB}
NO_DEFAULT_PATH)
find_library(open_blas_lib blas ${OPEN_BLAS_LIB} NO_DEFAULT_PATH)
endif ()
### Source files
set(NEO2_SRC_FILES
./neo2.f90
)
set_source_files_properties(${NEO2_SRC_FILES} PROPERTIES COMPILE_FLAGS
"${DEBUGFLAG} ${FFLAG_DEBUG} ${FFLAGS}")
### Define executable
add_executable(${EXE_NAME}
${NEO2_SRC_FILES}
)
if (CMAKE_Fortran_COMPILER_ID MATCHES "GNU")
set(LINALG_LIBRARIES ${lapack_lib} ${f77blas_lib} ${atlas_lib})
endif ()
### Add libraries linked with executable
target_link_libraries(${EXE_NAME}
${fgsl_lib} ${gsl_lib}
${LINALG_LIBRARIES}
)
program neo2
use, intrinsic :: iso_c_binding
use fgsl ! Fortran interface of the GSL Library
implicit none
real(fgsl_double), parameter :: x_low = 1.93d1
real(fgsl_double), parameter :: x_up = 2.0d1
! absolute error
real(fgsl_double), parameter :: epsabs = 1.0d-13
! relative error
real(fgsl_double), parameter :: epsrel = 1.0d-13
integer(fgsl_size_t), parameter :: n = 100_fgsl_size_t
integer(fgsl_size_t) :: neval ! number of function evaluations
real(fgsl_double) :: ra, rda ! result and absolute error
integer(fgsl_int) :: statusval
type(fgsl_function) :: stdfunc
type(fgsl_integration_cquad_workspace) :: workspace
!~ type(fgsl_integration_workspace) :: workspace
neval = 0
stdfunc = fgsl_function_init(func_to_integrate, c_null_ptr)
workspace = fgsl_integration_cquad_workspace_alloc(n)
!~ workspace = fgsl_integration_workspace_alloc(1000_fgsl_size_t)
statusval = fgsl_integration_cquad(stdfunc, x_low, x_up, &
& epsabs, epsrel, workspace, ra, rda, neval)
!~ statusval = fgsl_integration_qags(stdfunc, x_low, x_up, &
!~ & epsabs, epsrel, 1000_fgsl_size_t, workspace, ra, rda)
call fgsl_function_free(stdfunc)
call fgsl_integration_cquad_workspace_free(workspace)
!~ call fgsl_integration_workspace_free(workspace)
write(*,*) ra, rda, neval
stop
contains
function func_to_integrate(x, params) bind(c)
real(c_double), value :: x
type(c_ptr), value :: params
real(c_double) :: func_to_integrate
func_to_integrate = exp(-x**2)
end function func_to_integrate
end program neo2
/**
*
* g++ main.cpp -o main.x -L`gsl-config --libs`
*
*/
#include <iostream>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_integration.h>
double f(double x, void * params) {
return std::exp(-x*x);
}
int main (void)
{
size_t const n(100);
double const a(19.3);
double const b(20.0);
double const epsabs(1.0e-13);
double const epsrel(1.0e-13);
double result(0.0);
double abserr(0.0);
size_t nevals(0);
gsl_integration_cquad_workspace * workspace(NULL);
gsl_function F;
F.function = &f;
F.params = NULL;
workspace = gsl_integration_cquad_workspace_alloc(n);
gsl_integration_cquad(&F, a, b, epsabs, epsrel, workspace, &result, &abserr, &nevals);
std::cout<<"result: " << result << "\n";
std::cout<<"abserr: " << abserr << "\n";
std::cout<<"nevals: " << nevals << "\n";
gsl_integration_cquad_workspace_free(workspace);
return 0;
}
