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
        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…

References

Posted in Programming Tagged with: , ,
One comment on “Multidimensional Linear Interpolation
  1. Stefano says:

    I am tuned …. 🙂

Leave a Reply

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

*