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

Data fit dfdInterpolateEx1D function warning code

jelich__christopher
753 Views

Hello everybody, I am having some problems with the dfdInterpolateEx1D function of the MKL Data Fit module. I am calling the above function, togheter with all the relative initialization/setting/clean function, within a loop and the results are sometimes corrects and sometimes wrong. Particularly, the wrong results provided by the above function are always equal to 0 and in correspondence to these events the error/warning status is equal to 10. As described in the library manual, positive status are associated to warning messages however I can not find the warning status 10 in the mkl_df_defines.h header file which, according to the manual itself, contains all the defined error/warning status. What can I do to understand what is wrong with my code?

The version of the library is defined as INTEL_MKL_VERSION = 110101, I am using the 32 bit MKL interface and I am compiling/linking my code with g++ 4.6.3 on an Ubuntu machine.

Many thanks in advance for your help

Christopher

0 Kudos
10 Replies
Andrey_N_Intel
Employee
753 Views

Hi Christopher,

Can you please send the example which demonstrates the issue? Also, specify your CPU and type of linking against Intel MKL, static or dynamic.

Thanks,

Andrey

0 Kudos
jelich__christopher
753 Views

Hi Andrey, thank you for your reply. I am dinamically linking against MKL and my CPU is an Intel Core2 Quad CPU Q9550 @2.83GHz. Concerning the example: I am using the Data Fit library to perform a linear interpolation of a 1D complex function by making 2 interpolation steps, one for the real part and one for the imaginary part. It would be better, I guess, making a single interpolation over a 2D real function whose components are defined as the real and imaginary part of my original complex function. I will try this after solving my problems on the 1D case. The code is the following:

void interp(double *X, complexType<double> *Y, double *XI, complexType<double> *YI, int Nx, int Nxi, const complexType<double> *extrapval)
{

int status; 
DFTaskPtr task;

MKL_INT nx = Nx;
MKL_INT xhint = DF_NO_HINT;

MKL_INT ny = 1;
MKL_INT yhint = DF_NO_HINT;   

MKL_INT s_order = DF_PP_LINEAR;
MKL_INT s_type = DF_PP_DEFAULT;
MKL_INT  bc_type = DF_BC_NOT_A_KNOT;
double* bc = NULL;
MKL_INT  ic_type = DF_NO_IC;
double* ic = NULL;
double* scoeff = new double[(nx-1)* SPLINE_ORDER];
MKL_INT scoeffhint = DF_NO_HINT;

MKL_INT compType = DF_INTERP;
MKL_INT siteHint = DF_NO_HINT;
MKL_INT ndorder = 1;
MKL_INT dorder[] = {0};
double* datahint = NULL;
MKL_INT rhint = DF_NO_HINT;
MKL_INT* cell = NULL;

dfdInterpCallBack extraFunc;
extraFunc = reinterpret_cast<dfdInterpCallBack>(linear_extrap_func);

double* complexComponentY = new double[nx];
for(int complexIndex=0;complexIndex<nx;++complexIndex)
{
  complexComponentY[complexIndex] = Y[complexIndex].re;
}

status = dfdNewTask1D(&task,
		      nx,
		      X,
		      xhint,
		      ny,
		      complexComponentY,
		      yhint);

status = dfdEditPPSpline1D(task,
			   s_order,
			   s_type,
			   bc_type,
			   bc,
			   ic_type,
			   ic,
			   scoeff,
			   scoeffhint);

status = dfdConstruct1D(task,
                        DF_PP_SPLINE,
		        DF_METHOD_STD);


double* complexComponentYI = new double[Nxi];
status = dfdInterpolateEx1D(task,
		            compType,
		    	    DF_METHOD_PP,
		    	    Nxi,
		    	    XI,
		    	    siteHint,
		    	    ndorder,
		    	    dorder,
			    datahint,
			    complexComponentYI,
			    rhint,
			    cell,
			    extraFunc,
			    (void*)extrapval,
			    extraFunc,
			    (void*)extrapval,
			    0,
			    0,
			    0,
			    0);

for(int complexIndex=0;complexIndex<Nxi;++complexIndex)
{
  YI[complexIndex].re = complexComponentYI[complexIndex];
}

for(int complexIndex=0;complexIndex<nx;++complexIndex)
{
  complexComponentY[complexIndex] = Y[complexIndex].im;
}

status = dfDeleteTask(&task);

status = dfdNewTask1D(&task,
                      nx,
		      X,
		      xhint,
		      ny,
		      complexComponentY,
		      yhint);

status = dfdEditPPSpline1D(task,
                           s_order,
		           s_type,
			   bc_type,
			   bc,
			   ic_type,
			   ic,
			   scoeff,
			   scoeffhint);

status = dfdConstruct1D(task,
                        DF_PP_SPLINE,
			DF_METHOD_STD);

status = dfdInterpolateEx1D(task,
                            compType,
		    	    DF_METHOD_PP,
		    	    Nxi,
		    	    XI,
		    	    siteHint,
		            ndorder,
		    	    dorder,
			    datahint,
			    complexComponentYI,
			    rhint,
			    cell,
			    extraFunc,
			    (void*)extrapval,
			    extraFunc,
			    (void*)extrapval,
			    0,
			    0,
			    0,
			    0);

for(int complexIndex=0;complexIndex<Nxi;++complexIndex)
{
  YI[complexIndex].im = complexComponentYI[complexIndex];
}

delete [] complexComponentY;
delete [] complexComponentYI;
delete [] scoeff;

status = dfDeleteTask(&task);

return;

}    

The complex function is defined by the variable pairs (X,Y), with X ordered increasingly, the interpolation points are strored in the XI array and the interpolated value in the YI one.

UPDATE: I made a mistake somewhere, the returned status is the one returned by my auxiliary function used for extrapolation. However, there is still something not working properly in my code since I am not trying to interpolate outside the definition range of my complex function. Concerning the above auxiliary function, I am using the one defined here: https://software.intel.com/fr-fr/forums/topic/517073.

Thank you for your help.

Christopher

 

 

0 Kudos
Andrey_N_Intel
Employee
753 Views

Hi Christopher, thanks for the code and the additional details. We will have a look. Andrey

0 Kudos
VictoriyaS_F_Intel
753 Views

Hello Christopher,

I have tried to reproduce the issue using the code you have provided. But that code contains only the interpolation function and does not contain data structures that are used in the code, values of the input parameters of the function interp(), etc.

Could you please provide the reproducer that could be compiled and run on our side, i.e. the one that will contain the main() function?

Best regards,

Vika

0 Kudos
jelich__christopher
753 Views

Hello Victoriya,

many thanks for your reply and for your help. I apologize for not having reported a complete example, you can find it here below:

#include <fstream>
#include <complex>
#include <iostream>
#include "mkl_df.h"
#define SPLINE_ORDER DF_PP_LINEAR

int linear_extrap_func( MKL_INT64* n, MKL_INT64 cell[], double site[], double r[], void *params )
{
	r[0] = 0.;
	return 1;
}

void interp1(double *X, std::complex<double> *Y, double *XI, std::complex<double> *YI, int Nx, int Nxi, const std::complex<double> *extrapval)
{

	int status;
	DFTaskPtr task;

	MKL_INT nx = Nx;
	MKL_INT xhint = DF_NO_HINT;

	MKL_INT ny = 1;
	MKL_INT yhint = DF_NO_HINT;

	MKL_INT s_order = DF_PP_LINEAR;
	MKL_INT s_type = DF_PP_DEFAULT;
	MKL_INT bc_type = DF_BC_NOT_A_KNOT;
	double* bc = NULL;
	MKL_INT ic_type = DF_NO_IC;
	double* ic = NULL;
	double* scoeff = new double[(nx-1)* SPLINE_ORDER];
	MKL_INT scoeffhint = DF_NO_HINT;

	MKL_INT compType = DF_INTERP;
	MKL_INT siteHint = DF_NO_HINT;
	MKL_INT ndorder = 1;
	MKL_INT dorder[] = {0};
	double* datahint = NULL;
	MKL_INT rhint = DF_NO_HINT;
	MKL_INT* cell = NULL;

 	dfdInterpCallBack extraFunc = linear_extrap_func;

	double* complexComponentY = new double[nx];
	for(int complexIndex=0;complexIndex<nx;++complexIndex)
	{
	  complexComponentY[complexIndex] = Y[complexIndex].real();
	}

	status = dfdNewTask1D(&task,
			      nx,
			      X,
			      xhint,
			      ny,
		              complexComponentY,
			      yhint);

	status = dfdEditPPSpline1D(task,
				   s_order,
			           s_type,
				   bc_type,
				   bc,
			           ic_type,
				   ic,
				   scoeff,
				   scoeffhint);

	status = dfdConstruct1D(task,
				DF_PP_SPLINE,
				DF_METHOD_STD);


	double* complexComponentYI = new double[Nxi];
	status = dfdInterpolateEx1D(task,
				    compType,
				    DF_METHOD_PP,
				    Nxi,
				    XI,
				    siteHint,
				    ndorder,
				    dorder,
				    datahint,
				    complexComponentYI,
				    rhint,
				    cell,
				    extraFunc,
				    (void*)extrapval,
				    extraFunc,
				    (void*)extrapval,
				    0,
				    0,
				    0,
				    0);

	for(int complexIndex=0;complexIndex<Nxi;++complexIndex)
	{
	  YI[complexIndex].real() = complexComponentYI[complexIndex];
	}

	for(int complexIndex=0;complexIndex<nx;++complexIndex)
	{
	  complexComponentY[complexIndex] = Y[complexIndex].imag();
	}

	status = dfDeleteTask(&task);

	status = dfdNewTask1D(&task,
			      nx,
			      X,
		              xhint,
			      ny,
			      complexComponentY,
			      yhint);

	status = dfdEditPPSpline1D(task,
			           s_order,
				   s_type,
				   bc_type,
				   bc,
				   ic_type,
				   ic,
				   scoeff,
				   scoeffhint);

	status = dfdConstruct1D(task,
				DF_PP_SPLINE,
				DF_METHOD_STD);

	status = dfdInterpolateEx1D(task,
				    compType,
				    DF_METHOD_PP,
			            Nxi,
				    XI,
				    siteHint,
				    ndorder,
				    dorder,
				    datahint,
				    complexComponentYI,
				    rhint,
				    cell,
				    extraFunc,
				    (void*)extrapval,
				    extraFunc,
				    (void*)extrapval,
				    0,
				    0,
				    0,
				    0);

	for(int complexIndex=0;complexIndex<Nxi;++complexIndex)
	{
	  YI[complexIndex].imag() = complexComponentYI[complexIndex];
	}

	delete [] complexComponentY;
	delete [] complexComponentYI;
	delete [] scoeff;

	status = dfDeleteTask(&task);

	return;

}


int main()
{
	const unsigned int Nxi(1);
	const unsigned int Nx(582);

	double* X = new double[Nx];
	std::complex<double>* Y = new std::complex<double>[Nx];

	double* XI = new double[Nxi];
	XI[0] = 64732.867582857019;
	std::complex<double>* YI = new std::complex<double>[Nxi];

	// Correct interpolation result:
	//-1.5086179120445398e-06 + 3.9879957076769775e-07*i

	std::ifstream data;
	data.open("data.dat");
	for(unsigned int i=0;i<Nx;++i){
		data>>X>>Y.real()>>Y.imag();
	}
	data.close();

	std::complex<double> extrapValue(0.,0.);
	interp1(X, Y, XI, YI, Nx, Nxi, &extrapValue);

	std::cout<<YI[0]<<std::endl;

	delete [] X;
	delete [] Y;
	delete [] XI;
	delete [] YI;

	return 0;
}

The (ASCII) data.dat file is attached to the present message. Hardcoded in the main function you can find an example of value for the interpolation point (XI) and the expected result. As you can see the XI value is within the interpolated function (X,Y) definition range, however the extrapolation function is called. The same behavior is observed also changing the value of XI.

Many thanks again for your help!

Best regards

Christopher

0 Kudos
jelich__christopher
753 Views

Here is the data file containing the interpolation function definition, I did not uploaded it with my previous message.

Best regards

Christopher

0 Kudos
VictoriyaS_F_Intel
753 Views

Hello Christopher,

I have run the example you have provided with MKL 11.1 Update 1 and got (0,0) as a result. Indeed, there was a bug in MKL 11.1 Update 1 in df?InterpolateEx1D function.

But it was fixed in MKL 11.2 Update 1. When I have run the example with MKL 11.2 Update 1, I have got the correct results: (-1.50862e-06,3.988e-07).

Could you please check that the issue is fixed on your side with MKL 11.2 Update 1?

Best regards,

Vika

0 Kudos
jelich__christopher
753 Views

Hello Victoriya,

thank you for your help, I can confirm that with MKL 11.2 (from Composer XE 2015) I have got the correct results too. Now I will try to make some more extensive tests with 2D functions but in the meanwhile everything seems to work fine.

Many thanks again.

Best regards

Christopher

0 Kudos
jelich__christopher
753 Views

Hello everybody again,

I have a new question concerning the Data Fit module, in this case for the interpolation of N-dimensional functions (2-d in my case). I have tried to modify the already posted code to perform the interpolation of a 2-d function defined as the (real,imag) component pairs of the previously used complex function. The modified code is the following: 

#include <fstream>
#include <iostream>
#include "mkl_df.h"
#define SPLINE_ORDER DF_PP_LINEAR

const unsigned int FUNCTION_SIZE(1);

int linear_extrap_func( MKL_INT64* n, MKL_INT64 cell[], double site[], double r[], void *params )
{
	r[0] = 0.;

	if(FUNCTION_SIZE == 2)
	{
		r[1] = 0.;
	}

	return 1;
}

void interp1(double *X, double *Y, double *XI, double *YI, int Nx, int Ny, int Nxi, double *extrapval)
{

	int status;
	DFTaskPtr task;

	MKL_INT nx = Nx;
	MKL_INT xhint = DF_NO_HINT;

	MKL_INT ny = Ny;
	MKL_INT yhint = DF_NO_HINT;

	MKL_INT s_order = SPLINE_ORDER;
	MKL_INT s_type = DF_PP_DEFAULT;
	MKL_INT bc_type = DF_BC_NOT_A_KNOT;
	double* bc = NULL;
	MKL_INT ic_type = DF_NO_IC;
	double* ic = NULL;
	double* scoeff = new double[(nx-1)*SPLINE_ORDER];
	MKL_INT scoeffhint = DF_NO_HINT;

	MKL_INT compType = DF_INTERP;
	MKL_INT siteHint = DF_NO_HINT;
	MKL_INT ndorder = 1;
	MKL_INT dorder[] = {0};
	double* datahint = NULL;
	MKL_INT rhint = DF_NO_HINT;
	MKL_INT* cell = NULL;

 	dfdInterpCallBack extraFunc = linear_extrap_func;

	status = dfdNewTask1D(&task,
			      nx,
			      X,
			      xhint,
		              ny,
			      Y,
			      yhint);

	std::cout<<"dfdNewTask1D status: "<<status<<std::endl;

	status = dfdEditPPSpline1D(task,
				   s_order,
				   s_type,
				   bc_type,
				   bc,
				   ic_type,
			           ic,
				   scoeff,
				   scoeffhint);

 	std::cout<<"dfdEditPPSpline1D status: "<<status<<std::endl;

	status = dfdConstruct1D(task,
				DF_PP_SPLINE,
				DF_METHOD_STD);

        std::cout<<"dfdConstruct1D status: "<<status<<std::endl;

 	status = dfdInterpolateEx1D(task,
				    compType,
				    DF_METHOD_PP,
				    Nxi,
				    XI,
				    siteHint,
				    ndorder,
				    dorder,
				    datahint,
				    YI,
				    rhint,
				    cell,
				    extraFunc,
				    (void*)extrapval,
				    extraFunc,
				    (void*)extrapval,
				    0,
				    0,
				    0,
				    0);

 	std::cout<<"dfdInterpolateEx1D status: "<<status<<std::endl;

	status = dfDeleteTask(&task);

 	std::cout<<"dfDeleteTask status: "<<status<<std::endl;

	return;
}


int main()
{
	if(FUNCTION_SIZE != 2 && FUNCTION_SIZE != 1)
	{
		std::cout<<"Error in function dimension"<<std::endl;
		return 1;
	}

	const unsigned int Nxi(1);
	const unsigned int Nx(582);
	const unsigned int Ny(FUNCTION_SIZE);

	double* X = new double[Nx];
	double* Y = new double[Nx*Ny];

	double* XI = new double[Nxi];
	XI[0] = 64732.867582857019;
	double* YI = new double[Nxi*Ny];

	// Correct interpolation result:
	//(-1.5086179120445398e-06, 3.9879957076769775e-07)

	std::ifstream data;
	data.open("data.dat");
	if(FUNCTION_SIZE == 2)
	{
		unsigned int j(0);
		for(unsigned int i=0;i<Nx;++i,j+=Ny){
			data>>X>>Y>>Y[j+1];
		}
	}
	else
	{
		double dummy;
		for(unsigned int i=0;i<Nx;++i){
			data>>X>>Y>>dummy;
		}

	}
	data.close();

	double extrapValue(0.);
	interp1(X, Y, XI, YI, Nx, Ny, Nxi, &extrapValue);

	if(FUNCTION_SIZE == 2)
	{
		std::cout<<YI[0]<<" "<<YI[1]<<std::endl;
	}
	else
	{
		std::cout<<YI[0]<<std::endl;
	}

	delete [] X;
	delete [] Y;
	delete [] XI;
	delete [] YI;

	return 0;
}

The same data of the previous test is used and the function dimension (1 or 2 only) can be modified making use of the FUNCTION_SIZE global variable.

Running the code with FUNCTION_SIZE=1 the correct result is obtained. Instead, running the interpolation with FUNCTION_SIZE=2 I observe the following behavior: without any modification to the code above there is a segmentation fault within the dfdConstruct1D function. With a larger memory allocation for the scoeff array, particularly changing

double* scoeff = new double[(nx-1)*SPLINE_ORDER];

into

double* scoeff = new double[ny*(nx-1)*SPLINE_ORDER];

the segmentation fault disappears but a wrong result (0.) is obtained for the second component of the interpolated function. Most likely there is something wrong with my parameters however I can not find the error. Particularly, if I am not wrong, my 2-d function definition is in the row-major format and therefore in the default format for the Data Fit function used by my code. Moreover, according to the library documentation, the correct size for the scoeff array should be s_order*(nx-1).

Can you see something wrong in my code?

Many thanks in advance for your attention!

Best regards

Christopher

P.S.: the topic of the discussion is quite different than the initial one, please let me know if I should open a new discussion!

 

0 Kudos
VictoriyaS_F_Intel
753 Views

Christopher,

You are right that the size of scoeff array should be ny*(nx-1)*SPLINE_ORDER. There is probably a bug in Data Fitting documentation. Thank you for noticing this.

Can you see something wrong in my code?

Yes, there is a problem with Y array that is passed to dfdNewTask1D function. As the default storage format for Y array is row-major, following loop should be changed.

		unsigned int j(0);
		for(unsigned int i=0;i<Nx;++i,j+=Ny) {
    			data>>X>>Y>>Y[j+1];
		}

The correct version is:

		for(unsigned int i=0;i<Nx;++i){
			    data>>X>>Y>>Y[Nx+i];
		}

Another problem is that there is a known bug in Data Fitting dfdInterpolateEx1D function that causes incorrect results when ny > 1. This bug will be fixed in the one of the upcoming MKL releases. Until then you could use either 
dfdInterpolate1D function where it is possible, or call dfdInterpolateEx1D with ny = 1 multiple times to obtain interpolation results for each dimension separately.

Best regards,
Vika

0 Kudos
Reply