Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

Random numbers with vdRngGaussianMV.

Fiori
Beginner
513 Views

Hi!

I want to generate random numbers from a multivariate normal distribution with parameters mu = [3.0 5.0 2.0] and sigma = [ 16.0 8.0  4.0;  8.0 13.0 17.0; 4.0 17.0 62.0].

I wrote a mex-code using the the function vdRngGaussianMV and the results of the simulation are

>> mean(out)
ans =
    3.4558    2.9934    3.3013
>> cov(out)
ans =
  146.5081   -1.9068   -2.1787
   -1.9068  144.7461   -0.7771
   -2.1787   -0.7771  142.3955

Could you please tell me what I am doing wrong.

Thank you very much.

The code is the following:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <mkl.h>
#include "mkl_vml.h"
#include "mex.h"
#include "matrix.h"
#include "mkl_vsl.h"
#include <time.h>

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

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


/* main fucntion */
void mexFunction(int nlhs,  mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    
    int npars,  N;
    double *cov, *mean, *out;
    
    /* make pointers to input data */
    npars = (int)mxGetScalar(prhs[0]);
    N = (int)mxGetScalar(prhs[1]);
    cov = mxGetPr(prhs[2]);
    mean = mxGetPr(prhs[3]);
    
    /* make pointers to output data */
    plhs[0] = mxCreateDoubleMatrix( N, npars, mxREAL);
    out = mxGetPr(plhs[0]);
    
    
    /* call */
    normalMVN(npars, N, cov, mean, out);
    
    return 0;
    
}


double normalMVN(int npars, int N, double *cov, double *mean, double *out)
{
    
    /* This we need it for distributions */
    VSLSSTaskPtr task;
    VSLStreamStatePtr stream;
    
    char   uplo;
    int i, j, n, lda,  info;
    double *T, fiori[3];
    
    T = (double *)mxMalloc( npars*npars*sizeof( double* ) );
    
    uplo = 'L';
    n = npars;
    lda = npars;
    
    cblas_dcopy(npars*npars, cov, 1, T, 1);
    
    dpotrf( &uplo, &n, T, &lda, &info );
    if(info != 0){mexPrintf("c++ error: Cholesky failed\n\n");}
    
    vslNewStream( &stream, BRNG, SEED );
    
    vdRngGaussianMV( VSL_METHOD_DGAUSSIANMV_BOXMULLER2, stream, N, out, npars, VSL_MATRIX_STORAGE_FULL, mean, T );
  
    vslDeleteStream( &stream );
    
    /* Free memory */
    mxFree(T);
    
    return 0;
}

 

 

0 Kudos
1 Solution
Andrey_N_Intel
Employee
513 Views

Hi,

Thanks for the generated data. I reproduced covariance and mean estimates with your dataset. This indicates the bug in that version of Intel MKL. I attach C test case that was based on the your code, and the output produced by this test case. For the experiments I used Intel MKL from Intel Parallel Studio 2016 for Linux. Can you please download the latest version of the Studio and run the experiments with it?

Thanks, Andrey

 

 

View solution in original post

0 Kudos
5 Replies
Andrey_N_Intel
Employee
513 Views

Hi,

Per quick analysis of your code, it makes sense to apply a few corrections:

- replace uplo = 'L'; with uplo = 'U';

-remove '*' in sizeof( double* ) in T = (double *)mxMalloc( npars*npars*sizeof( double* ) );

For debugging purposes can you please set SEED  to a fixed number, provide value of N, and the data set you get from the generator for those parameters (I assume N is not huge, something like ~1000) ? It would help to compare the results on your and our side.

Also, which library version do you use?

Andrey

 

0 Kudos
Fiori
Beginner
513 Views

Hello.

I corrected the code as you proposed.  I set the  SEED equal to 777,  npars =3, N=5000, sigma = [ 16.0 8.0  4.0;  8.0 13.0 17.0; 4.0 17.0 62.0];,  m=[3.0 5.0 2.0]; I have installed on my pc  "Intel Parallel Studio XE 2013 SP1, Update 2 for Windows".

I run the code with the command  out = fio(npars, N, sigma, m); and the results are

>> mean(out)
ans =
    3.3099    3.4070    3.3614
>> cov(out)
ans =
   32.5436   -0.5491   -0.9751
   -0.5491   29.8680   -0.6707
   -0.9751   -0.6707   32.0381
>>

Please see the attached txt for the simulated data.

The corrected code is:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <mkl.h>
#include "mkl_vml.h"
#include "mex.h"
#include "matrix.h"
#include "mkl_vsl.h"
#include <time.h>

#define SEED 777
#define BRNG  VSL_BRNG_MCG31 // VSL basic generator to be used

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


/* main fucntion */
void mexFunction(int nlhs,  mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    
    int npars,  N;
    double *cov, *mean, *out;
    
    /* make pointers to input data */
    npars = (int)mxGetScalar(prhs[0]);
    N = (int)mxGetScalar(prhs[1]);
    cov = mxGetPr(prhs[2]);
    mean = mxGetPr(prhs[3]);
    
    /* make pointers to output data */
    plhs[0] = mxCreateDoubleMatrix( N, npars, mxREAL);
    out = mxGetPr(plhs[0]);
    
    
    /* call */
    normalMVN(npars, N, cov, mean, out);
    
    return 0;
    
}


double normalMVN(int npars, int N, double *cov, double *mean, double *out)
{
    
    char   uplo;
    int i, j, n, lda,  info;
    double *T;
    
    
    /* This we need it for distributions */
    VSLSSTaskPtr task;
    VSLStreamStatePtr stream;
    vslNewStream( &stream, BRNG, SEED );
    
    
    /* Following variables are used in Cholesky factorization subroutine */
    T =  mxMalloc( npars*npars*sizeof( double ) );
    
    uplo = 'U';
    n = npars;
    lda = npars;
    
    cblas_dcopy(npars*npars, cov, 1, T, 1);
    
    /* MKL Choelsky factorization routine call */
    dpotrf( &uplo, &n, T, &lda, &info );
    if(info != 0){mexPrintf("c++ error: Cholesky failed\n\n");}
    
    
    /*****************************************************
     * Generating random numbers
     * from multivariate normal distribution
     *****************************************************/
    vdRngGaussianMV( VSL_METHOD_DGAUSSIANMV_BOXMULLER2, stream, N, out, npars, VSL_MATRIX_STORAGE_FULL, mean, T );
    
    
    vslDeleteStream( &stream );
    
    /* Free memory */
    mxFree(T);
    
    return 0;
}

 

0 Kudos
Andrey_N_Intel
Employee
514 Views

Hi,

Thanks for the generated data. I reproduced covariance and mean estimates with your dataset. This indicates the bug in that version of Intel MKL. I attach C test case that was based on the your code, and the output produced by this test case. For the experiments I used Intel MKL from Intel Parallel Studio 2016 for Linux. Can you please download the latest version of the Studio and run the experiments with it?

Thanks, Andrey

 

 

0 Kudos
Fiori
Beginner
513 Views

Hi!

I really appreciate your help. I installed on my pc   'Intel® Parallel Studio XE Cluster Edition , 2016 for Windows'  but  I don't know how to run my mex-code. With  "Intel Parallel Studio XE 2013 SP1" I used the command

'mex -v -g phi1c.c -I"C:\Program Files (x86)\Intel\Composer XE 2013 SP1\mkl\include" -L"C:\Program Files (x86)\Intel\Composer XE 2013 SP1\mkl\lib\intel64"  -lmkl_intel_lp64 -lmkl_core -lmkl_sequential -lmkl_cdft_core -lmkl_intel_thread'

Could you please tell me if you know  what to do?

Thank you very much.

 

0 Kudos
Ying_H_Intel
Employee
513 Views

Hi GF G.

Are you asking the eqivalent folder in Intel® Parallel Studio XE Cluster Edition , 2016 for Windows? 

We changed the folder a little, for example, 

 C:\Program Files (x86)\Intel\Composer XE 2013 SP1\mkl

in 

C:\Program Files (x86)\Intel_sw_development_tools\parallel_studio_xe_2016.0.0.xx\compilers_and_libraries_2016\windows\mkl

You may replace them in your command 

Best Regards

Ying 

0 Kudos
Reply