Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6956 Discussions

Using MKL DataFitting cubic spline routines

geerits
New Contributor I
2,396 Views

Dear reader,

I'm in the process of evaluating whether it is a good idea to switch from IMSL to MKL. Since I have been assisted quite well with my previous queries regarding MKL routines, I will give it another try with my next problem:

 

I'm trying to replace the IMSL routine DCSIEZ(..) which does the most basic cubic spline interpolation using the "not-a-knot" boundary condition. I'm trying to mimic this routine using the MKL DataFitting routines. To this purpose I wrote a little wrapper routine (MKL_C_SPLINE(..)) that exposes exactly the same input/output as my DCSIEZ() routine. In this routine I use a modified version of your cubic spline-based C example as presented in the MKL reference manual. You find the routine at the bottom of this mail. When I replace DCSIEZ() by MKL_C_SPLINE(), the code still compiles and links without any problems. However, at run-time the code crashes. I have a strong impression that this might have something to do with a possible data type conflict. I'm using the standard C data types (int, float and double) but your DataFitting example uses generic types (e.g., MKL_INT). Frankly, I do not know/understand how this works exactly, even after reading some of the documentation. Can you please have a look at my wrapper routine below and give some explanations as to what I'm doing wrong??? Your help is very much appreciated!

 

Best rgrds,

Tim.

 

#include "mkl.h"

/* This routine replaces IMSL cubic spline routine DCSIEZ()
This (wrapper) routine calls the appropriate MKL DataFitting
routines as explained in the MKL reference manual.
*/
#define SPLINE_ORDER DF_PP_CUBIC /* A cubic spline to construct */

void MKL_C_SPLINE(int *N, double *F_TMP, double *S_TMP, int *N_OUT, double *F_NEW, double *S_NEW)
{
int NX, NSITE;
int status; /* Status of a Data Fitting operation */

double *scoeff;

DFTaskPtr task; /* Data Fitting operations are task based */

/* Parameters describing the partition */
MKL_INT nx; /* The size of partition x */
/*double x[NX]; Partition x; This refers to the input array F_TMP of size N */
MKL_INT xhint; /* Additional information about the structure of breakpoints */

/* Parameters describing the function */
MKL_INT ny; /* Function dimension */
/*double y[NX]; Function values at the breakpoints; This refers to the input array S_TMP of size N */
MKL_INT yhint; /* Additional information about the function */

/* Parameters describing the spline */
MKL_INT s_order; /* Spline order */
MKL_INT s_type; /* Spline type */
MKL_INT ic_type; /* Type of internal conditions */
double* ic; /* Array of internal conditions */
MKL_INT bc_type; /* Type of boundary conditions */
double* bc; /* Array of boundary conditions */

 

/*double scoeff[(NX - 1)* SPLINE_ORDER]; Array of spline coefficients */
MKL_INT scoeffhint; /* Additional information about the coefficients */

/* Parameters describing interpolation computations */
MKL_INT nsite; /* Number of interpolation sites */
/*double site[NSITE]; Array of interpolation sites */
MKL_INT sitehint; /* Additional information about the structure of
interpolation sites */

MKL_INT ndorder, dorder; /* Parameters defining the type of interpolation */

double* datahint; /* Additional information on partition and interpolation sites */

/* double r[NSITE]; Array of interpolation results */
MKL_INT rhint; /* Additional information on the structure of the results */
MKL_INT* cell; /* Array of cell indices */

NX = *N;
NSITE = *N_OUT;

scoeff = (double*)calloc((NX - 1)* SPLINE_ORDER, sizeof(double));

/* Initialize the partition (The partition is input F_TMP) */
nx = NX;
xhint = DF_NON_UNIFORM_PARTITION; /* The partition is non-uniform. */

/* Initialize the function */
ny = 1; /* The function is scalar. */

/* Set function values (The function is input S_TMP) */
yhint = DF_NO_HINT; /* No additional information about the function is provided. */

/* Create a Data Fitting task */
status = dfdNewTask1D(&task, nx, F_TMP, xhint, ny, S_TMP, yhint);

/* Check the Data Fitting operation status (TO BE IMPLEMENTED!!) */

 

/* Initialize spline parameters */
s_order = DF_PP_CUBIC; /* Spline is of the fourth order (cubic spline). */
s_type = DF_PP_BESSEL; /* Spline is of the Bessel cubic type. */

/* Define internal conditions for cubic spline construction (none in this example) */
ic_type = DF_NO_IC;
ic = NULL;

/* Use not-a-knot boundary conditions. In this case, the is first and the last
interior breakpoints are inactive, no additional values are provided. */
bc_type = DF_BC_NOT_A_KNOT;
bc = NULL;
scoeffhint = DF_NO_HINT; /* No additional information about the spline. */

/* Set spline parameters in the Data Fitting task */
status = dfdEditPPSpline1D(task, s_order, s_type, bc_type, bc, ic_type,
ic, scoeff, scoeffhint);

/* Check the Data Fitting operation status (TO BE IMPLEMENTED!!) */

/* Use a standard method to construct a cubic Bessel spline: */
/* Pi(x) = ci,0 + ci,1(x - xi) + ci,2(x - xi)2 + ci,3(x - xi)3, */
/* The library packs spline coefficients to array scoeff: */
/* scoeff[4*i+0] = ci,0, scoef[4*i+1] = ci,1, */
/* scoeff[4*i+2] = ci,2, scoef[4*i+1] = ci,3, */
/* i=0,...,N-2 */

status = dfdConstruct1D(task, DF_PP_SPLINE, DF_METHOD_STD);

/* Check the Data Fitting operation status (TO BE IMPLEMENTED!!) */

/* Initialize interpolation parameters */
nsite = NSITE;

/* Set site values, i.e., values at which spline needs to be evaluated. These are contained in F_NEW (input of this routine) */

sitehint = DF_NON_UNIFORM_PARTITION; /* Partition of sites is non-uniform */

/* Request to compute spline values */
ndorder = 1;
dorder = 1;
datahint = DF_NO_APRIORI_INFO; /* No additional information about breakpoints or
sites is provided. */
rhint = DF_MATRIX_STORAGE_ROWS; /* The library packs interpolation results
in row-major format. */
cell = NULL; /* Cell indices are not required. */

/* Solve interpolation problem using the default method: compute the spline values
at the points site(i), i=0,..., nsite-1 and place the results to array r */
status = dfdInterpolate1D(task, DF_INTERP, DF_METHOD_STD, nsite, F_NEW,
sitehint, ndorder, &dorder, datahint, S_NEW,
rhint, cell);

/* Check Data Fitting operation status */

/* De-allocate Data Fitting task resources */
status = dfDeleteTask(&task);

/* Check Data Fitting operation status(TO BE IMPLEMENTED!!) */

free(scoeff);


}

 

0 Kudos
1 Solution
Gennady_F_Intel
Moderator
2,142 Views

This issue has been resolved and we will no longer respond to this thread. If you require additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only. 



View solution in original post

17 Replies
Gennady_F_Intel
Moderator
2,375 Views

Tim,

at the very first glance, this code looks correct. Could You give us the whole reproducer code which e could check and investigate on our end?


0 Kudos
geerits
New Contributor I
2,358 Views

Hi Gennady,

Thanks for your quick reply and happy to hear that my code seems to be correct. That by itself is already a hint :-).

The call to this routine happens in a medium-large size code with several calls to different IMSL routines. That is, I currently do not have a simple 'reproducer' code for you. As a very first step I just tried to replace DCSIEZ() by MKL_C_SPLINE() to see what happens. Unfortunately the code crashes. I will now compile it in sequential mode (it was running in parallel mode using OpenMP) to debug it. I hope this will reveal some clues as to the problem(s). I will isolate the problem step by step and will get back with you if new questions pop up. Hope that's ok?

 

Best regards,

Tim.

0 Kudos
geerits
New Contributor I
2,349 Views

Hi Gennady,

As mentioned in my previous reply I have started debugging my code in sequential mode. The code crashes during execution at the very beginning, i.e., during the task initiation in  MKL_C_SPLINE():

 

status = dfdNewTask1D(&task, nx, F_TMP, xhint, ny, S_TMP, yhint);

 

 I have double checked the type, size and content of F_TMP and S_TMP. Both have size nx (11 in this case) and all contain valid double values. Any useful comments/hints from your side??

 

With kind regards,

Tim.

0 Kudos
Gennady_F_Intel
Moderator
2,344 Views

Tim,

Could you show how do you compile and link this example?

-Gennady



0 Kudos
geerits
New Contributor I
2,335 Views

Hi Gennady,

 

I use Visual Studio but with the Intel compiler (Version 19). Below you find the compiler command line (copied from the C/C++-->Command Line tab):

 

/MP /GS /Qopenmp_stubs /GA /W3 /Gy /Zc:wchar_t /I"C:\MATLAB\R2021a\extern\include" /I"C:\Program Files (x86)\VNI\imsl\fnl701\Intel64\include\dll" /ZI /Od /Fd"x64\Release\vc140.pdb" /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MATLAB_MEX_FILE" /D "_WINDLL" /D "_UNICODE" /D "UNICODE" /Qipo /Zc:forScope /Oi /MDd /Fa"x64\Release\" /EHsc /Fo"x64\Release\" /Qprof-dir "x64\Release\" /Fp"x64\Release\WV_DAMA_SLOW_INV_OpenMP_MEX64.pch"

 

The linker command line looks as follows:

 

/MP /GS /Qopenmp_stubs /GA /W3 /Gy /Zc:wchar_t /I"C:\MATLAB\R2021a\extern\include" /I"C:\Program Files (x86)\VNI\imsl\fnl701\Intel64\include\dll" /ZI /Od /Fd"x64\Release\vc140.pdb" /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MATLAB_MEX_FILE" /D "_WINDLL" /D "_UNICODE" /D "UNICODE" /Qipo /Zc:forScope /Oi /MDd /Fa"x64\Release\" /EHsc /Fo"x64\Release\" /Qprof-dir "x64\Release\" /Fp"x64\Release\WV_DAMA_SLOW_INV_OpenMP_MEX64.pch"

 

Please let me know in case you need additional info. Looking forward hearing from you.

 

With kind regards,

Tim.

 

 

0 Kudos
Gennady_F_Intel
Moderator
2,312 Views

Hi Tim,

We don't see which MKL's libs you use while linking this example.

Please take a look at the MKL LInker Adviser ( https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl-link-line-advisor.html#gs.m4f0om ) to see the recommended list of libraries you need to use.

For example, for win64 and static sequential case, You need to link against these mkl's libs:

 mkl_intel_lp64.lib mkl_sequential.lib mkl_core.lib


-Gennady



0 Kudos
geerits
New Contributor I
2,301 Views

Hi Gennady,

Ok, I see that the command line, i.e., C/C++-->Command Line tab under the "properties" tab is not the full story. Sorry for that! I have attached a copy of the compiler/linker report (from VS) in pdf format. It is very long but since I do not need what the problem is I thought I better give you the whole report :-). As you can see, the three libraries you are referring to are included (i.e., their dll versions).

 

Your help is very much appreciated!

 

Best regards,

Tim.

 

0 Kudos
mecej4
Honored Contributor III
2,276 Views

Tim,

The attached file is a working version of the MKL example code that you have been using. It fits a cubic spline with not-a-knot end conditions to 11 data points (x, y) with x = 0, 0.1, ..., 0.9, 1.0 and y = sin(x), and interpolates at 0.05, 0.15, ..., 0.95. The output is:

T:\LANG\MKL>splfit
New task status = 0
EditPPSpline status = 0
Construct status = 0
Interpolate task status = 0
 0.05 0.050041 0.049979
 0.15 0.149438 0.149438
 0.25 0.247403 0.247404
 0.35 0.342897 0.342898
 0.45 0.434965 0.434966
 0.55 0.522686 0.522687
 0.65 0.605185 0.605186
 0.75 0.681637 0.681639
 0.85 0.751279 0.751280
 0.95 0.813377 0.813416
Delete task status = 0
0 Kudos
mecej4
Honored Contributor III
2,252 Views

Tim,

Here is a Fortran driver for the C subroutine that you listed in the first post of this thread. That C code has one bug, which will become evident if you print out the value of the  status variable after each MKL routine returns. The returned value for the dfdInterpolate1D routine is -1002, which means "unsupported method". As documented in the MKL reference for that routine, the supported method is denoted DF_METHOD_PP, in contrast to your choice of DF_METHOD_STD. Fix that, and your code should work properly.

 

 

program tspline
implicit none
integer :: i, N,Ni, Nout
double precision :: x(11),y(11),xi(10),yi(10)
do i=1,11
   x(i) = (i-1)*0.1d0
   y(i) = sin(x(i))
end do
do i = 1,10
   xi(i) = i*0.1d0-0.05d0
end do
N = 11
NI = 10
call MKL_C_SPLINE(N, x, y, Ni, xi, yi)
print 10,(i,xi(i),yi(i),sin(xi(i)),i=1,10)
10 format(i3,3ES12.5)
end program

 

 

0 Kudos
geerits
New Contributor I
2,237 Views

Thanks for all your efforts! However, it still doesn't solve my problem. Although you correctly noted the faulty method I used in "dfdInterpolate1D" (i.e., DF_METHOD_STD) I would like to recall that my code crashes on:

 

status = dfdNewTask1D(&task, nx, F_TMP, xhint, ny, S_TMP, yhint);

 

in MKL_C_SPLINE() and not on "dfdInterpolate1D". It doesn't come that far. So I cannot check any status, not even that of "dfdNewTask1D". Despite of that I have changed the method in "dfdInterpolate1D" from DF_METHOD_STD to DF_METHOD_PP, as suggested by you and reran the code in DEBUG mode. The result was exactly the same as prior to the change. The code crashes during execution of:

 

status = dfdNewTask1D(&task, nx, F_TMP, xhint, ny, S_TMP, yhint);

 

Did you have a look at the compiler/linker report I attached in an earlier mail? I attached it again to this mail. Does that show anything suspicious to you?

 

PS: By the way, in the documentation "C-code example of cubic spline-based interpolation", the method used in "dfdInterpolate1D" is DF_METHOD_STD. This is inconsistent with the documentation of "dfdInterpolate1D".

 

With kind regards,

Tim.

 

 

 

0 Kudos
mecej4
Honored Contributor III
2,227 Views

The linker report is very complex and provides an incomplete story regarding your project. You are attempting to build a DLL (MEX file) for use with Matlab, using Matlab 2021, IMSL 7.0 (very old) and MKL2019. One problem with such a task is that Matlab and IMSL come with their own bundled versions of MKL, and calling routines from more than one MKL from the same EXE or DLL is quite troublesome. You have to make sure that you are consistent with the integer arguments (i.e., LP64 vs ILP64 models), and avoid causing the wrong versions of MKL to be used. On top of all this, it appears that your DLL could use an SQL database through ODBC!

Try the self-contained C example program that I gave above. Next, write a small driver in C or Fortran to call the spline fitting subprogram code that you posted, passing to it the same kind of x-y data that your IMSL calls used. At this stage, there should be no connection to IMSL, Matlab or SQL databases. When debugging your code, do not use compiler options that cause IPO to be performed.

After you get these small test projects to work, you may attempt to make similar changes and build the complex MEX file.

0 Kudos
geerits
New Contributor I
2,213 Views

Your comment

 

"You are attempting to build a DLL (MEX file) for use with Matlab, using Matlab 2021, IMSL 7.0 (very old) and MKL2019. One problem with such a task is that Matlab and IMSL come with their own bundled versions of MKL, and calling routines from more than one MKL from the same EXE or DLL is quite troublesome. You have to make sure that you are consistent with the integer arguments (i.e., LP64 vs ILP64 models), and avoid causing the wrong versions of MKL to be used."

 

is not what I like to hear BUT is extremely useful. Thanks!! I didn't know that. I was hoping I could use IMSL and MKL simultaneously as this would minimize the amount of work I would have to do in testing the MKL versus the IMSL version. On the other hand my code is working perfectly fine while using IMSL and MATLAB (even when linking to IMSL libraries AND using them. But not the data fitting routines but the SVD routine zgesvd()). Coincidence??

So, from your response I take that creating a new project while excluding IMSL still does not guarantee it to work because of the MATLAB libraries I have to link with?? Moreover, the problem you describe appears to be MKL routine dependent. You may recall my earlier problem, where I replaced IMSL routine DLSVCR(..) with MKL routine zgesvd() for the purpose of SVD (you helped me with that as well :-)). This worked just fine while linking to the IMSL, MKL and MATLAB libraries.

 

Your help is very much appreciated!

 

Best regards,

Tim.

 

 

0 Kudos
mecej4
Honored Contributor III
2,202 Views

What little I know about your project is not enough to provide a broad understanding of which pieces of IMSL, MKL and Matlab are being used and how those pieces interact.

It is possible to use versions of the IMSL libraries that do not use MKL; this may be a good choice if the IMSL routines are used only a few times in each run and consume a minor part of the execution time. Look at the LINK_FNL_SHARED_IMSL option.

If that is not the case, you may consider removing the IMSL dependency altogether. This choice may require reprogramming (as you have tried for the spline interpolation) to use MKL or some other replacement code in C or Fortran that provides equivalent capability to what IMSL is providing now.

If you have to use IMSL, MKL, Matlab and ODBC together, you may need to do considerable work to coordinate the interfaces and arrange so that "crossed wires" do not happen.

0 Kudos
geerits
New Contributor I
2,191 Views

Thanks again for your prompt reply but you bring me in a difficult position, especially because there can still be 'cross-talk' between MKL and MATLAB libraries. I'm not only talking about this project but also upcoming projects. We, scientists (i.e., Geophysicists/Physicists) typically prototype in MATLAB. The 'heavy' computational stuff is written in C/C++ and/or Fortran (the latter mostly legacy code) and wrapped as a mex file (i.e., MATLAB dll). We call mex files from Matlab and we write additional high level MATLAB code for more simple mathematical tasks and plotting. As to the heavy computational tasks we like to use professional libraries, developed by mathematicians (experts), so we do not have to worry about any details. From my 20+ years of experience I must say that IMSL has never resulted in 'unpleasant' surprises. Unfortunately, I cannot say that of MKL.

 

PS: My last and final step will be to create a new project, remove all IMSL dependencies and substitute the MKL. I will have to keep the MATLAB library dependencies though. If that doesn't work I think that the total time spent on this will justify us simply buying a company wide IMSL license:-).

 

Best regards,

Tim.

mecej4
Honored Contributor III
2,270 Views

Note to Gennady F. or other MKL engineers:

On Windows, using a 32-bit command window for OneAPI, the command

     icl /Qmkl splfit.c

works fine, but the similar command with icx fails. It appears that the correct MKL library names are not passed to the linker. If the OBJ file is generated separately using icl, however, icx is then able to link it with the correct libraries and produce an EXE.

 

 

T:\LANG\MKL>icx /Qmkl splfit.c
Intel(R) oneAPI DPC++/C++ Compiler for applications running on IA-32, Version 2022.0.0 Build 20211123
Copyright (C) 1985-2021 Intel Corporation. All rights reserved.

LINK : fatal error LNK1104: cannot open file 'mkl_intel.lib'
clang-cl: error: linker command failed with exit code 1104 (use -v to see invocation)

T:\LANG\MKL>icl /Qmkl /c splfit.c
Intel(R) oneAPI DPC++/C++ Compiler for applications running on IA-32, Version 2022.0.0 Build 20211123
Copyright (C) 1985-2021 Intel Corporation.  All rights reserved.

T:\LANG\MKL>icx /Qmkl splfit.obj
Intel(R) oneAPI DPC++/C++ Compiler for applications running on IA-32, Version 2022.0.0 Build 20211123
Copyright (C) 1985-2021 Intel Corporation. All rights reserved.

 

 

This problem is not seen in a 64-bit command window configured for OneAPI and MKL.

 

0 Kudos
Gennady_F_Intel
Moderator
2,163 Views

thanks, Mecej4. That's true. the icx driver for 32 bit doesn't recognize the /Qmkl option.  The issue is escalated versus the compiler team.

In the case, if we could use the "standard" linking line, as follows: icx C:\..\mkl\latest\include testmkl.c MKLROOT32\mkl_core.lib MKLROOT32\mkl_intel_c.lib MKLROOT32\mkl_sequential.lib  then it works as expected.

Gennady_F_Intel
Moderator
2,143 Views

This issue has been resolved and we will no longer respond to this thread. If you require additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only. 



Reply