# Modern Fortran Quadpack **URL:** \ **Category:** Announcements\ **Created:** [June 10, 2022, 3:22pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714 "2022-06-10T15:22:32Z")\ **Posts on this page:** 16\ **Page:** 1
**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\ **Post date:** [June 10, 2022, 3:22pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/1 "2022-06-10T15:22:32Z")
I think it’s been mentioned in other posts, but I wanted to officially announce the new, Modern Fortran Quadpack library for numerical integration (quadrature): > **[GitHub - jacobwilliams/quadpack: Modern Fortran QUADPACK Library for 1D numerical...](https://github.com/jacobwilliams/quadpack)** > > Modern Fortran QUADPACK Library for 1D numerical quadrature This is a complete modernization and refactoring of the old library from [Netlib](http://www.netlib.org/quadpack/). Changes include: - It has been converted from FORTRAN 77 fixed form to modern free form syntax. This includes elimination of all GOTOs and other obsolescent language features. - It is now a single stand-alone module, and has no dependencies on any other code from SLATEC or LINPACK. - It is a Fortran Package Manager package. - The separate routines in the original library for single and double precision have been eliminated. The library now exports a single (`real32`), double (`real64`) and quadruple (`real128`) precision interface using the same code by employing a preprocessor scheme. - The coefficients have been regenerated with full quadruple precision. This was done using [a new program](https://github.com/jacobwilliams/kronrod) that employed the [MPFUN2020](https://github.com/jacobwilliams/mpfun2020-var1) arbitrary precision Fortran library. Quad-precision versions of these routines are available in no other library that I’ve been able to find. - Some minor bugs have been fixed in the original code that were in there for decades (see [here](https://github.com/scipy/scipy/issues/14807), [here](https://github.com/scipy/scipy/pull/14836), [here](https://github.com/jacobwilliams/quadpack/issues/1), and [here](https://github.com/jacobwilliams/quadpack/issues/7) ). - Other procedures not present in the original QUADPACK have been added (new routines, and modernized ones from old libraries such as [SLATEC](http://www.netlib.org/slatec/) and the [NSWC Library](https://github.com/jacobwilliams/nswc)): `QUAD`, `AVINT`, `QNC79`, `GAUSS8`, `SIMPSON`, and `LOBATTO`. - The SLATEC docstrings have been converted to [Ford](https://github.com/Fortran-FOSS-Programmers/ford) style, which allows for auto-generation of the [API docs](https://jacobwilliams.github.io/quadpack/). - Some typos, there for decades, have been corrected in the comments. - It’s unit-tested with GitHub Actions CI The goal here is to restart Quadpack development where it left off 40 years ago. Bugs can be fixed, and new routines can be added. This can be the standard state-of-the-art library for numerical quadrature once again, second to no other library in any other programming language. There’s no reason to be stuck with the old Netlib FORTRAN 77 code anymore. Check the GitHub issues to see some ideas on new methods that we want to add. Join us!! Just like MINPACK, maybe we eventually move this under the fortran-lang umbrella, and also maybe try to get it into SciPy? @certik what do you think? ### See also - [Packs - The Classic Fortran Libraries](https://fortran-lang.discourse.group/t/packs-the-classic-fortran-libraries/2751) ---
**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\ **Post date:** [June 11, 2022, 5:06am UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/2 "2022-06-11T05:06:43Z")
Thanks! I think this is great. Yes, we should move under fortran-lang. I think minpack is more ahead in terms of modernization, so we should try that first into SciPy. Then the other libraries. ---
**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\ **Post date:** [June 12, 2022, 2:28pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/3 "2022-06-12T14:28:14Z")
Just curious re: the following: 1. Is `Modern Quadpack` part of Fortran `stdlib`? If not, can it not be? 2. What are the implications of inquiry 1 above? 3. Say `Modern Quadpack` is further refactored to use the `kinds`, `math` constants, etc. from Fortran `stdlib`: prima facie, it appears a good step forward to me but are there cons about which I am unaware other than dependencies on `stdlib` which is growing to be quite big? ---
**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\ **Post date:** [June 13, 2022, 2:42am UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/4 "2022-06-13T02:42:26Z")
It is not currently part of `stdlib`. I’m not sure it needs to be. With FPM, we can have a bunch of standalone libraries and you can just pull down the ones you need, no problem. I did not use the `fypp` thing that `stdlib` is using so I could stick with standard fortran, so my IDE/linter/etc will work, which I have come to depend on. But, if you look at the code, you can see that it publishes `real32`, `real64`, and `real128` versions of all the routines. So, from a user point of view, you get the same thing. And from a developer’s point of view, you have no code duplication (just a little boilerplate) and a file that editors/linters can understand as Fortran. Also, it’s not clear if SciPy would accept an `stdlib` dependency. I don’t know. Putting it under fortran-lang I think just gives it (maybe) an air of legitimacy. That’s really the only reason for that. I agree about waiting to see if we can get Minpack in, and then go from there. ---
**Author:** ![sshestov](https://avatars.discourse-cdn.com/v4/letter/s/bcef8e/32.png) [@sshestov](https://fortran-lang.discourse.group/u/sshestov)\ **Post date:** [June 29, 2026, 3:57pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/5 "2026-06-29T15:57:02Z")
A naive question: I’m trying to use `qawo` while working with double precision (`use iso_fortran_env, only: wp=>real64` ). My g(x) function is defined as ```fortran real(wp) function qawo_f(x) real(wp), intent(in) :: x qawo_f = BESSEL_J0(2 * dpi * rAplane*sqrt(x)*lamz0_1) end function qawo_f ``` dpi, rAplane and lamz0\_1 are defined in my host program. The question is: should I use `dqawo(qawo_f,a,b...) ` and not `qawo(qawo_f,a,b...)` Is it always like this when using procedures as arguments? ---
**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\ **Post date:** [June 29, 2026, 6:59pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/6 "2026-06-29T18:59:27Z")
As long as the compiler can disambiguate the interface, then the generic interface name should suffice. If it cannot, you will get a compiler error. Since the argument is a function, the compiler should be able to determine what specific routine “qawo” to use from the given arguments and the return value of the function. Are you getting an error message? ---
**Author:** ![sshestov](https://avatars.discourse-cdn.com/v4/letter/s/bcef8e/32.png) [@sshestov](https://fortran-lang.discourse.group/u/sshestov)\ **Post date:** [June 30, 2026, 7:19am UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/7 "2026-06-30T07:19:23Z")
Yes, the error is here (I’m running it with fpm) ```bash 229 | call qawo(qawo_f,0._wp,(710.0_wp)**2,dpi*lamz0_1,1 ,1d-2 ,1d-4 ,val , eps, neval, ierr, leniw, maxp1, lenw, last, iwork, work) | 1 Error: Interface mismatch in dummy procedure ‘f’ at (1): Type mismatch in function result (REAL(4)/REAL(8)) ``` in the invocations everything is double. `qawo_f` is defined in my post above. Switching to `dqawo` immediately solves the problem of compilation. If can use `qawo` if I create another routine `real function qawo_f4(x)` and change the arguments to real in the call. That’s why I have this question ---
**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\ **Post date:** [June 30, 2026, 7:40am UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/8 "2026-06-30T07:40:05Z")
Hm, is your function qawo\_f in a module? I have written a small sample program that shows this should work (see the attached file) [use\_func.f90](https://fortran-lang.discourse.group/uploads/short-url/1CE2dPcHDel5CEptNBe8g9NAfur.f90) (907 Bytes) , but it relies on the function’s actual interface to be visible to the compiler. Can you show us a minimal example? If your function has an _implicit_ interface, then the return type will be default (single-precision) real and that could cause the mismatch. ---
**Author:** ![sshestov](https://avatars.discourse-cdn.com/v4/letter/s/bcef8e/32.png) [@sshestov](https://fortran-lang.discourse.group/u/sshestov)\ **Post date:** [June 30, 2026, 8:46am UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/9 "2026-06-30T08:46:26Z")
Thanks for your comment and example. In my case qawo\_f is contained within the main program. I think it means explicit interface. I will check yours and provide a MWE. ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [June 30, 2026, 4:10pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/10 "2026-06-30T16:10:36Z")
> [@sshestov](#): > > Switching to `dqawo` immediately solves the problem of compilation. It looks like `qawo()` is expecting `real32` arguments and `dawo()` is expecting `real64` arguments. If there is a generic interface `qawo()` for both, then you can use that and get the right specific routine just based on the argument types. Otherwise, you will need to use manually the correct specific routine to match your arguments. ---
**Author:** ![sshestov](https://avatars.discourse-cdn.com/v4/letter/s/bcef8e/32.png) [@sshestov](https://fortran-lang.discourse.group/u/sshestov)\ **Post date:** [July 1, 2026, 9:52am UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/11 "2026-07-01T09:52:41Z")
That is exactly (generic interface) what I was expecting and trying to understand ---
**Author:** ![sshestov](https://avatars.discourse-cdn.com/v4/letter/s/bcef8e/32.png) [@sshestov](https://fortran-lang.discourse.group/u/sshestov)\ **Post date:** [July 1, 2026, 9:54am UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/12 "2026-07-01T09:54:39Z")
Here I provide a MWE. Switching from real32 to real64 prevents compilation; if I use real64, I have to switch to dqawo. ```fortran program main_test use iso_fortran_env, only: wp=>real32, real64, real32 use, intrinsic :: ieee_arithmetic use quadpack, only: qawo, dqawo implicit none real(kind=wp),PARAMETER :: dpi=3.141592653589793238462643383279502884197169399375105d0 real(kind=wp) :: lamz0_1 real(kind=wp) :: eps, val, r_Aplane integer :: ierr !! **********************variables for QAWO method********************** integer, parameter :: limit=2000 integer, parameter :: leniw = limit*2 integer, parameter :: maxp1 = 41 integer, parameter :: lenw = limit*4 + maxp1*25 real(kind=wp) :: work(lenw) integer :: iwork(leniw), last, neval, nevalM lamz0_1=1.2596d-2 r_Aplane=1.5d0 write(*,*) " ****************test of QuadPack's QAWO method:*******************" call qawo(qawo_f,0._wp,100.0_wp,real(dpi*lamz0_1,kind=wp), 1, 1e-2_wp, 1e-4_wp, val, eps, neval, ierr, leniw, maxp1, lenw, last, iwork, work) print*, "Ierr (0==success):", ierr print*, "neval:", neval print*, "Working precision wp: ", wp write(*,*) " " contains real(wp) function qawo_f(x) real(wp), intent(in) :: x qawo_f = BESSEL_J0(2.*dpi*r_Aplane*sqrt(x)*lamz0_1) end function qawo_f end program main_test ``` ---
**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\ **Post date:** [July 1, 2026, 10:58am UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/13 "2026-07-01T10:58:09Z")
I think I understand it. The package does not define generic interfaces, instead it uses separate names. You could solve this by simply defining your own generic interface. Something along these lines: ```auto module my_quadpack use quadpack_single, only: sqawo => dqawo use quadpack_double, only: dqawo interface qawo module prodedure :: sqawo, dqawo end interface end module myquadpack ``` and use that instead of the original quadpack. ---
**Author:** ![sshestov](https://avatars.discourse-cdn.com/v4/letter/s/bcef8e/32.png) [@sshestov](https://fortran-lang.discourse.group/u/sshestov)\ **Post date:** [July 1, 2026, 12:19pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/14 "2026-07-01T12:19:01Z")
Many thanks, now it’s clear! (I’m just starting with fpm; was suspecting I was doing some stupid things) ---
**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\ **Post date:** [July 1, 2026, 12:26pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/15 "2026-07-01T12:26:40Z")
The way this package has been set up is a bit unusual in my view, I had to examine the source code to unravel the problem 😇. A welcome procrastination … ---
**Author:** ![sshestov](https://avatars.discourse-cdn.com/v4/letter/s/bcef8e/32.png) [@sshestov](https://fortran-lang.discourse.group/u/sshestov)\ **Post date:** [September 17, 2026, 2:22pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/16 "2026-09-17T14:22:41Z")
Dear all, I’m having troubles with OpenMP-parallelized do-loops with `QAWO` or `QAGS`. My minimal (not working) example is below ```fortran !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(r, r_Aplane, int_cos, int_sin, f_cos, f_sin, val, eps, neval, ierr, last, iwork, work) do j=1, 500 r=A_axis(j) !! A-plane coord. !! ************ my straightforward sampling of the 1st plane works well with OpenMP f_cos = 2. * dpi * r_axis * cos(dpi*r_axis**2 *lamz0_1) * BESSEL_J0(2.*dpi*r*r_axis*lamz0_1) f_sin = 2. * dpi * r_axis * sin(dpi*r_axis**2 *lamz0_1) * BESSEL_J0(2.*dpi*r*r_axis*lamz0_1) int_cos = integrate_trap(real(r_axis,kind=8), f_cos) int_sin = integrate_trap(real(r_axis,kind=8), f_sin) !! ********but I want to implement it via QUAD_PACK routines************************** !! **************variable substitution in the integral ksi <=> r_axis***************** !! Int ( ksi*cos( .. ksi^2) * J0(2pi r ksi/lamz0) dksi ) => !! 1/2 Int ( cos(ksi^2) * J0(2pi r sqrt(ksi^2) / lamz0) dksi^2 r_Aplane=r call dqags(qags_fc,0._wp,710.0_wp,1d-2, 1d-4 ,val, eps, neval, ierr, limit, lenw, last, iwork, work) int_cos = val*0.5*2. ! 2 to conform the definition in my original code r_Aplane=r call dqags(qags_fs,0._wp,710.0_wp,1d-2, 1d-4 ,val, eps, neval, ierr, limit, lenw, last, iwork, work) int_sin = val*0.5*2. !!!! temporarily put newly computed real part to imaginary or vise versa !!psi_imag(j) = 1. - ( cos(dpi*r**2*lamz0_1)*int_sin + sin(dpi*r**2*lamz0_1)*int_cos ) * lamz0_1 psi_real(j) = ( cos(dpi*r**2*lamz0_1)*int_cos - sin(dpi*r**2*lamz0_1)*int_sin ) * lamz0_1 if (mod(j,100) == 0) write(*,'(I4,A)', advance='no') int(j/100)," " end do !$OMP END PARALLEL DO ... contains real(wp) function qags_fc(x) real(wp), intent(in) :: x qags_fc = 2. * dpi * x * cos(dpi * x**2 * lamz0_1) * BESSEL_J0(2.*dpi*r_Aplane*x*lamz0_1) end function qags_fc end ``` When I use self-implemented trapezoid implementation and explicit sampling (in r\_axis) it works flawlessly. But I wanted to switch to QUAD\_PACK routines or a separate Levin method, with a hope that presence of cos/sin terms can speed up. I’ve started from Levin method - it was working but giving spurious outliers; then QAWO was working only until some j. I guess this is because Bessel J0 is also oscillating. Finally I switched to QAGS, which seems to work, but I struggling with OpenMP. It works, but gives completely wrong results. Do you have any ideas/suggestions? PS: this is to calculate diffraction behind a circular occulter.