Multidimensional Linear Interpolation (Part 2)

In an earlier post, I mentioned that we needed an object-oriented modern Fortran library for multidimensional linear interpolation. Well, here it is. I call it finterp, and it is available on GitHub. It can be used for 1D-6D interpolation/extrapolation of data on a regular grid (i.e., not scattered). It has a similar object-oriented interface as my bspline-fortran library, with initialize, evaluate, and destroy methods like so:

type(linear_interp_3d) :: interp
call interp%initialize(xvec,yvec,zvec,fmat,iflag)
call interp%evaluate(x,y,z,f)
call interp%destroy()

For example, the low-level 3D interpolation evaluate routine is:

pure subroutine interp_3d(me,x,y,z,fxyz)

implicit none

class(linear_interp_3d),intent(inout) :: me
real(wp),intent(in) :: x
real(wp),intent(in) :: y
real(wp),intent(in) :: z
real(wp),intent(out) :: fxyz  !! Interpolated \( f(x,y,z) \)

integer,dimension(2) :: ix, iy, iz
real(wp) :: p1, p2, p3
real(wp) :: q1, q2, q3
integer :: mflag
real(wp) :: fx11, fx21, fx12, fx22, fxy1, fxy2

call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag)
call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag)
call dintrv(me%z,z,me%iloz,iz(1),iz(2),mflag)

q1 = (x-me%x(ix(1)))/(me%x(ix(2))-me%x(ix(1)))
q2 = (y-me%y(iy(1)))/(me%y(iy(2))-me%y(iy(1)))
q3 = (z-me%z(iz(1)))/(me%z(iz(2))-me%z(iz(1)))
p1 = one-q1
p2 = one-q2
p3 = one-q3

fx11 = p1*me%f(ix(1),iy(1),iz(1)) + q1*me%f(ix(2),iy(1),iz(1))
fx21 = p1*me%f(ix(1),iy(2),iz(1)) + q1*me%f(ix(2),iy(2),iz(1))
fx12 = p1*me%f(ix(1),iy(1),iz(2)) + q1*me%f(ix(2),iy(1),iz(2))
fx22 = p1*me%f(ix(1),iy(2),iz(2)) + q1*me%f(ix(2),iy(2),iz(2))
fxy1 = p2*( fx11 ) + q2*( fx21 )
fxy2 = p2*( fx12 ) + q2*( fx22 )

fxyz = p3*( fxy1 ) + q3*( fxy2 )

end subroutine interp_3d

The finterp library is released under a permissive BSD-3 license. As far as I can tell, it is unique among publicly-available Fortran codes in providing linear interpolation for up to 6D data sets. There used to be a Fortran 77 library called REGRIDPACK for regridding 1D-4D data sets (it had the option to use linear or cubic interpolation independently for each dimension). Written by John C. Adams at the National Center for Atmospheric Research in 1999, it used to be hosted at the UCAR website, but the link is now dead. You can still find it on the internet though, as part of other projects (for example here, where it is being used via a Python wrapper). But the licensing is unclear to me.

The dearth of easy-to-find, easy-to-use, high-quality open source modern Fortran libraries is a big problem for the long-term future of the language. Fortran users don’t seem to be as interested in promoting their language as users of some of the newer programming languages are. There is a group of us working to change this, but we’ve got a long way to go. For example, Julia, a programming language only four years old, already has a ton of libraries, some of which are just wrappers to Fortran 77 code that no one’s ever bothered to modernize. They even have a yearly conference. An in depth article about this state of affairs is a post for another day.

References

Tagged with: , ,
5 comments on “Multidimensional Linear Interpolation (Part 2)
  1. Lil JD says:

    I note there are six people listed at that “group of us” link and I don’t recognize any of the names

    • Stefano says:

      Dear Lil JD,

      indeed our group is bigger than 6 people. Currently we are 49. However, not all of us have made the group-affiliation public, so other cannot directly see who did not made public her/his affiliation. However, 6 or 50 makes little difference: we should reach the “critical mass” of Julia-like communities with order of 100 or 1000 users to make some difference, but Fortran has a different history with respect all other languages…

      The work of Jacob is making a difference: I using his libraries and I am promoting them to my colleagues, but Jacob alone cannot do all. If you are interested on Fortran, please join us.

      My best regards.

       

  2. Stefano says:

    Dear Jacob,

    I totally agree, as it often happens. I follow Julia and Python (and Nim) with the hope that Fortran will be inspired by them.

    I was not aware about the Julia-con: do you think it could be helpful to have something similar for “modern” Fortran? Maybe, Damian can have the “power” to organize something similar (I know he done such for CAF)…

    Cheers

  3. Jacob says:

    There should certainly be a modern Fortran conference. The closest thing I’ve seen is the BSC Fortran Specialist Group. Maybe there are others (certainly there are HPC conferences). I don’t know who has the “power” to make that happen.

  4. Stefano says:

    Jacob,

    thank you very much, I was not aware about BSC Fortran group.

    I meant “power” as an alias for “visibility, authority, connections…”, Damian is surely a “power Fortraner” 🙂

    Cheers

Leave a Reply

Your email address will not be published. Required fields are marked *

*