#include <math.h>
#include <stdio.h>
#include <mkl.h>
#include <time.h>

//#define SEED time(NULL)
#define SEED 777
#define BRNG    VSL_BRNG_MCG31 // VSL basic generator to be used

int normalMVN(int npars, int N, double *cov, double *mean, double *out);

#define DIM 3
#define NN  5000


/* main fucntion */
int main()
{
	int i, j;
	double mean[] = { 3.0, 5.0, 2.0 };
	double cov[] = { 16.0, 8.0, 4.0, 8.0, 13.0, 17.0, 4.0, 17.0, 62.0 };
	double out[DIM*NN];

	normalMVN(DIM, NN, cov, mean, out);

	for (i = 0; i < NN; i++)
	{
		for (j = 0; j < DIM; j++)
		{
			printf("  %1.4lf", out[i*DIM+j]);
		}
		printf("\n");
	}

	return 0;
}


int normalMVN(int npars, int N, double *cov, double *mean, double *out)
{
	/* This we need it for distributions */
	VSLStreamStatePtr stream;

	char   uplo;
	int i, j, n, lda, info;
	double *T, fiori[3];

	T = (double *)mkl_malloc(npars*npars*sizeof(double), 64);

	uplo = 'U';
	n = npars;
	lda = npars;

	cblas_dcopy(npars*npars, cov, 1, T, 1);

	dpotrf(&uplo, &n, T, &lda, &info);
	if (info != 0) { printf("c++ error: Cholesky failed\n\n"); }

	vslNewStream(&stream, BRNG, SEED);
	vdRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER2, stream, N, out, npars, VSL_MATRIX_STORAGE_FULL, mean, T);
	vslDeleteStream(&stream);

	/* Free memory */
	mkl_free(T);

	return 0;
}