Degenerate Conic

Algorithms • Modern Fortran Programming • Orbital Mechanics

Feb 28, 2016

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