MKL Interpolation gives strange results (fortran)

[Erik_T_1]

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)

[mecej4]

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]

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]

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]

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]

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]

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]

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]

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]

Hi Ying,

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

Erik

[Ying_H_Intel · Accepted solution]

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]

Hi Ying,

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