Modern Fortran Quadpack

Announcements
1

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):

This is a complete modernization and refactoring of the old library from Netlib. 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 that employed the MPFUN2020 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, here, here, and here ).

  • Other procedures not present in the original QUADPACK have been added (new routines, and modernized ones from old libraries such as SLATEC and the NSWC Library): QUAD, AVINT, QNC79, GAUSS8, SIMPSON, and LOBATTO.

  • The SLATEC docstrings have been converted to Ford style, which allows for auto-generation of the API docs.

  • 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

29 Likes
2

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.

5 Likes
3

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?
1 Like
4

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.

3 Likes
5

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

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?

6

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?

7

Yes, the error is here (I’m running it with fpm)

  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

8

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 (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.

9

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.

10

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.

11

That is exactly (generic interface) what I was expecting and trying to understand

12

Here I provide a MWE. Switching from real32 to real64 prevents compilation; if I use real64, I have to switch to dqawo.

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
13

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:

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.

14

Many thanks, now it’s clear! (I’m just starting with fpm; was suspecting I was doing some stupid things)

15

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 :innocent:. A welcome procrastination …