MKL library (Data Fitting Component)- substitute of P-chip function in MKL library

Arnab_D_ — Fri, 17 Mar 2017 11:37:16 GMT

Dear Experts,

I want to implement the p-chip MATLAB function using the Intel MKL library. The details of the pchip function can be found at :https://in.mathworks.com/help/matlab/ref/pchip.html

I tried using the Hermite Cubic Spline but it returns zero values in the interpolated function. Why does it return zero values.

Can someone point to such an implementation, if you have tried previously ?

The Code

// This is the main DLL file.

#include <algorithm>
#include <iostream>

#include "mkl_df.h"

#include "MKLWrapper64.h"
#include "Utils.h"

using namespace MKLWrapper64;

/// Performs interpolation of up to order 3,
/// with free-end boundary condition (i.e. second derivative of both ends = 0)
/// INPUT PARAMETER RESTRICTIONS:
/// --> x[] must be in ascending order
/// --> r[] must be of the format {lowerBound,upperBound} if isUniform
/// or contain outsize number of x-values if !isUniform
/// --> order should be within the range [1, 3]
/// OUTPUT PARAMETERS:
/// --> out[] contains output interpolated y-values
/// --> return int signifies no error (0) or error (non-zero status number)

int Interpolator::Interpolate(int nx, array<float>^ x, array<float>^ y,

   array<float>^ r, int outsize, array<float>^ out, int order, bool isUniform)
{
   int status;
   int ny = 1; // y-axis is 1 dimensional
   DFTaskPtr taskPtr;
   pin_ptr<float> xPtr = &x[0];
   pin_ptr<float> yPtr = &y[0];
   pin_ptr<float> interpolationRange = &r[0];

   // create new data fitting task
   status = dfsNewTask1D(&taskPtr, nx, xPtr, DF_UNIFORM_PARTITION, ny, yPtr, DF_NO_HINT);
   if (Utils::checkError(status, "dfsNewTask1D") != 0) return status;

   // specify type of spline
   int splineOrder;
   switch (order) {
   case 1: splineOrder = DF_PP_LINEAR;   break;
   case 2: splineOrder = DF_PP_QUADRATIC; break;
   case 3: splineOrder = DF_PP_CUBIC;       break;
   default: splineOrder = DF_PP_CUBIC;       break;
   }

   int sorder = DF_PP_CUBIC;
   int stype = DF_PP_HERMITE;
   int bc_type = DF_BC_1ST_LEFT_DER | DF_BC_1ST_RIGHT_DER;
   int ic_type = DF_IC_1ST_DER;

   float* splineCoeff = new float[ny*splineOrder*(nx - 1)];
   status = dfsEditPPSpline1D(taskPtr, sorder, stype,
       bc_type, NULL, ic_type, NULL, splineCoeff, DF_NO_HINT);
   if (Utils::checkError(status, "dfsEditPPSpline1D") != 0) return status;

   // construct spline
   status = dfsConstruct1D(taskPtr, DF_PP_SPLINE, DF_METHOD_STD);
   if (Utils::checkError(status, "dfsConstruct1D") != 0) return status;

   // interpolate
   int intUniform = isUniform ? DF_UNIFORM_PARTITION : DF_NON_UNIFORM_PARTITION;
   MKL_INT evalOrder = DF_PP_STD + 1; // evaluate 0th derivative of spline (i.e. spline value itself)
   MKL_INT* dataOrder = new MKL_INT[outsize]; // instructs dfdInterpolate1D to evaluate 0th derivative at all points
   std::fill_n(dataOrder, outsize, evalOrder);
   float* resultArr = new float[outsize];
   status = dfsInterpolate1D(taskPtr, DF_INTERP, DF_METHOD_PP, outsize, interpolationRange,
       intUniform, evalOrder, dataOrder, nullptr, resultArr, DF_MATRIX_STORAGE_ROWS, NULL);
   if (Utils::checkError(status, "dfsInterpolate1D") != 0) return status;

   // delete data fitting task
   status = dfDeleteTask(&taskPtr);
   if (Utils::checkError(status, "dfDeleteTask") != 0) return status;

   for (int i = 0; i < outsize; i++)
       out = resultArr;
   return 0;
}

Hello Arnab,

VictoriyaS_F_Intel — Mon, 20 Mar 2017 16:02:46 GMT

Hello Arnab,

Per our analysis of your test case, it requires several modifications:

As it was mentioned in the previous response, to construct Hermite cubic spline coefficients one is required to provide the array of internal conditions (ic). You can find the implementation details in dfshermitecubicspline.c example available in Intel MKL.
When the boundary condition is set to DF_BC_1ST_LEFT_DER | DF_BC_1ST_RIGHT_DER you should provide an array of two elements which specifies the values of the function’s first derivatives in the endpoints of the interpolation interval.
If you create Data Fitting task by providing uniform partition of the interpolation interval [a,b]:

status = dfsNewTask1D(&taskPtr, nx, xPtr, DF_UNIFORM_PARTITION, ny, yPtr, DF_NO_HINT);

you are required to provide the array x that defines the partition as a two-element array that represents the endpoints of the interpolation interval [a, b].

Please find the additional details about setting internal and boundary conditions for the different types of splines at https://software.intel.com/en-us/node/522222#D69BA4E7-E8BD-4413-A2A9-769D5002F832

Extra details about using different types of partitions for splines construction and spline-based interpolation are available at https://software.intel.com/en-us/node/522220

Please, let us know, if it helps address your questions on your test case. Also, let us know, if you have more comments or questions on the Data Fitting component of Intel MKL.

topic MKL library (Data Fitting Component)- substitute of P-chip function in MKL library in Intel® oneAPI Math Kernel Library

MKL library (Data Fitting Component)- substitute of P-chip function in MKL library

Hello Arnab,