Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- MKL Interpolation gives strange results (fortran)

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Erik_T_1

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-09-2015
11:06 PM

157 Views

MKL Interpolation gives strange results (fortran)

With the help of the reference manual for MKL data fitting i have tried to put together a I have a piece of fortran code for interpolation using MKL. It seems to me to give somewhat strange results (I don't know what I am missing here).:

include 'mkl_df.f90' include 'mkl_dfti.f90' program TestMyFortranCode USE MKL_DF_TYPE USE MKL_DF use ISO_C_BINDING TYPE (DF_TASK) dftask integer(c_int),parameter :: NSITES = 1024 integer(c_int) :: f_index(27) real(c_double) :: f(27) real(c_double) :: x(29) integer(c_int) :: nx integer(c_int) :: xhint real(c_double) :: y(29) integer(c_int) :: ny integer(c_int) :: yhint integer(c_int) :: sorder integer(c_int) :: stype real(c_double),target :: site(NSITES) integer(c_int) :: nsite integer(c_int) :: sitehint real(c_double) :: scoeff((size(x)-1)*DF_PP_CUBIC) integer(c_int) :: scoeffhint integer(c_int) :: bc_type integer(c_int) :: ndorder

integer(c_int) :: nder

integer(c_int) :: derorder(3) integer(c_int) :: calc_type integer(c_int) :: method real(c_double),target :: r(size(site)) integer(c_int) :: rhint real(c_double),POINTER :: site_ptr(:) real(c_double),POINTER :: r_ptr(:) integer(c_int) :: error !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! error = 0 site = (/(i*1.0d0, i=0,NSITES-1)/)/(NSITES) ! f_index=(/ (ii,ii=14,40)/) f=10d0**(f_index/10.0d0) x = (/-1.d0, f, 22100.d0/)/22050.d0 nx = size(x) xhint = DF_NO_HINT y = (/0.d0, 5.d0, 1.d0, -3.d0, -2.d0, -6.d0, -9.d0, -10.d0, -14.d0, -14.d0, & -15.d0, -15.d0, -16.d0, -20.d0, -24.d0, -19.d0, -10.d0, -13.d0, -15.d0, -16.d0, & -17.d0, -18.d0, -21.d0, -28.d0, -29.d0, -33.d0, -45.d0, -66.d0, -77.d0/) ny = 1 yhint = DF_NO_HINT nsite = NSITES sitehint = DF_NO_HINT sorder = DF_PP_CUBIC stype = DF_PP_NATURAL scoeffhint = DF_UNIFORM_PARTITION bc_type = DF_BC_FREE_END rhint = DF_MATRIX_STORAGE_ROWS ndorder = 3 derorder = (/1,0,1/) nder = 2 error = dfdNewTask1D(dftask, nx, x, xhint, ny, y, yhint ) error = dfdEditPPSpline1D( dftask, sorder, stype, bc_type, & & scoeff=scoeff, scoeffhint=scoeffhint ) calc_type = DF_PP_SPLINE method = DF_METHOD_STD error = dfdConstruct1D( dftask, calc_type, method ) calc_type = DF_INTERP method = DF_METHOD_PP site_ptr => site r_ptr => r error = dfdInterpolate1D( dftask, calc_type, method, & & nsite, site_ptr, sitehint, ndorder, derorder, r=r_ptr, & & rhint=rhint ) error = dfDeleteTask( dftask ) write (*,*) r_ptr end program TestMyFortranCode

y are the function values of x and site are the interpolation sites (x-values) that I want interpolated values for put in r.

The results for uneven indices looks somewhat as they should, closer to the y-values, but for even indices the results are more random (i.e. -68820.abc or 1013.xyz)

After running the function "dfdNewTask1D" the second value of dftask seems always to be 0 (zero) and the first value seems to get a large random value (I.e. 7424768) (from what I can read out from debug-mode in msvs). Don't know if this is as it should be??

Does anyone know what is going wrong, or if the results are normal, and then what are the results for even indices describing?

(I use MSVS 2013)

1 Solution

Ying_H_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-13-2015
08:13 PM

157 Views

Hi Erik,

we found the reason of the fail.

according to the line derorder = (/1,0,1/), results array should contain function values as long as function second derivative values also. Hence you should allocate result space twice by 1024:

Error line was: real(c_double),target :: r(size(site))

Corrected line: real(c_double),target :: r(size(site)***2**)

if fix this, you will get expected result.

PS: Also we found why we have 10 arguments for dfdinterpolate1d() in user example instead of 12 in **mkl_df.f90** interface file.

Answer: Optional arguments can be omitted.

You can see in mkl\__release_win\mkl\include\mkl_df.f90:

INTERFACE INTEGER FUNCTION dfdinterpolate1d(task, type, method, nsite, &

& site, sitehint, ndorder, dorder, datahint, r, rhint, cell)

USE MKL_DF_TYPE

TYPE(DF_TASK) :: task

INTEGER,INTENT(IN) :: type

INTEGER,INTENT(IN) :: method

INTEGER,INTENT(IN) :: nsite

REAL(KIND=8),DIMENSION(*),INTENT(IN) :: site

** INTEGER,INTENT(IN),OPTIONAL :: sitehint**

** INTEGER,INTENT(IN),OPTIONAL :: ndorder**

** INTEGER, DIMENSION(*),INTENT(IN),OPTIONAL :: dorder**

** REAL(KIND=8),DIMENSION(*),INTENT(IN),OPTIONAL :: datahint**

REAL(KIND=8),DIMENSION(*),INTENT(OUT) :: r

** INTEGER,INTENT(IN),OPTIONAL :: rhint**

** INTEGER,DIMENSION(*),INTENT(OUT),OPTIONAL :: cell**

END FUNCTION

END INTERFACE

that only first 5 + middle 1 arguments are necessary, and 6 ones are optional and can be omitted; and some arguments even **can go in any order** if their name present in function call like in myfunc(arg1,arg1,arg3name=arg3value). That is why

dfdInterpolate1D( dftask, calc_type, method, nsite, site_ptr, sitehint, ndorder, derorder, **rhint=rhint**, r=r_ptr )

is correct call, in spite of order of last 2 arguments.

Best Regards,

Ying

Link Copied

11 Replies

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-10-2015
11:45 AM

157 Views

Which version of MKL are you using? The manual page at https://software.intel.com/en-us/forums/topic/567633 shows the calling sequence as

status = dfdinterpolate1d(task, type, method, nsite, site, sitehint, ndorder, dorder, datahint, r, rhint, cell)

There are no optional arguments and the number of arguments is 12, whereas your call has only 10. Do you have an overloaded Fortran 95 module procedure with the same name? I find this mysterious.

There is a working example, dfdcubicspline_interp.f, that is included with the MKL distribution in the examples/datafittingf/source sub-directory.

Erik_T_1

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-10-2015
09:55 PM

157 Views

Thanks for the reply mecej4. I have the MKL 11.2.3 and I have used the dfdcubicspline_interp.f formulation which only have 10 arguments. (=> I have no overloded module). The formulation in dfdfifthorderspline_interp.f and dfdlookup_interp.f is also only using 10 arguments. None of the interpolation examples makes use of the datahint optional parameter.

According to mkl_df.f90 both the datahint and cell (as well as sitehint, ndorder, dorder and rhint) are optional.

Unfortunately, providing the parameters

real(c_double) :: datahint(4) integer(c_int) :: cell(0) datahint = (/1.,DF_NO_APRIORI_INFO,0.,0./)

makes no significant change to the result (r). The result is very similar, with values for even indices differing quite a bit from the values in y, as without these parameters.

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-11-2015
03:40 AM

157 Views

Here is what I think is happening; if my guess is correct, what is involved is an implementation error that the MKL team needs to correct. I used the Windows Version 15.0.3.208 32-bit compiler, and the version of MKL included with it, namely, 11.2.3.

The source files for the modules MKL_DF_TYPE and MKL_DF contain only interfaces and type definitions. The actual Fortran 95 code needed to redirect the call to the interface name to the actual F77 routines in the libraries is not to be found! Compounding the problem is the choice of the same name(s), e.g., dfdinterpolate1d, for the F95 interface routine and the F77 routine in the library. The linker ends up directing your F95 call to the F77 routine, but the F77 routine expects 12 arguments and knows nothing about how to treat optional arguments. As a result, the last two arguments are pulled from the stack and the values of those are unpredictable because the caller did not put those arguments on the stack.

Here are details of a test that I ran to check those conjectures. I ran the MKL example **dfdcubicspline_interp.f** and the following program, which I know must fail:

program tst call dfdinterpolate1d() end

I compiled both examples with **/Zi /Qmkl**, opened each in the VS debugger and set a breakpoint to the line where **dfdinterpolate1d()** is called. After reaching that breakpoint, I pressed F11, which took me to the entry point **_DFDINTERPOLATE1D** in the MKL DLL. I then looked at the EIP value and the first few lines of the disassembly.

I found that the same MKL entry point was reached in either case. There is a problem: where is the F95 wrapper/glue code?

Erik_T_1

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-11-2015
03:58 AM

157 Views

Once again, thanks a lot mecej4.

What can I do about it at this time? Report the findings somewhere?

As I wrote previously, I tried the function again including the datahint and cell arguments in the routine call, with similar results as if I don't include them. Either I don't understand way the input arguments should be specified or else there is something wrong.. ?

BTW I use Intel(R) Visual Fortran Compiler XE 15.0.4.221 [IA-32]

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-11-2015
04:33 AM

157 Views

Erik T. wrote:

What can I do about it at this time? Report the findings somewhere?

Usually one of the MKL/Intel people will respond in a day or two with a reply in this forum. If you have Premier Support, you can file a formal bug report.

I think that if in your code you make only the F77-style calls (including all arguments as listed in the documentation, even if their values may never be used in the MKL function/subroutine), things may work correctly. If you follow this route, you may need to remove the include lines that make interfaces available, but you need to keep the declarations of types and named constants available for the compiler to use while building your application.

Ying_H_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-11-2015
11:25 PM

157 Views

Hi Erik,

Yes, we will investigate the issue. I just quick run the code and is the below result right or not?

-5270.06285751546 -15.3311452042090 -4371.38789948286

-15.4034790145797 -3472.71294145026 -15.4791246620909

-2574.03798341765 -15.5572251035164 -1675.36302538505

-15.6369232956298 -776.688067352448 -15.7173621952050

121.986890680154 -15.7976847590155 1020.66184871276

-15.8770339438352 1919.33680674536 -15.9545527064378

2818.01176477796 -16.0293944534202 3301.02463352778

-16.1011243940587 3154.54491028549 -16.1698459262363

3008.06518704319 -16.2356987439030 2861.58546380089

-16.2988225410086 2715.10574055859 -16.3593570115032

2568.62601731630 -16.4174418493365 2422.14629407400

-16.4732167484587 2275.66657083170 -16.5268214028195

2129.18684758941 -16.5783955063690 1982.70712434711

-16.6280787530570 1836.22740110481 -16.6760108368335

1689.74767786252 -16.7223314516484 1543.26795462022

-16.7671802914516 1396.78823137792 -16.8106970501931

1250.30850813562 -16.8530214218229 1103.82878489333

-16.8942931002907 957.349061651030 -16.9346517795467

810.869338408733 -16.9742371535406 664.389615166436

-17.0131886368992 533.094098647316 -17.0516296843778

431.247674276086 -17.0896594620255 329.401249904855

-17.1273750981614 227.554825533625 -17.1648737211045

125.708401162395 -17.2022524591742 23.8619767911641

-17.2396084406894 -77.9844475800663 -17.2770387939694

Best Regards,

Ying

Erik_T_1

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-11-2015
11:52 PM

157 Views

Ying H: Yes, the results looks correct, from about index 130 and forward. With diffenent NSITES I get somewhat different numbers but the values of the even indices are similar (too large and shifting between positive and negative sign).

Thank you very much, I am looking forward to a solution!

Ying_H_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-12-2015
07:13 PM

157 Views

Hi Erik,

We can see the problem on seveval machine. It seems works fine on Linux SSE 4.2 machine. we will investigate these and update you if any solution.

Thanks

Ying

Erik_T_1

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-12-2015
11:34 PM

157 Views

Hi Ying,

Thanks. I am running on a lenovo T440 windows 7 machine.

Erik

Ying_H_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-13-2015
08:13 PM

158 Views

Hi Erik,

we found the reason of the fail.

according to the line derorder = (/1,0,1/), results array should contain function values as long as function second derivative values also. Hence you should allocate result space twice by 1024:

Error line was: real(c_double),target :: r(size(site))

Corrected line: real(c_double),target :: r(size(site)***2**)

if fix this, you will get expected result.

PS: Also we found why we have 10 arguments for dfdinterpolate1d() in user example instead of 12 in **mkl_df.f90** interface file.

Answer: Optional arguments can be omitted.

You can see in mkl\__release_win\mkl\include\mkl_df.f90:

INTERFACE INTEGER FUNCTION dfdinterpolate1d(task, type, method, nsite, &

& site, sitehint, ndorder, dorder, datahint, r, rhint, cell)

USE MKL_DF_TYPE

TYPE(DF_TASK) :: task

INTEGER,INTENT(IN) :: type

INTEGER,INTENT(IN) :: method

INTEGER,INTENT(IN) :: nsite

REAL(KIND=8),DIMENSION(*),INTENT(IN) :: site

** INTEGER,INTENT(IN),OPTIONAL :: sitehint**

** INTEGER,INTENT(IN),OPTIONAL :: ndorder**

** INTEGER, DIMENSION(*),INTENT(IN),OPTIONAL :: dorder**

** REAL(KIND=8),DIMENSION(*),INTENT(IN),OPTIONAL :: datahint**

REAL(KIND=8),DIMENSION(*),INTENT(OUT) :: r

** INTEGER,INTENT(IN),OPTIONAL :: rhint**

** INTEGER,DIMENSION(*),INTENT(OUT),OPTIONAL :: cell**

END FUNCTION

END INTERFACE

that only first 5 + middle 1 arguments are necessary, and 6 ones are optional and can be omitted; and some arguments even **can go in any order** if their name present in function call like in myfunc(arg1,arg1,arg3name=arg3value). That is why

dfdInterpolate1D( dftask, calc_type, method, nsite, site_ptr, sitehint, ndorder, derorder, **rhint=rhint**, r=r_ptr )

is correct call, in spite of order of last 2 arguments.

Best Regards,

Ying

Erik_T_1

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-13-2015
09:49 PM

157 Views

Hi Ying,

Thank you for that clarification. (The answer seems so simple now :) )

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

For more complete information about compiler optimizations, see our Optimization Notice.