Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
28497 Discussions

## MKL interpolation function in Fortran?

Beginner
2,298 Views

I need to do some simple linear interpolation in a Fortran code. I struggle to find the right subroutine and a simple example on how to do it.

For example, Let x(i) = 1.0, 2.0 ,3.0 ,4.0 ,5.0  and y(i) = 1.0, 4.0, 9.0, 16.0, 25.0

I want to plug in x= 3.75 and return a value for y.

Yes, I know the example above is y=x**2, but my real problem is more complicated.

Can someone point me to the correct mkl subroutine and the syntax for calling it?

1 Solution
Moderator
2,037 Views

Hi William,

Has the information provided helped? Could you please let us know if we could close this case on our end?

Best Regards,

Shanmukh.SS

15 Replies
Honored Contributor III
2,271 Views

Linear interpolation in 1-D is a trivial task, so it is unlikely that anyone would create a library routine for it.

1. Step through the sorted array x and find the index i such that x(i) <= t < x(i+1). If no such index is found, either do extrapolation or declare failure. If the array is large, you may use binary search instead of stepping through the x array sequentially.
2. Calculate the interpolated value s at x = t using s = y(i) + [y(i+1) - y(i)] / [x(i+1) - x(i)] * [t - x(i)]

Moderator
2,233 Views

This is more a question for the MKL Forum but

you can try the MKL function for interpolation.  It uses splines, so it's more complicated than a linear interpolation.

https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/data-fitting-functions/data-fitting-computational-routines/df-interpolate1d-df-interpolateex1d.html

Beginner
2,215 Views

I did not know there was a mkl forum. I will look for it.

Elaborating on the complexity of the problem, the equation I am trying to operate on is

y= x/a + b*(x/c)**n

I need to solve for the slope (dy/dx) as a function of y only. I was unable to generate an analytical solution. So my plan was to generate data points that has y as a function of x, then calculate dy/dx as a function of y, then interpolate that. I thought there might be a mkl function to make life easier. But I can just use the routine suggested by the previous respondent. In any event, thanks for the help.

New Contributor III
2,196 Views

Using interpolation to calculate derivatives amplifies any interpolation errors.  You should use a method that guarantees continuous derivatives and discourages oscillation from overfitting.  Cubic splines can be a good choice.

Fortunately your expression y= x/a + b*(x/c)**n has an analytic derivative.  (I cheated and used the symbolic algebra package maxima.)  dy/dx is the final expression below.  Checked it with a few random values.

Maxima 5.46.0 https://maxima.sourceforge.io
using Lisp CLISP 2.49+ (2010-07-17)
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) ex:x/a + b*(x/c)^n;
x n   x
(%o1)                             b (-)  + -
c     a
(%i2) diff(ex,x);
x n
b n (-)
c     1
(%o2)                            -------- + -
x       a

New Contributor III
2,173 Views

maxima can generate the Fortran expression, too.  Trivial in this case, but useful for longer expressions.  The fortran() function generates fixed format.  There is also an f90() function that generates free-format.

(%i1) y:x/a + b*(x/c)^n;
x n   x
(%o1)                             b (-)  + -
c     a
(%i2) dydx:diff(y,x);
x n
b n (-)
c     1
(%o2)                            -------- + -
x       a
(%i3) fortran(dydx);
(b*n*(x/c)**n)/x+1/a
(%o3)                                done

Beginner
2,163 Views

way cool. thanks.

Moderator
2,207 Views

Here's the link to the MKL Forum, if you need more help, although there are knowledgeable MKL users who hang out here.

Honored Contributor III
2,198 Views

Unless the values of a, b, c and n, and the range of interest in x conspire to make things simple, the expression that you revealed is not amenable to linear interpolation. If the nonlinearity is high, linear interpolation can be quite inaccurate.

The variable p = dy/dx may not be monotonic in x, so you have to be in possession of a rule to select one value based on some external reason such as physics, etc. Similarly, if the range of x contains 0 and n is not an integer, more study is needed.

Much can be learned by making a phase-plane plot of p versus y and inspecting the plot. This should be done before selecting a numerical routine (such as one of the spline interpolation routines available in MKL).

Honored Contributor I
2,148 Views

You could also use automatic differentiation. It is a fairly simple technique - see for instance https://sourceforge.net/p/flibs/svncode/HEAD/tree/trunk/src/computing/automdiff.f90. But an even simpler technique is to use complex numbers. Here is a small program I wrote to experiment with this technique:

! derivative.f90 --
!     Use the complex step method to determine the derivative of a real function
!
!     Method described by Nicholas Higham
!
!     With simple expressions the result is quite intriguingly accurate. Note though the
!     last example - there you get a considerable fluctuation even with h = 1.0e-20 to 1.0e-30
!
program derivative
implicit none

real(kind=kind(1.0d0))    :: h
complex(kind=kind(1.0d0)) :: x, y
integer :: i

do i = 1,10
h = 0.1d0 ** i
x = cmplx(1.0d0,h,kind=kind(1.0d0))

write(*,*) imag( sin(x) ) / h,  cos(real(x)), (imag( sin(x) ) / h - cos(real(x))) / cos(real(x))
write(*,*) imag( cos(x) ) / h, -sin(real(x)), (imag( cos(x) ) / h + sin(real(x))) / (-sin(real(x)))
enddo
do i = 1,60
h = 0.1d0 ** (3*i)
x = cmplx(1.0d0,h,kind=kind(1.0d0))

y = x + x**2 + x**3

write(*,*) imag( y ) / h, h
enddo
do i = 1,60
h = 0.1d0 ** (3*i)
x = cmplx(1.0d0,h,kind=kind(1.0d0))

y = cos(x)/(1.0d0 - sin(x))

write(*,*) imag( y ) / h, h
enddo
end program derivative

Moderator
2,038 Views

Hi William,

Has the information provided helped? Could you please let us know if we could close this case on our end?

Best Regards,

Shanmukh.SS

Valued Contributor III
2,009 Views

One quick way to work out what function may fit, before you get stuck into Fortran is to look at the data in EXCEL.  It has many fine interpolation functions.

Then you use Fortran.   of course if the data is complex, like the COVID death data, you use a FFT routine to look for the waves in the data.  A French group of Dr missed a peak in one of their published data sets because they did not use FFT.

It is an interesting field in stats, but you have to know what you are going to program if you want to save a lot of time.

Good luck.  PS You got excellent help in the posts above.

Beginner
1,982 Views

Thanks. I actually tried that approach and got close, but not close enough as the tail of the data trailed off. Could just not find the right functional form

Honored Contributor III
1,926 Views

>> not close enough as the tail of the data trailed off.

That may be an artifact of the signal (or simulation) not stabilizing and/or your sampling rate being too slow to catch the higher frequencies that may be present at the time on an event (equivalent to simulation integration step size being to long).

Jim Dempsey

Valued Contributor III
1,922 Views

Give a data sample, you have, excluding myself, some of the best Fortran programmers and excellent applied math people.  I just like to hang out with the bright people at the nerds table.

Not all data fits one of the few patterns we understand, most Gaussian data is not really Gaussian it is just close enough that we accept the difference and get on with fun programming.

There is a book - How round is my circle, I was thinking about it in the library the other day as I watched a lady cut our Easter eggs, I suddenly thought - do a linear transformation and then linear regression.  For @mecej4 - solving the implicit problem in the title of the book.

Or as my pure math teacher once said if you look at some obscure figure in the 23rd dimension through the 22nd dimension you will see a donut.  If that guy had said the moon was made of blue cheese I would have accepted it directly, I failed his exam five times and finally gave up only to be told the exam was out of his textbook which was not on the syllabus.  Never do it alone.  As I explained to my 16 yr old who failed a French exam at school when she thought I would be upset, Darling I failed exams, it is not the knocking down that is important it is the getting back up.

And I have a theory that @mecej4 is really the Fonz.

Moderator
1,949 Views

Hi William,

Glad to know that your issue is resolved. If you need any additional information, please post a new question as this thread will no longer be monitored by Intel and further discussion within the forum is considered community only.

Have a great day!

Best Regards,

Shanmukh.SS