topic Here is what I think is in Intel® oneAPI Math Kernel Library

MKL Interpolation gives strange results (fortran)

Erik_T_1 — Mon, 10 Aug 2015 06:06:53 GMT

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)

Which version of MKL are you

mecej4 — Mon, 10 Aug 2015 18:45:13 GMT

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.

Thanks for the reply mecej4.

Erik_T_1 — Tue, 11 Aug 2015 04:55:00 GMT

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.

Here is what I think is

mecej4 — Tue, 11 Aug 2015 10:40:00 GMT

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?

Once again, thanks a lot

Erik_T_1 — Tue, 11 Aug 2015 10:58:00 GMT

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]

Quote:Erik T. wrote:

mecej4 — Tue, 11 Aug 2015 11:33:00 GMT

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.

Hi Erik,

Ying_H_Intel — Wed, 12 Aug 2015 06:25:52 GMT

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

Ying H: Yes, the results

Erik_T_1 — Wed, 12 Aug 2015 06:52:30 GMT

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!

Hi Erik,

Ying_H_Intel — Thu, 13 Aug 2015 02:13:00 GMT

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

Hi Ying,

Erik_T_1 — Thu, 13 Aug 2015 06:34:21 GMT

Hi Ying,

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

Erik

Hi Erik,

Ying_H_Intel — Fri, 14 Aug 2015 03:13:21 GMT

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

Hi Ying,

Erik_T_1 — Fri, 14 Aug 2015 04:49:20 GMT

Hi Ying,

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