# Fortran, making a more efficient bilinear interpolation

I'm trying to write an efficient bilinear (2D)-interpolation, after reading some recipes, as a fortran-mex for Matlab that is used extensively throughout a long algorithm of solar image processing, and therefore is one of my main bottlenecks. My 1st attempt work but is slow. See code below. Any advice that can speed up the code is most welcome. I'm using gcc version 4.3.6. I cannot use a later one as it the only one fully compatible with my version of Matlab. Any later version prevented the use of some embedded matlab-mex functions necessary for interfacing / displaying error messages in the matlab prompts. My mex-setup can compile Fortran 90 and 95 (so far...)

I already have a C++ mex equivalent used so far and was hoping that a well-written Fortran code would be faster, i'm using the C++ as reference for the speedup that Fortran can or cannot offer.

SUBROUTINE L2DINTERPOL3(IntIm,Image,x,y,NPts,M,N)

INTEGER*8        :: NPts,M,N,jj
INTEGER*8        :: x1(NPts),y1(NPts),x2(NPts),y2(NPts)
DOUBLE PRECISION :: IntIm(NPts)
DOUBLE PRECISION :: Image(M,N)
DOUBLE PRECISION :: f1,f2
DOUBLE PRECISION :: x(NPts),y(NPts)
DOUBLE PRECISION :: wx1(NPts),wx2(NPts),wy1(NPts),wy2(NPts)

c Take the nearest integers around the input.

x1 = FLOOR(x)
x2 = CEILING(x)
y1 = FLOOR(y)
y2 = CEILING(y)

wx1 = (x2-x)/(x2-x1)
wx2 = (x-x1)/(x2-x1)
wy1 = (y2-y)/(y2-y1)
wy2 = (y-y1)/(y2-y1)

c Whenever the input are already integers, weights go to infinity,
c so set each pair to 1 and 0.

WHERE(x1 .EQ. x2)
wx1 = 1.
wx2 = 0.
END WHERE

WHERE (y1 .EQ. y2)
wy1 = 1.
wy2 = 0.
END WHERE

c Main calculation, loop over each element of the input.

DO 10 jj=1,NPts

f1 = wx1(jj)*Image(y1(jj),x1(jj)) + wx2(jj)*Image(y1(jj),x2(jj))
f2 = wx1(jj)*Image(y2(jj),x1(jj)) + wx2(jj)*Image(y2(jj),x2(jj))

IntIm(jj) = wy1(jj)*f1 + wy2(jj)*f2

10   CONTINUE

RETURN
END

[EDIT] As i only have regular grid, x2-x1 = 1 or 0, i've made another version that does not require division. This got rid of the where constructs. In addition, it seems that replace the DO loop with a "FORALL" loop improved by about 15% on 2 cpus. My understanding is that FORALL uses multithread when present, and each computation at each jj is independent from the other at jj+k or jj-k which makes the FORALL statement sensible. Doesn't it ?

SUBROUTINE L2DINTERPOL(IntIm,Image,x,y,NPts,M,N)
implicit none

mwSize, PARAMETER             :: dp = kind(0.d0) ! Double precision
mwSize                        :: NPts, M,N       ! Input
REAL(dp),DIMENSION(Npts)      :: x,y             ! Input
REAL(dp),DIMENSION(M,N)       :: Image           ! Input
mwSize                        :: jj
mwSize, DIMENSION(NPts)       :: x1,y1,x2,y2
REAL(dp),DIMENSION(Npts)      :: wx1,wx2,wy1,wy2
REAL(dp),DIMENSION(Npts)      :: IntIm           ! Output

x1 = FLOOR(x)
x2 = CEILING(x)
y1 = FLOOR(y)
y2 = CEILING(y)

wx1 = x2-x
wy1 = y2-y

wx2 = 1 - wx1
wy2 = 1 - wy1

FORALL( jj=1 : NPts)

IntIm(jj) = wy1(jj)*(wx1(jj)*Image(y1(jj),x1(jj))+wx2(jj)*
$Image(y1(jj),x2(jj))) + wy2(jj)*(wx1(jj)*$        Image(y2(jj),x1(jj))+wx2(jj)*Image(y2(jj),x2(jj)))

END FORALL

RETURN
END
• Can you post a full program, that can be compiled and it returns a timing? How many points (NPts) do you have? What are the dimensions (M, N) of the image? So that we can make sure that we benchmark your actual problem. I wrote a trilinear interpolation (in 3d) here: github.com/certik/hfsolver/blob/master/src/interp3d.f90, so you can compare the speed against it. How are the points in x, y distributed? Is it a regular mesh? If so, then you can exploit this feature to possibly speed things up. Commented Dec 2, 2013 at 17:57
• Npts can be of the order 10^5, 10^6 tops. This is a function for Matlab mex files only. The goal here is to make it faster than the c++ version that i also have as Matlab mex c++ file. So i need to stick to that for direct comparison. But i can put the mex gateway if you're experienced with it. I'm updating my original question with a more compact code, as i indeed have a regular grid. Commented Dec 3, 2013 at 9:17
• (M, N) can varies from (512, 512) up to (1000,1000), they can be rectangular, length of each dimension can be odd or even. Commented Dec 3, 2013 at 9:32
• forall doesn't use threads, as far as I know. But you should be able to use openmp to parallelize the loop. The next improvement would be to use the fact that you have a 2d regular grid, so instead of using x1, y1 1D array, you can use a single 2D array, or possibly get rid of arrays completely and just use increments (i.e. it's regular, so it's easy to calculate the coordinates). Typically what I do is I try to write the algorithm in couple different ways, optimize each and see which one is the fastest. Commented Dec 3, 2013 at 19:32
• I read again about forall. it is parallelized by some compilers. I guess GCC is not one of them, it's not mentioned anywhere. mun.ca/hpc/hpf_pse/manual/hpf0031.htm As for your suggestion, i don't see what you mean i tried a 1D array with converting (x1,y1) coords to its 1D equivalent. I don't see how i'd get rid of arrays. I have to extract the neighboring points for each interpolated point, that will still make me deal with vector indexing, at least ? Can you put a simple example ? I'm not following. Commented Dec 4, 2013 at 12:06

Your solution looks pretty optimal to me (assuming it works correctly, which I haven't checked). I don't think there is a lot that you can optimize -- switching language, compiler, optimization flags may give you 20% speedup, but nothing that will change the fact that this function is your bottleneck. If you need to make your code faster, then you will need to change the algorithm that calls this function so often.

• I've updated the code a bit. Getting rid of the WHERE construct. I've also converted the Matlab algorithm that calls this function to another Fortran mex file. This did not improve computing time. I'm not experienced with the compiler options. I've changed the optimisation from -O to -O3, it seemed to improve a bit, but ~10 %. Not bad, but i don't know what i'm doing when i tweak this. I read a few things of what difference there's between -O, -O2, -03,... I also tried a 1D access of the arrays, instead of 2D. That did not improve speed at all. Commented Dec 3, 2013 at 9:39
• Use these compiler options for production runs (benchmarking): fortran90.org/src/… Commented Dec 3, 2013 at 19:25

Did you try to sort the $\tt x$ array before invoking L2DINTERPOL3 ? I mean only the $\tt x$ array.

What I am thinking about is based on the fact that two following arrays are neighbor of each other with $M$- length of increment (i.e. column-major in Fortran array model),

Image( * , x1(jj) ) and Image( * , x2(jj) )

due to ${\tt x2(jj)} = {\tt x1(jj)} +1$, since ${\tt x1(jj)} = [x_{\rm jj}]$ and ${\tt x2(jj)} = [x_{\rm jj}] +1$. Sorting $\tt x$ before invoking L2DINTERPOL3 may help $\tt x1(jj+1)$ to be near $\tt x1(jj)$ when we proceed $\tt jj = jj+1$ in the loop. Hence, it may optimize access to the storage of $\tt Image(:,:)$.

This idea is simply described as follows:

call SORT( NPts, x, iorder )

do j = 1, NPts
x(j) = x( iorder(j) )
y(j) = y( iorder(j) )
enddo

call L2DINTERPOL3(IntIm,Image,x,y,NPts,M,N)

do j = 1, NPts
IntIm( iorder(j) ) = IntIm(j)
enddo

where the last loop is to recover $\tt IntIm(:)$ based on $\tt iorder(:)$. Can you implement this idea and let me know its effect?