Update to Bicubic Interpolation Code (TOMS 760)

Dr. Antony Lewis spotted a typo in the Fortran 90 translation of TOMS algorithm 760. The file in question (toms760.f90) is part of the archive of Alan Miller’s Fortran Software that I host.

Although I do not actively maintain Alan’s collection of algorithms, in this case I verified that the definition of the matrix idlt in the previous Fortran 90 version did not correspond to the original Fortran 77 implementation. As Dr. Lewis noted in his email to me, this is a simple Fortran column-major ordering mistake.

That is to say, the following definitions of idlt are not equivalent:

integer :: idlt(3,4)
data :: ((idlt(i,j),j=1,4),i=1,3) / -3, -2, -1, 1, -2, -1, 1, 2, -1, 1, 2, 3 /

and

integer :: idlt(3,4) = reshape( (/ -3, -2, -1, 1, -2,-1, 1, 2, -1, 1, 2, 3 /), (/ 3, 4 /) )

The updated version of the file is here. If you’d like to patch your local copy, this is the diff output:

diff --git a/mirror/amiller/toms760.f90 b/mirror/amiller/toms760.f90
index 970f453..24e9d80 100644
--- a/mirror/amiller/toms760.f90
+++ b/mirror/amiller/toms760.f90
@@ -17,6 +17,9 @@ SUBROUTINE rgbi3p(md, nxd, nyd, xd, yd, zd, nip, xi, yi, zi, ier)
 ! Code converted using TO_F90 by Alan Miller
 ! Date: 2003-06-11  Time: 10:11:03

+! Typo in definition of IDLT in SUBROUTINE RGPD3P fixed by Jason Blevins
+! Date: 2015-06-21  Time: 00:22:26
+
 ! Rectangular-grid bivariate interpolation
 ! (a master subroutine of the RGBI3P/RGSF3P subroutine package)

@@ -472,7 +475,7 @@ REAL    :: b00xa(4), b00ya(4), b01a(4), b10a(4), cxa(3,4), cya(3,4),   &
 ! Data statements
 ! DATA ((idlt(jxy,jpexy),jpexy=1,4),jxy=1,3)/-3,-2,-1,1,-2,-1,1,2,-1,1,2,3/
 INTEGER, SAVE  :: idlt(3,4) = RESHAPE(   &
-                  (/ -3,-2,-1, 1,-2,-1, 1,2,-1, 1,2,3 /), (/ 3, 4 /) )
+                  (/ -3,-2,-1, -2,-1,1, -1,1,2, 1,2,3 /), (/ 3, 4 /) )
 !     ..

 ! Calculation