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 elseif ( x==ent(locc) ) then ishift = (locc-lmin)*istep do k = 1 , knots index(k) = index(k) + ishift enddo istep = istep*ndim cycle main else loca = locc endif if ( locb-loca<=1 ) exit enddo 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 endif ishift = istep if ( x-ent(lmin+1)==0.0_wp ) then do k = 1 , knots index(k) = index(k) + ishift enddo istep = istep*ndim cycle main endif ishift = 0 eta = h/(ent(lmin+1)-ent(lmin)) endif 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) enddo knots = 2*knots istep = istep*ndim enddo main fint = 0.0_wp do k = 1 , knots i = index(k) fint = fint + weight(k)*table(i) enddo 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…