Access violation in spline function

[netphilou31]

Hi,

In an effort to remove IMSL library from our projects (because of the new licensing policy) I am giving a try to the MKL library. I mainly use cubic spline functions using Intel Visual Fortran Compiler XE on IA-32, version 12.1.1 (Package ID: w_fcompxe_2011.7.258). I have started with from an existing project that compiles and links fine with IMSL (the project exports functions used in Excel through VBA to provide spline features). I have added a new function that uses the spline possibilities of the MKL (10.3 update 10). The call to the function from VBA works fine, but I get an Access violation when calling the function dfdConstruct1D(). It seems that the problem is located in the mkl at location : _mkl_df_kernel_n8_dDFAkimaCubicSpline1D() + 0x1765c octets. The source code was created from an example provided in the documentation and I don't understand what I am doing wrong. The source file is attached.

Best regards.

[mecej4]

It would help matters if you can post a small Fortran driver routine that calls your DLL with argument values that cause the access violation. Presumably, not every call to the DLL causes access violatons, so we need to have a sample case to reproduce the problem.

There is a Matlab implementation of one-dimensional Akima interpolation which you may find helpful.

[netphilou31]

Hi mece,

Thanks for your reply. Please findattached a ZIP file containingthe IVF project as well as the Excel file that calls the dll with the data used in sheet1 (you will have to change the location of the library in the VBA code). I did not include the module file (MKL_DF) whichis directly accessed from its original location.

Best regards,

[netphilou31]

Hi,

Just for giving some additional information. After the access violation message, if I stop into the debugger and loork at the contents of the scoeff array (coefficients of the piecewise cubic spline), it shows correct information: i.e. if I copy these coefficients, paste them into the excel worksheet and try to calculate the spline function, it gives correct results.

Best regards,

[Sergey_M_Intel2]

Hi,

Many thanks for additional detail. Please give us some time to reproduce & root cause the issue.

Sergey

[VictoriyaS_F_Intel]

Hi netphilou31,

There're several issues with the subroutine AkimaSpline_MKL that you're using for Akima cubic interpolation.

1. Before calling dfdNewTask1D() function you set ny=nx. This is wrong, because setting ny=nx means that you will construct coefficients for nx Akima cubic splines on the same partition at once. But you only need coefficients for 1 spline. So, to fix the issue please make following change in your code:

! Create Data Fitting task

ny = 1

errcode = dfdNewTask1D(task, nx, xx, xhint, ny, yy, yhint)

This was the cause of access violation in you application. There were no enough memory for construction of nx splines.

2. datahint parameter is initialized incorrectly when calling dfdInterpolate1D. You should either skip it, because it's unnecessary in your case, or make following change to the code:

datahint(1) = 1.d0 ! Dimension of the task

datahint(2) = dble(DF_NO_APRIORI_INFO) ! No additional information about breakpoints or sites is provided.

3. Interpolation results are packed into r array - parameter of dfdInterpolate1D function, they aren't returned as the functions' result.

It will be better to use following code sequence to obtain interpolation result:

errcode = dfdInterpolate1D(task, DF_INTERP, DF_METHOD_PP, nsite, site, sitehint, ndorder, dorder, datahint, r, rhint)

call CheckDfError(errcode)

Y = r(1)

After applying those three changes described above the code should start working as expected.

Please seefixed file in attachment.

[netphilou31]

Hi Victoriya,

Thanks a lot ! I did very basic mistakes that I should not have done !
I have changed the code according to your suggestions but unfortunately the dfdInterpolate1D function returns always 0 in r(1) whatever the value stored in site(1). If I have correctly understood, the site() array contains the x values for which one wants an interpolation and the r() array contains the results ?

Best regards,

[netphilou31]

Ok,

I have found that I have to set the dorder(1) value to 1 to get the interpolated valued in r() instead of 0 as mentionned in the documentation :

dorder: Array of size ndorder that defines the order of the derivatives to be computed at the interpolation sites. If all the elements in dorder are zero, the library computes the spline values only. If you do not need interpolation computations, set ndorder to zero and pass a NULL pointer to dorder.

Moreover, but it is probably not a problem with the example, I get very different results from the DCSAKM function from the IMSL library ! The DCSAKM function gives a very good spline but the one I have used gives a bad one. Probably IMSL is using constraints as first or second derivatives at breakpoints to get better results. Il will check that and try, if possible, to fill the gap between both spline results.

Best regards,

[Sergey_M_Intel2]

Hello again,

Please let us know if you still see significant deviation between libraries. If you get stuck then we might have a deeper look at the issue.

Sergey

[netphilou31]

Hi Sergey,

Please find attached some comparisons I have made between the MKL and IMSL libraries Akima spline functions. As mentioned in my previous message, it is possible that IMSL uses some kind of constrained spline. I have also attached a PDF describing another spline method that you can have a look on. When programmed this method gives very ggod results (even better than the IMSL Akima spline). If you need the Excel workbook and the IVF project I have used please tell me.

Best regards,

Phil.

P.S: For your information, the natural splines give almost the same results in both libraries (and MKL natural spline gives better results than Akima in MKL)

[VictoriyaS_F_Intel]

Hi Phil,

Thank you for your time spent to analyse the issue!
The data you've provided is enough to continue investigation of the differences in Akima splines computations on our side.

Thanks,
Vika

[netphilou31]

Hi Victoriya,

Please find attached another couple of comparisons I have made using the example provided in the IMSL documentation for the CSAKM/DCSAKM function. It seems that the problem lies at the end-points (the Akima spline is similar to the one from IMSL except at the end-points, i.e. in the first two and last two intervals).
Thanks again.

Best regards,

[mecej4]

Akima interpolation is a rather simple calculation, and it is a bit of overkill to use MKL or IMSL when a few lines of code (the file akima.m referred to in #1 has 21 lines of code) will calculate the Akima interpolant. Here are the lines of code (in Matlab) to fit and plot the result for your data using akima.m:

[bash]cd s:/akima x=0:0.1:1; y=sin(15*x); % 11 data points xi=linspace(0,1,101); yc=sin(15*x); % calculated yi=akima(x,y,xi); % interpolated plot(xi,yi,'-',x,y,'o',xi,yc) xlabel('x'); ylabel('y'); grid legend('Akima','Data','calculated') orient landscape print -dpng 's:/akima/akm.png' [/bash]

[Sergey_M_Intel2]

It is not always about simplicity of the code. General matrix multiplication is also less than 20 lines of the code but if you want to get closer to the peak performance then you'll likely use GEMM from MKL, right?

Likewise, if the spline construction and/or evaluation represents performance sensitive part of the application it's certainly make sense touse one fromMKL.

Regards,
Sergey

[netphilou31]

It is not really a question of performance of the algorithm but rather a question of fitting quality. When I used the IMSL library, I have found that the Akima spline function was very useful because it preserves the shape of the curve and avoids overshooting problems (as the ones described by CJC Kruger in the article attached in a previous reply). The same good fitting was done using the CSCON/DCSCON routine but slightly different from Akima because it preserves the concavity of the curve rather than the shape.
For my job it is important to avoid overshooting problems so I can use the algorithm proposed by CJC Kruger ; but, because I also want to remove other calls to IMSL functions, I was trying the MKL library. However, I have found that its use was very different from IMSL thus I ran into difficulties for implementing it.
Anyway, are you sure that there is no problems in the MKL for the fitting of the first and last two intervals (as shown on the graphs I have attached previously), because for the internal intervals the results are very close between both libraries. Moreover if the polynomial of the first interval is erroneous, the second is probably also erroneous (because of the continuity of the derivative at the breakpoint).

Best regards,

Phil.

[VictoriyaS_F_Intel]

Phil, hello again,

The examples provided by you disclosed the issue in MKL Akima spline construction, which wasn't caught by our internal testing. Thank you very much for your efforts!
Indeed, coeficients of Akima spline are calculated incorrectly for two beginning and two ending intervals. This issue isplanned to be fixed in further MKL releases. Hopefully you will try one ofthem.

Thanks,
Vika

[netphilou31]

Hi Victoriya,

Could it be possible to post a reply on this thread to be informed that the problem is fixed ?
Thanks.

Best regards,

Phil.

[mecej4]

Agreed. Is performance, indeed, an issue for the OP?

In any case, the Matlab code can be useful in providing test results to validate the results from MKL.

[Gennady_F_Intel]

Hello Phil, Thiis the real issue for Data Fitting component of MKL.The issue has been submitted to our internal development tracking database for further fixing. I will inform you once a new update becomes available. Here is a bug tracking number for your reference: DPD200288598

regards, Gennady

[netphilou31]

Hi,

Many thanks to all of you.
I will wait impatiently for the next release.

Best regards,

[netphilou31]

Hi,

To complete and to answer to mecej4, the Matlab code is very useful for testing this function because it shows that the routine dfdInterpolate1D returns NaN for a portion of the flat part of the curve when testing the Akima examples (see the attached plot in which I have replaced the NaN values by #VALUE into Excel).

Best regards,