- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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*.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page