#include <iostream>
#include <vector>
#include <mkl.h>
#include <mkl_df.h>

using namespace std;

int main () {
    const int nx = 11; // #rows in x
    const int ny = 1; // #cols in y
    
    vector<double> x(nx), y(nx);
    for(int i=0; i<nx; ++i){
        x[i] = i;
        y[i] = i;
    }

    DFTaskPtr task;
    auto xhint = DF_NON_UNIFORM_PARTITION;
    auto yhint = DF_MATRIX_STORAGE_COLS;
    int status = dfdNewTask1D(&task, nx, &x[0], xhint, ny, &y[0], yhint);


    auto s_order = DF_PP_LINEAR;
    auto s_type  = DF_PP_DEFAULT;
    auto bc_type = DF_BC_NOT_A_KNOT;
    auto ic_type = DF_NO_IC;
    auto scoeffhint = DF_NO_HINT;
    vector<double> coeffs(ny*(nx-1)*s_order);
    status = dfdEditPPSpline1D(task, s_order, s_type, bc_type, nullptr, ic_type, nullptr, &coeffs[0], scoeffhint);
    status = dfdConstruct1D(task, DF_PP_SPLINE, DF_METHOD_STD);
    
    const MKL_INT nsite = 1;
    vector<double> xi(nsite), yi(nsite);
    MKL_INT ndorder = 1;
    MKL_INT dorder = 1;
    double* datahint;
    datahint = DF_NO_APRIORI_INFO;
    MKL_INT rhint = DF_MATRIX_STORAGE_ROWS;
    MKL_INT sitehint = DF_NON_UNIFORM_PARTITION;
    xi[0] = 10;

    status = dfdInterpolate1D(task,
                              DF_INTERP,
                              DF_METHOD_PP,
                              nsite,
                              &xi[0],
                              sitehint,
                              ndorder,
                              &dorder,
                              datahint,
                              &yi[0],
                              rhint,
                              NULL);


    status = dfDeleteTask(&task);

    for(int i=0; i<nsite; ++i){
        cout << yi[i] << ", ";
    }

    return 0;
}