A modernized SPLPAK

This post brought NCL to my attention, which I’d not heard of before. In the ngmath folder you can find a lot of old Fortran code for interpolation. See here for an ngmath overview. I guess at one point it was a standalone library. It looks like development of this code stopped decades ago (some of it looks like Fortran 66). But, there are some nice routines in there. I noticed some that indicated they were from something called SPLPAK (a Pack I had never heard of before). It has code for computing N-D cubic spline fitting using a least squares method. The comments indicate it was started in 1972.

So, I present a modernized SPLPAK. Improvements include:

  • Converted the Fortran 66 code to modern Fortran. Eliminated all GOTOs and other ancient nonsense. The code is now much more readable.
  • Eliminated the hard-coded limit of only 4 dimensions. Now it will work with any number of dimensions. This may now be the only publicly-available least squares cubic spline Fortran code that has this (does anybody know of another?)
  • You can now specify the real kind using a compiler directive. So, if you want real128, you got it.
  • Made it threadsafe. The original code was not, since it included save statements and common blocks. There are no more global variables, and all data is stored in a splpak_type, which includes initialize and evaluate methods for computing and evaluating the splines.
  • It can be used with the Fortran Package Manager.
  • Documentation is now auto generated from the code using FORD.

It’s great to see this code polished up and ready for another 50 years of service, rather than hidden away in a deep subdirectory of a mostly-abandoned software package.

15 Likes

fitpack does have a hardcoded limit of 10 dimensions. Extending it too is actually a great idea!

1 Like

I was considering in approaching this problem myself. I think it was necessary.

Many thanks!

Nice work Jacob. Looking through the NCL code there appears to be some of the old NCAR graphics routines (in Fortran) that might be useful. I’m retireing in a few months so maybe I’ll give it a look (can only play so much golf and catch so many fish).

1 Like

Fantastic.

This supersedes it beautifully, but for reference the abstract (1991 version) of the previous SPLPAK read …


SPLPAK - Carl deBoor's cubic spline package.

Splines are  functions used to represent user data.  Typically
the spline will be built up out of polynomial "pieces", with some of
the higher derivatives not continuous at the "joins".  The most common
spline used is the piecewise cubic.  Splines are efficient and accurate.
In particular, they are far better than using high order polynomial
approximation.  

Applications included in the package will handle two-dimensional data, 
two-point boundary value problems, least squares approximation, integration, 
differentiation, smoothing, and other topics.
 
Reference:    Carl de Boor, 
              A Practical Guide to Splines, 
              Springer Verlag
              ISBN: 0-387-90356-9

I would have thought it deserved a new name like SPLINEPAK to avoid confusion with the original but looking around that does not seem to be an issue, so I welcome back SPLPAK, new and improved.

I think that SPLPAK is different from this one. This one’s says it’s from Dave Fulker at NCAR.

NCAR seems to have once had a lot of Fortran libraries. There are traces of them on the internet but also a lot of dead links. So it’s hard to know the provenance of this.

The NCL git history is fascinating. It has all the commits all the way back to 1990 (looks like it went through cvs/svn/git conversion over the years). It predates the first release of Python. In all that time, there was Fortran code in there, totally unchanged. SPLPAK is even older, from 1972. When FORTRAN 77 came along and added if/else/endif constructs, it was not updated. When Fortran 90/95 came along and added free-form source, it was not updated. When Fortran 2003 came along and added classes, it was not updated. The same story we have seen for many classic Fortran codes.

If it ain’t broke, don’t fix it?

And now they are throwing it all in the trash and replacing it with Python.

It took me a few days to do all the upgrades that could have been done over the course of 40 years if this code had been properly maintained. Not really a lot of work actually. The new version still works great, and now has features available that were not possible in the old version. The code is also no longer a scene of ancient horrors that would scare away a potential Fortran users. :slight_smile:

1 Like

NCL was used (and still is by some) a lot in my field, including my old PhD lab. When I came in as a student I replaced a MATLAB-based graphics pipeline that was running operationally during hurricane seasons. The pipeline was both slow and producing ugly plots–I replaced it with a set of Fortran programs and GrADS scripts that were glued with shell. This ran very fast and was a bit less ugly. However GrADS is quite minimal feature-wise and takes a lot of work/code to make plots look good. Next year a student in the lab volunteered (to my dismay) to implement the graphics pipeline in NCL. It was slower than GrADS, faster than MATLAB, and the plots looked OK (still not great). A few years later I re-implemented the operational system in Python and replaced some of the NCL graphics with Matplotlib. Some, but not all. For NCL had one feature that no other graphics software that I know of had–curly vectors. Curly vectors are streamlines with arrowheads whose length and density was scaled by the magnitude of the vector field. I implemented these in Matplotlib once but couldn’t nearly meet the speed of NCL.

NCL was terrible and wonderful. It’s the kind of software that people either loved or hated. I didn’t love it but know many who did. It sure had its time and place.

2 Likes

For fast Python graph, typically in user interfaces, I use pyqtgraph that’s much faster than matplotlib.

1 Like