#define N   10      // number of break points in partition
#define M   4       // number of interpolation sites in block
#define N_SITE_BL	10  // number of blocks of interpolation sites

#include <stdlib.h>
#include <stdio.h>

#include "omp.h"
#include "mkl_df.h"


void ThrInterp(int i_block, MKL_INT n, double *x, double *y, MKL_INT order,
               double *coeff, MKL_INT m, double *site, double *r);

int main()
{
    int i, errcode = 0;
    MKL_INT n = N;      // number of break points in partition
    MKL_INT m = M;      // number of interpolation sites in block
    MKL_INT order = DF_PP_LINEAR;   // spline order

    double *x, *y, *coeff, *site, *r;

    // allocate memory
    x     = (double *)malloc(n * sizeof(double));
    y     = (double *)malloc(n * sizeof(double));
    coeff = (double *)malloc((n-1) * order * sizeof(double));
    site  = (double *)malloc(m * N_SITE_BL * sizeof(double));
    r     = (double *)malloc(m * N_SITE_BL * sizeof(double));

    // generate break points and function values
	for (i = 0; i < n; i++)
    {
        x[i] = (double)i;
        y[i] = 2.0 * x[i];
    }

    // generate interpolation sites
    for (i = 0; i < m*N_SITE_BL; i++)
    {
        site[i] = (double)n * ((double)rand() / (double)RAND_MAX);
    }

    // calculate interpolation results in parallel
#pragma omp parallel for
	for (i = 0; i < N_SITE_BL; i++)
	{
		ThrInterp(i, n, x, y, order, coeff, m, site, r);
	}

    // print results
    printf("Interpolation results\n");
    for (i = 0; i < m*N_SITE_BL; i++)
    {
        printf("%lf \t => \t %lf\n", site[i], r[i]);
    }

    // release memory
    free(x);
    free(y);
    free(coeff);
    free(site);
    free(r);
}

void ThrInterp(int i_block, MKL_INT n, double *x, double *y, MKL_INT order,
               double *coeff, MKL_INT m, double *site, double *r)
{
    int errcode = 0;
    double * site_ptr = site + i_block*m;
    double * r_ptr = r + i_block*m;
    MKL_INT dorder = 1;
    DFTaskPtr task;

    // create Data Fitting task
    errcode = dfdNewTask1D( &task, n, x, DF_NON_UNIFORM_PARTITION, 1, y, DF_NO_HINT );
    if (errcode != 0) { printf("Error in dfdNewTask1D, code = %d\n", errcode); }

    // edit Data Fitting task
    errcode = dfdEditPPSpline1D( task, order, DF_PP_DEFAULT, DF_NO_BC, 0, DF_NO_IC, 0,
                                 coeff, DF_NO_HINT );
    if (errcode != 0) { printf("Error in dfdEditPPSpline1D, code = %d\n", errcode); }

    // construct spline
    errcode = dfdConstruct1D( task, DF_PP_SPLINE, DF_METHOD_STD );
    if (errcode != 0) { printf("Error in dfdConstruct1D, code = %d\n", errcode); }

    // calculate interpolation results
    errcode = dfdInterpolate1D( task, DF_INTERP, DF_METHOD_PP,
                                    m, site_ptr, DF_NON_UNIFORM_PARTITION, 1,
                                    &dorder, 0, r_ptr, DF_NO_HINT, 0 );
    if (errcode != 0) { printf("Error in dfdInterpolate1D, code = %d\n", errcode); }

    // delete Data Fitting task
    errcode = dfDeleteTask(&task);
    if (errcode != 0) { printf("Error in dfDeleteTask, code = %d\n", errcode); }
}