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

MKL interpolation function in Fortran?

BillAvery
Beginner
3,841 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?

0 Kudos
1 Solution
ShanmukhS_Intel
Moderator
3,580 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


View solution in original post

0 Kudos
15 Replies
mecej4
Honored Contributor III
3,814 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)]

You stated that your real problem is more complicated. Please elaborate.

0 Kudos
Ron_Green
Moderator
3,776 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

 

0 Kudos
BillAvery
Beginner
3,758 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.

0 Kudos
David_Billinghurst
New Contributor III
3,739 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)
Distributed under the GNU Public License. See the file COPYING.
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

 

0 Kudos
David_Billinghurst
New Contributor III
3,716 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

 

0 Kudos
BillAvery
Beginner
3,706 Views
0 Kudos
Barbara_P_Intel
Employee
3,750 Views

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

0 Kudos
mecej4
Honored Contributor III
3,741 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).

0 Kudos
Arjen_Markus
Honored Contributor I
3,691 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

See also this thread on Fortran discourse: https://fortran-lang.discourse.group/t/fortran-code-snippets/2150

 

0 Kudos
ShanmukhS_Intel
Moderator
3,581 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


0 Kudos
JohnNichols
Valued Contributor III
3,552 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.  

Here is your example:

Screenshot 2023-03-21 084928.png

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.  

0 Kudos
BillAvery
Beginner
3,525 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 

0 Kudos
jimdempseyatthecove
Honored Contributor III
3,469 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

 

0 Kudos
JohnNichols
Valued Contributor III
3,465 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. 

0 Kudos
ShanmukhS_Intel
Moderator
3,492 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


0 Kudos
Reply