#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"

extern void thr_interp_(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++)
	{
        thr_interp_(&n, x, y, &order, coeff, &m, site+i*m, r+i*m);
	}

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