Add QUADPACK to `fpm` registry

[dhermes]

trafficstars

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 https://github.com/scipy/scipy/issues/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?

[ghost]

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

[ghost]

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

[awvwgk]

Repository transfers (and creation) can currently only be handled by @certik, @milancurcic and @LKedward, either by temporarily raising member privileges or by a two step transfer. I don't think we have yet established a workflow for this process.

[certik]

@LKedward, @milancurcic I would be interested in your feedback on maintaining this under fortran-lang.

[LKedward]

@LKedward, @milancurcic I would be interested in your feedback on maintaining this under fortran-lang.

I have no fundamental objections to moving this and other legacy packages to fortran-lang as long as we have one or more persons identified to spearhead the initial modernisation and maintenance required to bring the packages up to standard. (Which it looks like we do here!) In particular, we should strive to maintain a high quality of packages maintained under fortran-lang which I will clarify by quoting from our extra 'recommend criteria' in the package index guidance:

• Documentation: any form of written documentation aimed at users of the package. Ideally this should cover:
  • Supported / tested compilers, dependencies, build and install process, modules contained within the package, procedures made available and their interfaces, example code
• Contributing: details on how users may submit issues and contribute to the development of the package
• Tests: any form of executable test(s) that can be used to verify the functionality of the package
• Portability: no non-standard language extensions or proprietary dependencies
• fpm: support installation by the Fortran Package Manager fpm

[awvwgk]

Do we know of any users of QUADPACK, who are interested in devendoring their version in favor of a @fortran-lang one? Is anyone involved in this thread using QUADPACK in a production code and willing to offer their code base for testing?

I'm also wondering about the QUADPACK license, is it public domain since it is part of SLATEC? How would we handle and license contributions to this project?

[arjenmarkus]

I am not a user of QUADPACK, but I have experimented with it. In the various versions I have (one from Netlib and one by John Burkhardt) I see no regular test programs, like you find with LAPACK, but it might be a nice community project to set up a comprehensive set of such programs. Within the framework of fpm we can make them available. And perhaps en passant provide modern interfaces.

Op do 27 mei 2021 om 15:41 schreef Sebastian Ehlert < @.***>:

Do we know of any users of QUADPACK, who are interested in devendoring their version in favor of a @fortran-lang https://github.com/fortran-lang one? Is anyone involved in this thread using QUADPACK in a production code and willing to offer their code base for testing?

I'm also wondering about the QUADPACK license, is it public domain since it is part of SLATEC? How would we handle and license contributions to this project?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/fpm-registry/issues/43#issuecomment-849644531, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR5BQVL3LCLOZZSJDKLTPZDX3ANCNFSM446N37RQ .

[ivan-pi]

Is anyone involved in this thread using QUADPACK in a production code and willing to offer their code base for testing?

Well there are the possibly thousands of users using it via SciPy and Octave. Also the NAG quadrature routines are based upon QUADPACK. Also the IMSL library seems to be based off QUADPACK. Several papers discussing QUADPACK (with some tests) can also be found in

Keast, P., & Fairweather, G. (Eds.). (1987). Numerical integration: recent developments, software and applications. https://doi.org/10.1007/978-94-009-3889-2

SciPy has a set of tests for QUADPACK that can be adapted: test_quadpack.py.

The QUADPACK book (https://link.springer.com/book/10.1007%2F978-3-642-61786-7) contains error plots of 17 test integrals. There are also 10 demonstration programs. I would suggest to add these as the original tests.

[certik]

SciPy is the most obvious one. Yes, some of the SciPy maintainers are open to the idea of using fortran-lang's version, but of course it would depend on the details. I have not seen any objections like "absolutely not". So that means we should do our best, and if we do, there is a good chance of SciPy using it, one way or another.

[arjenmarkus]

I must say the book by Keast en Fairweather is a bit expensive (200 euros for the e-book version). The other one a bit less so.

By looking at the WIkipedia page I found an update of two of the automatic integrators in QUADPACK, namely QAG and QAGS, as well as an integrator for multidimensional functions (TOMS algorithms 691 and 824).

Should there be interest, I also have the implementaton of quadrature (cubature) formulae by Stroud - a vast collection of routines.

Op do 27 mei 2021 om 17:32 schreef Ivan Pribec @.***>:

Is anyone involved in this thread using QUADPACK in a production code and willing to offer their code base for testing?

Well there are the possibly thousands of users using it via SciPy and Octave. Also the NAG quadrature https://www.nag.com/numeric/cl/nagdoc_latest/html/d01/d01conts.html routines are based upon QUADPACK. Also the IMSL library http://envi.geoscene.cn/help/Subsystems/idl/Content/Reference%20Material/Functional%20List%20of%20IDL%20Routines/IMSLRoutinesQuadrature.htm seems to be based off QUADPACK. Several papers discussing QUADPACK (with some tests) can also be found in

Keast, P., & Fairweather, G. (Eds.). (1987). Numerical integration: recent developments, software and applications. https://doi.org/10.1007/978-94-009-3889-2

SciPy has a set of tests for QUADPACK that can be adapted: test_quadpack.py https://github.com/scipy/scipy/blob/4ec4ab8d6ccc1cdb34b84fdcb66fde2cc0210dbf/scipy/integrate/tests/test_quadpack.py .

The QUADPACK book ( https://link.springer.com/book/10.1007%2F978-3-642-61786-7) contains error plots of 17 test integrals. There are also 10 demonstration programs. I would suggest to add these as the original tests.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/fpm-registry/issues/43#issuecomment-849731558, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR7WH6HGC54LW3UI5Q3TPZQYHANCNFSM446N37RQ .

[certik]

So the next step, as @LKedward said, is to find a person who would be willing to take the lead on this. @dhermes and @brocolis, would you have time for that? The aim should be to maintaining this as a community, under fortran-lang. But it's good if at least initially there is one or two people who push the community to get it done, i.e., takes the lead on this.

One possible path forward is:

• Take the repository https://github.com/brocolis/quadpack by @brocolis.
• Figure out the license situation. If in doubt, we might want to just take the version from SciPy, which is under BSD license.
• Move it under fortran-lang.
• Improve the tests, per @ivan-pi's pointers
• Improve documentation
• Apply any fixes

What should be the long term goal? Should we include integration routines that we develop, such as by @arjenmarkus as indicated above? I would suggest to define what the long term objective is. I feel we should probably just "maintain it", i.e. ensure it compiles without warnings with all Fortran compilers, fix any bugs if they are discovered, and maintain modern Fortran and C API to it. But other than that, I think new algorithms should go into a new/separate fpm package. What do you think?

[jacobwilliams]

I would highly recommend not moving F77 code to fortran-lang. Or wrapping F77 with a modern interface. I think that is a bad idea. We need to modernize the code to show what can be done with modern Fortran. There is a lot of room for improvement in quadpack. We don't want to be stuck with the old F77 code like some sort of new netlib. Convert it to free-form and then go from there.

EDIT: For example, here's some modernized SLATEC integration code: https://github.com/jacobwilliams/quadrature-fortran

[certik]

I think we should discuss this at the next Fortran call.

@jacobwilliams I think we all agree that we want modern Fortran implementations of these basic algorithms, whether as standalone packages and/or part of stdlib. Many of such packages could be maintained under fortran-lang also.

The main problem that we are trying to solve in this thread is that codes depend on these old packages (fftpack, odepack, minpack, quadpack, ...) so having a maintained version helps. Also to compare against the "modern" packages, to ensure that they are actually better, or at least as good, i.e., did not regress on speed, accuracy or "simple API" (subjective).

The question is whether the modernized version should be part of the same package, replacing the old code (I think that is what you are suggesting?), or a separate package.

I think as long as performance and the API can be kept exactly the same, it is fine to be part of the same package. If the API changes, then I think it should probably be a separate package.

As another example, I took fftpack and modernized it here:

https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/fourier.f90#L201

Things to keep in mind when modernizing:

• performance: did it get slower? To answer it, I think we want to have the original package, easy to install, to compare against
• accuracy: is it as accurate?
• API: is the new API easier to use?

I have benchmarked my new "modern Fortran" fft implementation, and on the one benchmark I used (about million element 1D fft), it was about the same. But later I realized that for smaller fft, such as 128 (which is used for 3D 128^3, very common in electronic structure calculations), it might be slower. Based on my experience, I think it is actually really hard to ensure the new code is as fast.

For all these reasons, I think it would be valuable to "conserve" or "preserve" these old libraries, and only do changes that keep performance, accuracy and API. Moving from fixed to free form is fine I think. Keeping the performance is tough, for that we need good benchmarks, or not change the code significantly. Keeping the API is simple (we just need tests), but that does not allow to modernize the API.

[ivan-pi]

I must say the book by Keast en Fairweather is a bit expensive (200 euros for the e-book version). The other one a bit less so.

For the time-being, I have access to digital versions of both books through a university subscription. I am happy to extract the tests from the books.

One problem with modifying the API is that old codes that already use QUADPACK will require refactoring to use our "maintained" version. I think this is not constructive. Any older codes should continue to work, with the minor exception that we ask programmers to add a use quadpack, only: ... statement.

My personal goals for the modernization would include:

1. Provide an interface module (already done by @brocolis)
2. Provide the missing tests, build tools, and online documentation
3. Refactor the original sources to free form, removing any obsolescent language constructs, and following any modern computational/software practices.

Performing step 3 after adding tests, helps guarantee accuracy and performance are preserved. The results and timings of the tests, should be preserved as build artefacts.

Only after these steps are complete would I be in favor of adding new routines with simplified API's (through use of optional arguments, or derived types).

[arjenmarkus]

@Ivan Pribec @.***> I agree with your list of goals. And to prevent lots of work on old code, we might provide an old-style interface to the modern version of QUADPACK or any other venerable package.

Op do 27 mei 2021 om 23:59 schreef Ivan Pribec @.***>:

I must say the book by Keast en Fairweather is a bit expensive (200 euros for the e-book version). The other one a bit less so.

For the time-being, I have access to digital versions of both books through a university subscription. I am happy to extract the tests from the books.

One problem with modifying the API is that old codes that already use QUADPACK will require refactoring to use our "maintained" version. I think this is not constructive. Any older codes should continue to work, with the minor exception that we ask programmers to add a use quadpack, only: ... statement.

My personal goals for the modernization would include:

1. Provide an interface module (already done by @brocolis https://github.com/brocolis)
2. Provide the missing tests, build tools, and online documentation
3. Refactor the original sources to free form, removing any obsolescent language constructs, and following any modern computational/software practices.

Performing step 3 after adding tests, helps guarantee accuracy and performance are preserved. The results and timings of the tests, should be preserved as build artefacts.

Only after these steps are complete would I be in favor of adding new routines with simplified API's (through use of optional arguments, or derived types).

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/fortran-lang/fpm-registry/issues/43#issuecomment-849969361, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YRZLGHWNU7DTFS2HSPLTP26CJANCNFSM446N37RQ .

[jacobwilliams]

We need to modernize the API. We can't hold up something that requires you to stuff your problem data into iwork and rwork arrays as any kind of modern way to write code. This does great damage to Fortran anytime anyone sees code like this. It's just a fact. If somebody wants to use the old F77 code and never wants to change anything, that's fine, but those people should not be holding back progress (they have held back progress for decades now... it's time to move on). Make a modernized version, and let people see that it is better and they should use it, and if somebody is still using the old one they have an incentive to use the new one because it is better.

[certik]

@jacobwilliams, we all agree that we want to use modern API to these legacy libraries. That is not a question.

The question is how to get there. You said not to move the "F77 code to fortran-lang." You also suggested to rather work on a modern Fortran replacement, a "modernized version", and it is my understanding that you do not want to use the legacy code at all, but rewrite from scratch?

I think from a practical perspective, I want to have these legacy libraries maintained, because people use them. And I would like to be able to use them at least for performance and accuracy comparison against a "modernized version" that you propose.

To conclude, are you against the plan here: https://github.com/fortran-lang/fpm-registry/issues/43#issuecomment-849969361 ?

That plan still allows to have a separate library, written from scratch, in modern Fortran.

[awvwgk]

If the association of @fortran-lang with the legacy QUADPACK is a blocking issue here, than I might have a suggestion to move this forward. https://github.com/quadpack is currently free, I suggest we claim the @quadpack namespace by creating a new open-source organization and move the legacy version there. This allows to cooperatively develop the QUADPACK under a neutral namespace and once we have reached a point were we want to associate our quadpack-ng with @fortran-lang we can move it here.

I personally use this strategy a lot for my projects, once they mature I tend to move them into separate organization (see @toml-f or @dftd4), since it allows me to separate the stable version from my development branches in a fork and also gives a nice namespace to include additional example packages as well.

[jacobwilliams]

I absolutely support using the legacy code. There's no need to rewrite it from scratch. All I was warning against was being locked into the old api or coding style forever just to make sure that it could be used by somebody who refuses to make even the smallest modifications to their code to use it.

Consider how I made a modernized version of some these spline interpolation routines: https://github.com/jacobwilliams/bspline-fortran

At its core, are routines from CMLIB and SLATEC. I didn't have to rewrite them from scratch, but they have been changed. I had to change them in order to make the interfaces more modern, allow them to be reentrant, add features that weren't present in the original code, etc. At no point did I worry about keeping the old interface. Note that there are even two interfaces, an object oriented one (that I use), but also a subroutine one that is closer (but not identical) to the old interfaces.

So like all these sorts of codes my advice is:

• convert fixed to free form.
• put them in a module or modules
• remove obsolete features: goto, common, data, iwork(), rwork(), etc.

Then you can start expanding beyond what the old codes had (maybe adding new algorithms, maybe an object-oriented interface if it makes sense, etc.)

[certik]

@jacobwilliams perfect! I think the actual plan you just wrote is identical to https://github.com/fortran-lang/fpm-registry/issues/43#issuecomment-849969361, isn't it?

The only technical issue left is "put them in a module or modules", which I would very much prefer, but I don't know if this can cause issues in existing codes like SciPy. But we can help them move over to modules.

@jacobwilliams would you be against moving quadpack under fortran-lang, and we would do this work there?

What @awvwgk suggested above (quadpack organization) also works, but remember that there are tons of other libraries in the same category (fftpack, minpack, odepack, etc.), so we could create fortran-lang-legacy organization and put it there, but we might as well just put it under fortran-lang.