Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Arnab_D_
Beginner
65 Views

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

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;
}

0 Kudos
1 Reply
65 Views

Hello Arnab,

We have started to analyze your test case, and here are the preliminary results:

For the Hermite spline we require setting the internal conditions (ic) which are the function's 1-st derivative values at the internal points of the interpolation interval. The code you have provided does not set those internal conditions. This can be the reason why you get incorrect results.

Best regards,

Victoriya

Reply