Multidimensional Linear Interpolation
Some years ago, I needed a Fortran routine to do linear interpolation of a multidimensional (up to 6D) data set. Not wanting to reinvent the wheel, I ended up using a routine called FINT from CERNLIB (documentation here). It is written in Fortran 66, and has some hard-coded limits on the number of dimensions that can be used, but these are not fundamental to the algorithm. With a few minor changes, the routine can be used for interpolation of any number of dimensions. I present below a modern Fortran version that removes the limitations. Note that the original code is GPL licensed, thus so is this modification.
function fint(narg,arg,nent,ent,table)
use iso_fortran_env, only: wp => real64
implicit none
real(wp) :: fint
integer,intent(in) :: narg
integer,dimension(narg),intent(in) :: nent
real(wp),dimension(narg),intent(in) :: arg
real(wp),dimension(:),intent(in) :: ent
real(wp),dimension(:),intent(in) :: table
integer,dimension(2**narg) :: index
real(wp),dimension(2**narg) :: weight
logical :: mflag,rflag
real(wp) :: eta,h,x
integer :: i,ishift,istep,k,knots,lgfile,&
lmax,lmin,loca,locb,locc,ndim
!some error checks:
if (size(ent)/=sum(nent)) &
error stop 'size of ent is incorrect.'
if (size(table)/=product(nent)) &
error stop 'size of table is incorrect.'
if (narg<1) &
error stop 'invalid value of narg.'
lmax = 0
istep = 1
knots = 1
index(1) = 1
weight(1) = 1.0_wp
main: do i = 1 , narg
x = arg(i)
ndim = nent(i)
loca = lmax
lmin = lmax + 1
lmax = lmax + ndim
if ( ndim>2 ) then
locb = lmax + 1
do
locc = (loca+locb)/2
if ( x<ent(locc) ) then
locb = locc
else if ( x==ent(locc) ) then
ishift = (locc-lmin)*istep
do k = 1 , knots
index(k) = index(k) + ishift
end do
istep = istep*ndim
cycle main
else
loca = locc
end if
if ( locb-loca<=1 ) exit
end do
loca = min(max(loca,lmin),lmax-1)
ishift = (loca-lmin)*istep
eta = (x-ent(loca))/(ent(loca+1)-ent(loca))
else
if ( ndim==1 ) cycle main
h = x - ent(lmin)
if ( h==0.0_wp ) then
istep = istep*ndim
cycle main
end if
ishift = istep
if ( x-ent(lmin+1)==0.0_wp ) then
do k = 1 , knots
index(k) = index(k) + ishift
end do
istep = istep*ndim
cycle main
end if
ishift = 0
eta = h/(ent(lmin+1)-ent(lmin))
end if
do k = 1 , knots
index(k) = index(k) + ishift
index(k+knots) = index(k) + istep
weight(k+knots) = weight(k)*eta
weight(k) = weight(k) - weight(k+knots)
end do
knots = 2*knots
istep = istep*ndim
end do main
fint = 0.0_wp
do k = 1 , knots
i = index(k)
fint = fint + weight(k)*table(i)
end do
end function fint
As I recall, I wrote wrappers to this function in order to be able to input the data in a more reasonable format. For example, for a five dimensional data set it is more natural to input the five rank-1 abscissa vectors and the rank-5 ordinate matrix. The disadvantage of this routine is that you have to stuff the ordinate matrix into a vector (as well as stuffing all the abscissa vectors into a single vector). This could also be an advantage in some cases, since Fortran arrays have a maximum rank. For Fortran 90, it was 7, while Fortran 2008 expanded this to 15 (the current Intel Fortran Compiler allows up to 31 as a non-standard extension). If you ever, for example, needed to interpolate 32-dimensional data (assuming you had enough RAM), you could do it with FINT without having to declare any rank-32 matrices. In any event, an example 5D wrapper for FINT is given here:
function fint_5d_wrapper(x,y,z,q,r,&
xvec,yvec,zvec,qvec,&
rvec,fmat) result(f)
use iso_fortran_env, only: wp => real64
implicit none
integer,parameter :: narg = 5
real(wp) :: f
real(wp),intent(in) :: x,y,z,q,r
real(wp),dimension(:),intent(in) :: xvec
real(wp),dimension(:),intent(in) :: yvec
real(wp),dimension(:),intent(in) :: zvec
real(wp),dimension(:),intent(in) :: qvec
real(wp),dimension(:),intent(in) :: rvec
real(wp),dimension(:,:,:,:,:),intent(in) :: fmat
integer,dimension(narg) :: nent
nent(1) = size(xvec)
nent(2) = size(yvec)
nent(3) = size(zvec)
nent(4) = size(qvec)
nent(5) = size(rvec)
f = fint(narg,[x,y,z,q,r],nent,&
[xvec,yvec,zvec,qvec,rvec],&
reshape(fmat,[product(nent)]))
end function fint_5d_wrapper
This works well enough, but note that it isn't very computationally efficient since we are creating temporary copies of the potentially large data arrays (including having to reshape a rank-5 matrix) every time we do an interpolation. While the logic of FINT is a bit hard to fathom, the actual equations are fairly simple (see here, for example). What we really need now is an object-oriented modern Fortran library for multidimensional linear interpolation. Stay tuned...
References
- CERNLIB – CERN Program Library Short write-ups, June 1996.
- G. Hoffmann, Multidimensional Linear Interpolation