Re: Error message from random number generation

Lorna wrote:
Hello,

I am a relatively new user to Fortran, and have pretty basic
programming skills (my research is in statistics). I was looking for
some guidance with a program I'm working on, and came across this
forum. I'm hoping someone here may be able to point me in the right
direction with the problem I'm having.

I'm trying to sample numbers from a multivariate normal distribution,
to eventually be used in a Markov Chain Monte Carlo simulation. My
program is supposed to read in some model parameter estimates from a
previous MCMC, call a subroutine to calculate the covariance between
the model parameters, and then determine the Cholesky decomposition.
The Cholesky matrix is then supposed to be passed into another
subroutine that generates MVN values. It does all of this in a loop
(as eventually I'll want it to be able to update the parameter
estimates, update the covariance and Cholesky matrices, then resample
from a slightly different MVN distribution), and each time it is
supposed to write the variance and Cholesky matrices, along with the
sampled values, to a file.

When I run my program, I seem to get the desired output written to
file (i.e. it seems to be able to calculate the variance matrix,
Cholesky matrix, and even generate MVN random numbers), however, I
also get the following error message:

MKL ERROR : Parameter 1 was incorrect on entry to vslDeleteStr
MKL ERROR : Parameter 1 was incorrect on entry to vslDeleteStr
forrtl: severe (174): SIGSEGV, segmentation fault occurred

The error messages seem to be coming from a routine called
vsldeletestr, but your code below doesn't reference or contain
a routine with that name. Possible could it be a shortened version
of vsldelerestream?

Anyhow, the obvious problem is that something is wrong with the
first argument to vsldeletestr ;).

Two possibilities to consider. Do you need an interface to
vsldeletestream? The argument, stream, is a derived type
and they often require interfaces or at least common module
use to get the types right.

Also, most of your thing are of the form
status = some_function(argument list)
but you never check the status return. At least put a
PRINT *, status
is to see if something goes wrong earlier in the program.

The SIGSEGV error usually means you are either using too much
memory or you have an address that is out of bounds. You
can look for some way to make the stack bigger (operating system
dependent, often something like ulimit does it). Or, since
an argument is believed to be bad, it's likely that internal
pointers within STREAM are not pointing to what they should.
This can cause the sigsegv. The third possibility is that
you have a subscript that is out of bounds. Look at you
compiler documentation and try turning on all of the debugging
options you can find. In particular, subscript bounds checking.

I'd worry about the MKL error first. Often when you get a
cascade of error messages, only the first one makes much sense.
The rest are just a consequence of the first.

*** Hendrickson

I've tried running a smaller program that simply calculates the
Cholesky decomposition and calls the MVN random number generator, and
that program seems to run without a problem. I'm not sure what is
causing the problem when I move to a slightly more complex scenario.
I've included my program and the MVN random number generation
subroutine code below (there are additional subroutines to be added
in, so some of the variables are not used here). If anyone has any
suggestions on what I could try, it would be greatly appreciated.

Thanks,

Lorna

****************************
program trial

real :: harvest, ratio, accept_prob, accept_prob_zi, a0t, a1t, bt,
a0t1, a1t1, bt1, za0, za1, zb, lnlike_at, lnlike_at1, lnlike_zt,
lnlike_zt1, prior_a0t, prior_a0t1, prior_a1t, prior_a1t1, prior_bt,
prior_bt1
integer :: seed_in, n_rand_obs_in, n_mvn_obs_in, r_uni_count,
r_mvn_count, sigma_estimates
real, dimension(:), pointer :: r_uni_list, mvn_a0a1b_list
real, dimension(1:1000) :: a0MCMC, a1MCMC, bMCMC
real, dimension(1:3,1:3) :: Cholesky, cholesky_matrix, VarCovar

interface
subroutine uni_rng_mkl(seed_in, n_rand_obs_in, r_uni_list)
integer, intent(in) :: seed_in, n_rand_obs_in
real, dimension(:), intent(out) :: r_uni_list
end subroutine uni_rng_mkl
end interface

interface
subroutine mvn_rng_mkl(seed_in, n_mvn_obs_in, Cholesky,
mvn_a0a1b_list)
integer, intent(in) :: seed_in, n_mvn_obs_in
real, dimension(1:3,1:3), intent(in) :: Cholesky
real, dimension(:), intent(out) :: mvn_a0a1b_list
end subroutine mvn_rng_mkl
end interface

interface
subroutine covariance(a0MCMC, a1MCMC, bMCMC, Cholesky, VarCovar)
real, dimension(:), intent(in) :: a0MCMC, a1MCMC, bMCMC
real, dimension(1:3,1:3), intent(out) :: Cholesky, VarCovar
end subroutine covariance
end interface

n_rand_obs_in = 10
n_mvn_obs_in = 10
r_uni_count = 10
r_mvn_count = 10
sigma_estimates = 2

open(12, file = 'Study2_abestimates_Trial6.txt', status = 'old')
do i = 1, 1000
read(12,*) a0MCMC(i), a1MCMC(i), bMCMC(i)
end do
close(12)

a0t = 0.5e0
a1t = 0.5e0
bt = 2.0e0

call random_number(harvest)
seed_in = aint(harvest * 2.0e0**31)

allocate (r_uni_list(n_rand_obs_in))
call uni_rng_mkl(seed_in, n_rand_obs_in, r_uni_list)

open(13, file='Cholesky_estimates_X2.txt', status = 'old')

do m = 1, sigma_estimates

call covariance(a0MCMC, a1MCMC, bMCMC, Cholesky, VarCovar)

write(13, *) Cholesky
write(13, *) VarCovar

call random_number(harvest)
seed_in = aint(harvest * 2.0e0**31)

allocate (mvn_a0a1b_list(n_mvn_obs_in))
call mvn_rng_mkl(seed_in, n_mvn_obs_in, Cholesky, mvn_a0a1b_list)

do i = 1, 10
write(13, *) mvn_a0a1b_list(i)
end do

end do

close(13)

end program trial
***************************************
subroutine mvn_rng_mkl(seed_in, n_mvn_obs_in, Cholesky,
mvn_a0a1b_list)

use :: mkl_vsl_type
integer :: brng

integer, intent(in) :: seed_in, n_mvn_obs_in
real, dimension(1:3,1:3), intent(in) :: Cholesky
real, dimension(:), intent(out) :: mvn_a0a1b_list
real, dimension(1:3) :: mean_vector

type(VSL_STREAM_STATE) :: stream
integer :: status, methodrng, storagerng

external :: vslnewstream, vsrnggaussianmv,
vsldeletestream

mean_vector(1:3) = 0.0e0
methodrng = vsl_method_sgaussianmv_boxmuller
storagerng = vsl_matrix_storage_full

brng = vsl_brng_mcg31
status = vslnewstream(stream, brng, seed_in)
status = vsrnggaussianmv (methodrng, stream,
n_mvn_obs_in, mvn_a0a1b_list, 3, storagerng, mean_vector, Cholesky)

status = vsldeletestream(stream)

end subroutine
********************************************
.

Quantcast