Wrong result when calling quadpack from C

Help
1

Attachment

I am a beginner in Fortran. Recently, I would like to use QUADPACK(scipy/scipy/integrate/quadpack at main · scipy/scipy · GitHub) in C.

CMakeList.txt
cmake_minimum_required(VERSION 3.15)
project(quadpack Fortran C CXX)
enable_language(Fortran)
include_directories(include)
file(GLOB QUADPACK_SOURCES “src/*.f”)
add_library(quadpack STATIC ${QUADPACK_SOURCES} src/xerror.f90)

call in Fortran
add_executable(ex_dqng example/ex_dqng.f90)
target_link_libraries(ex_dqng quadpack)

call in C
add_executable(c_dqng test/c_dqng.c)
target_link_libraries(c_dqng quadpack)

Now I define a function f=\frac{2}{2+\sin(10\pi x)} and compute the integral

\int_0^1f(x)dx

The exact solution for the above integral is \frac{2}{\sqrt{3}}

However, the “call in Fortran” gives the correct solution rather the “call in C” gives the wrong result. Please see screenshots I post.

image
image1324×621 45.5 KB

image

image
image845×358 18.7 KB

image
image693×671 22.6 KB

Dev Environment

  • Windows 10 + MinGW 8.0.1
  • Ubuntu 16.4

image
image930×548 111 KB

In addition, I also notice that in Fortran
double precision, parameter :: ans = 2.0 / sqrt(3.0)
write(*, ‘(A, F20.17)’) 'Expected result: ', ans
! Expected result: 1.15470051765441895

image

In contrast, in C
printf(“2.0/sqrt(3.0) = %.20f\n”, 2.0 / sqrt(3.0));

image

image

1 Like
2

Welcome @Lavender (to the world of confusing precision…)

Apart from precision, I am wondering how integral_fun in the C code is declared. Also, the output of ier in both Fortran and C is 1. Is this OK…? (the manual seems to say it should be 0 for success.)

3

Thanks for your reply.

The integral_fun is defined as below:

typedef double (*integral_fun)(double);

In this case (call in Fortran), dqag indeed return ier=1. In my option, all the ouput should be same if I call it in C.
However, they give the different output which is confusing.

4

In Fortran, float literals are single precision by default (the opposite of C and C++). You need to write 2.0d0/sqrt(3.0d0) for double precision literals.

3 Likes
5

For comparison, this is the output from scipy.integrate.quad:

>>> ret = scipy.integrate.quad( lambda x: 2.0 / ( 2.0 + np.sin(10.0 * np.pi * x)), 0.0, 1.0, epsabs=0, epsrel=1.0e-3, full_output=True ); ret[0], ret[1], ret[2]['neval']
(1.1547005383813598, 1.555057807200999e-05, 315)

>>> ret = scipy.integrate.quad( lambda x: 2.0 / ( 2.0 + np.sin(10.0 * np.pi * x)), 0.0, 1.0, full_output=True ); ret[0], ret[1], ret[2]['neval']
(1.1547005383792537, 3.7604956188026965e-09, 483)

>>> 2.0 / np.sqrt(3.0)
1.1547005383792517

One reason for the difference will be that “1.0” etc in Fortran is single precision by default (so we need to write it as “1.0d0” etc), but I wonder why ier = 1 for the output…

6

This callback takes the input by value, but I’m guessing the Fortran callback takes it by reference. (Leaving aside other interoperability issues like name mangling, or the fact your Fortran compiler may use a different argument passing mechanism entirely.)

2 Likes
8

integrate.quad use dqagse, rather than dqng. Please see (scipy/scipy/integrate/__quadpack.h at main · scipy/scipy · GitHub)

1 Like
9

Thanks @ivanpribec, it works when I change 2.0 / sqrt(3.0) to 2.0d0 / sqrt(3.0d0))
image

10

Btw, the GNU Scientific Library also offers the algorithms from QUADPACK with C interfaces.

https://www.gnu.org/software/gsl/doc/html/integration.html

1 Like
11

Thanks for the info about quad()! And I also guess f() may need a pointer, like…

1 Like
12

Thanks a bunch! @septc your guess is right.

image
image594×700 31.6 KB

13

RE ier = 1, I’ve plotted the integrand, which seems pretty oscillatory in [0,1] (because of sin(10 pi x)). And, because qng() is a non-adaptive integrator (according to the following page), it may be simply that qng() does not provide sufficient accuracy for this integrand (unlike a globally adaptive one like qags()).

(I’ve recently used this quad() for my work, so just interested in the internal…)

test