GaussianMV parameter error

[jaureus]

Hello, I am trying to work with the random number generator of the mutlivariate Gaussian distribution in FORTAN. TO this end, I have compiled my MKL library, and am attempting to get the random numbers.


1. I have several questions. One it is hard to tell from the documentation. whether ytou need the Cholesky decompostion, or the full covariance matrix. which is required


2. When I attempt to run the below program, I recive the following error: Parameter 4 was incorrect on entry to vsRngGaussia

Does anyone have any suggestions as to what is wrong with my code. Please also note that I have tried using the Cholesky decomposition code from LAPACK but it does not seem to work.

program mcnet

implicit none

integer :: i, j

real :: harvest1, harvest2, sd

integer :: dimen, seed_in1, seed_in2, n_rand_obs_in

real , dimension (:, :), allocatable :: newidentity, sigmanew, meanvectnorm, meanvectnormtrans, meanvect, meanvecttrans, &

meanvectold, meanvecttransold, Q, r_mvgaussian_list1, r_mvgaussian_list2, sigmanewchol, newidentity2

interface

subroutine gaussianmv_rng_mkl(seed_in, n_rand_obs_in, r_mvgaussian_list, dimen, mean, tcv)

implicit none

integer, intent(in) :: seed_in, n_rand_obs_in, dimen

real, dimension(1:n_rand_obs_in, 1:dimen), intent(out) :: r_mvgaussian_list

real, dimension(1:1, 1:dimen), intent(in) :: mean

real, dimension(1:dimen, 1:dimen), intent(in) :: tcv

end subroutine gaussianmv_rng_mkl

end interface

! allowing for easy changes in dimensionality

dimen = 3

n_rand_obs_in = 3

allocate(newidentity(1:dimen, 1:dimen))

!allocate(sigmanew(1:dimen, 1:dimen))

allocate(meanvectnorm(1:dimen, 1:1))

allocate(meanvectnormtrans(1:1, 1:dimen))

!allocate(meanvect(1:dimen, 1:1))

!allocate(meanvecttrans(1:1, 1:dimen))

!allocate(meanvectold(1:dimen, 1:1))

!allocate(meanvecttransold(1:1, 1:dimen))

!allocate(sigmanewchol(1:dimen, 1:dimen))

allocate(newidentity2(1:dimen, 1:dimen))

allocate(r_mvgaussian_list2(1:n_rand_obs_in, 1:dimen) )

allocate(Q(1:dimen, 1:1))

!allocate(tcv(1:dimen, 1:dimen))

newidentity = 1.0

!scaling factor as per algoritm

sd = (2.38)**2/dimen

meanvect = 0.0

meanvectnorm = 0.0

do j = 1, dimen

! identity matrix times scaling factor as per Roberts and Rosthenthal (2009) gives innitial value of covariance matrix

newidentity( j, j) = 1.0*(0.01)/dimen

end do

print *, newidentity

call random_number(harvest1)

seed_in1 = aint(harvest1 * 2.0e0**31)

meanvectnorm(1,1) = 0.2

meanvectnorm(2,1) = 0.3

meanvectnorm(3,1) = 0.001

!print *, meanvectnorm

meanvectnormtrans = transpose(meanvectnorm)

! call spotrf('U', dimen, newidentity)

!call potrf( newidentity, 'U')

! note that since cholesky decomposition of newidentity does not change, since newidentity does not change one only needs to derive it once

newidentity2 = 0.0

do i = 1, dimen

do j = 1, dimen

if (i .ge. j) then

newidentity2(i,j) = newidentity(i,j)

end if

end do

end do

print *, newidentity2

print *, "test"

call gaussianmv_rng_mkl(seed_in1, n_rand_obs_in, r_mvgaussian_list2, dimen, meanvectnormtrans, newidentity2)

!subroutine gaussianmv_rng_mkl(seed_in, n_rand_obs_in, r_mvgaussian_list, dimen, mean, tcv)

do i = 1, n_rand_obs_in

do j = 1, dimen

print *, r_mvgaussian_list2(i, j)

end do

end do

end program


which in turns calls the subroutine that I have created to use mutlivariate random numbers,

subroutine gaussianmv_rng_mkl(seed_in, n_rand_obs_in, r_mvgaussian_list, dimen, mean, tcv)

use :: mkl_vsl_type

integer :: brng

integer, intent(in) :: seed_in, n_rand_obs_in

real, pointer, intent(out) :: r_mvgaussian_list

type(VSL_STREAM_STATE) :: stream

integer :: status, methodrng

external :: vslnewstream, vsrnggausianmv, vsldeletestream

methodrng = VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER

mstorage = VSL_MATRIX_STORAGE_FULL

brng = vsl_brng_mcg31

status = vslnewstream(stream, brng, seed_in)

status = vsrnggaussianmv( methodrng, stream, n_rand_obs_in, r_mvgaussian_list, dimen, mstorage, mean, tcv)

status = vsldeletestream(stream)

end subroutine

[yuriisig]

see ....\mkl\examples\vslf

[mecej4]

You have a large number of errors with respect to argument types in your sample code. To catch such errors at compile time rather than experiencing run-time aborts, use IMPLICIT NONE, use all available interface modules, and use /warn:all as a compiler option. Finally, note that one of the worst things you can do is to give the compiler an incorrect interface block.

The documentation clearly states that you need to supply the Cholesky factor T of the variance-covariance matrix as an argument.

There is also a functioning example, vsrnggaussianmv.f, provided in the MKL examples directory. The example shows the call the call to the Cholesky factorization routine spotrf.