Add QUADPACK to `fpm` registry

[dhermes]

I actually would be willing to do all of the work here but am not sure where to start.

AFAICT the original source has been converted to Fortran 90 by John Burkardt (single and double precision). The big concerns then are (A) LICENSE, this is the biggest (B) build tooling e.g. Make or CMake, not sure what fpm prefers (C) documentation.

[ivan-pi]

QUADPACK was included in SLATEC and is thereby in public domain (at least in the USA). See this comment from the Scipy repository. The website from Leuven University also mentions QUADPACK is included in Octave and that C translations are available in GSL. So open source usage seems to be permitted.

The version of QUADPACK included in SciPy has some fixes, like replacing calls to linpack routines with their respective lapack equivalents. Personally, I would prefer to use the Scipy version over those of Burkardt. The Intel Fortran compiler could be used to automatically generate the missing interfaces, and leave the source as F77 code. The automatically generated interfaces could then be manually compared to the documentation and published listings to get the intent attributes correct.

Concerning build tools, the idea of the fpm registry is to collect fpm compatible packages. The best approach is to allow fpm to do the building, no make or cmake is required. One hurdle will be how to provide both single and double precision versions. Edit: probably the single precision versions should be downloaded from netlib and prepended with s. Then the single and double precision versions can be placed under a generic interface. It will be necessary to cross-check against the SciPy modifications to bring the single precision versions up to date.

Currently, fpm offers no tools to automate documentation. Most fortran-lang projects are using FORD and GitHub actions to automate the site generation. Personally, I find FORD to be quirky. For the QUADPACK documentation it might be possible to extract the comments from the source files (either the SciPy or SLATEC versions) using a Python/Perl script and generate restructured text for some other documentation website generator (Sphinx, MkDocs, ...).

[certik]

We could consider using the quadpack from SciPy. That way we are fine with the license, we'll have the fixes and it can also be useful from the point of view of SciPy and other such libraries to eventually use the quadpack fpm package.

Also, I would recommend to maintain this under the fortran-lang github organization. Many of the issues facing here are identical to fftpack: #40, see the discussion there.

Let's figure out a good approach to these things.

[dhermes]

RE: SciPy source, the files in https://github.com/scipy/scipy/tree/v1.6.3/scipy/integrate/quadpack are all Fortran 77. The merits of the Burkardt code is that it is in Fortan 90, making it accessible to a wider set of compilers (e.g. lfortran). I am very open here on where the source comes from and I would say NOT very educated on historical forks since the original QUADPACK in 1981, just trying to see if there is something stable AND quasi-modern.

[certik]

Good point. The only reason I suggested SciPy is that I know they maintain it and fix bugs and the code should be correct and working. If we use some older unmaintained version, it might have bugs. Aren't there tools to convert fixed form to free form? We could start with the SciPy version, convert to free form, add to modules, etc.

P.S. LFortran will eventually have a parser for fixed form also (https://gitlab.com/lfortran/lfortran/-/issues/82).

[awvwgk]

I added a couple of formatters to the package index a while ago.

I made good experience with findent, which is primary targeted to parse legacy Fortran and handle dozens of non-standard extensions (packaged on conda-forge here). Camfort and fortran-src also have legacy Fortran parsers and there is fpt from the last monthly call.

So plenty of tools for the job out there, but most important is to have tests to verify the converted code against the original one.

[ivan-pi]

Good point. The only reason I suggested SciPy is that I know they maintain it and fix bugs and the code should be correct and working. If we use some older unmaintained version, it might have bugs.

Not sure if we can call it maintenance. The code has been untouched for the last 3 years and there are some open bug tickets however, only one seems to be an actual numerical issue, the rest look like compiler and build setup issues.

The SciPy version of QUADPACK would also benefit from some refactoring (see scipy/scipy#4455, great issue title!). The QUADPACK routines make use of the following three machine constants:

1. machine epsilon (d1mach(4))
2. the largest non-overflowing positive real number (d1mach(2))
3. the smallest non-underflowing positive real number (d1mach(1))

which are currently retrieved using the d1mach() function (a description can be found in the SLATEC documentation or BLAS faqs). The version used in SciPy looks like it is taken from Netlib, but from some third package (not SLATEC or QUADPACK). In the meantime (30 years ago...) Fortran 90 added intrinsic functions to retrieve these constants directly using intrinsic functions. I am aware of a few different refactored versions of d1mach:

• https://github.com/certik/fortran-utils/blob/master/src/legacy/amos/d1mach.f90
• http://computer-programming-forum.com/49-fortran/9d39e9771b0d8e20.htm
• https://people.sc.fsu.edu/~jburkardt/f_src/machine/machine.html
• one more approach is to make d1mach an array of size 5 (@jacobwilliams perhaps?), and thereby catching wrong inputs at compile time

Interestingly, the SciPy d1mach source has a commented out C implementation:

C implementation

 /* Standard C source for D1MACH -- remove the * in column 1 */
#include <stdio.h>
#include <float.h>
#include <math.h>
double d1mach_(long *i)
{
	switch(*i){
	  case 1: return DBL_MIN;
	  case 2: return DBL_MAX;
	  case 3: return DBL_EPSILON/FLT_RADIX;
	  case 4: return DBL_EPSILON;
	  case 5: return log10((double)FLT_RADIX);
	  }
	fprintf(stderr, "invalid argument: d1mach(%ld)\n", *i);
	exit(1); return 0; /* some compilers demand return values */
}

[LKedward]

Note that Carlos @brocolis already has an fpm package for the original quadpack here which also includes an f90 interface module.

[certik]

Note that Carlos @brocolis already has an fpm package for the original quadpack here which also includes an f90 interface module.

That's exactly what I had in mind! @brocolis would you be interested in moving the package under fortran-lang so that we can maintain it as a community?

Yes. Someone with appropriate permissions could move it to fortran-lang.

My original plan was to write an higher level function quad on top of QUADPACK, but since I read the GNU Octave implementation for quad, I believe I cannot re-implement it myself under a different license. However, I would be happy to transfer the repository to fortran-lang as it is now.

[ivan-pi]

The scipy.integrate.quad has a very similar interface to the Octave quad. The SciPy version is under a BSD-3 license.

But I would suggest we do it jointly under the fortran-lang repository.

[ivan-pi]

However, I would be happy to transfer the repository to fortran-lang as it is now.

I think if you follow the instructions on transferring repositories you could transfer to @milancurcic, and he could put it under the fortran-lang organization.

The other option would be to simply clone the repository, and push it into a new fortran-lang repository (not forking it).