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
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
********************************************
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.
Dick Hendrickson
> I've included my program and the MVN random number generation
> subroutine code below ...
...
> subroutine mvn_rng_mkl(seed_in, n_mvn_obs_in, Cholesky,
> mvn_a0a1b_list)
>
> use :: mkl_vsl_type
> integer :: brng
...
> external :: vslnewstream, vsrnggaussianmv, vsldeletestream
>
...
What you haven't included is the calling protocol/argument list for the
routine vslDeleteStr that raised an error. Something like you passed an
integer instead of real or vice versa would easily cause the problem.
It's also interesting there's the reference to "vsldeletestream" in the
subroutine mvn_rng_mkl--seems peculiar.
Overall, insufficient data...
--
>
> 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
vslDeleteStream is an Intel Math Kernel Library (MKL) routine. I
suggest that you ask for assistance at the Intel MKL user forum at
http://software.intel.com/en-us/forums/intel-math-kernel-library/
--
Steve Lionel
Developer Products Division
Intel Corporation
Nashua, NH
For email address, replace "invalid" with "com"
User communities for Intel Software Development Products
http://software.intel.com/en-us/forums/
Intel Software Development Products Support
http://software.intel.com/sites/support/
My Fortran blog
http://www.intel.com/software/drfortran
It's surprising that you code got that far.
The USE statement is incorrect, so it should have failed in compilation.
However, assuming that you correct that, it's telling you what the
error is, namely, that the first argument to vslDeleteStream is wrong.
Do you have an interface for it? And for the other external procedures
in that subroutine?
If you still can't see what the error is, turn on all compiler checks.
More than that is difficult to say, because you did not include the module
mkl_vsl_type, and (most importantly) the external procedure
vslDeleteStream which was the subject of your error.
BTW, it seems that you are using free source form, however,
you have not used proper Comment lines. Is it possible that
your long declaration line is too long?
Thanks for your responses. I apologize in advance for my ignorance on
some of the things you've mentioned - I'm not a programmer, just a
stats student who was advised to use Fortran for my research :)
The subroutine (called mvn_rng_mkl) I've used is similar to ones I
have used to sample from normal and uniform distributions, which I've
done without difficulty. The 'vsldeletestream' is not a (sub)routine
that I've created, but as a previous post mentions, is part of the MKL
routine (as far as I know!). I've always included it when generating
random numbers from other distributions, and have not had problems
with it. Interestingly, if I remove the entire line "status =
vsldeletestream(stream)" from my subroutine mvn_rng_mkl, the MKL error
messages I was receiving before disappear, and I am left with the
error:
forrtl: severe (174): SIGSEGV, segmentation fault occurred
If I include a DEALLOCATE statement in my program, to clear the list
'mvn_a0a1b_list' at the end of each loop (which I think I need/want to
do), I get the additional error:
*** glibc detected *** free(): invalid next size (fast):
0x600000000001ca90 ***
*** glibc detected *** free(): invalid next size (fast):
0x600000000001cae0 ***
forrtl: severe (174): SIGSEGV, segmentation fault occurred
While I suppose this is a step in the right(!?!) direction, I'm still
not sure why I'm having such difficulty sampling from a MVN
distribution, when I am able to sample from other distributions
without a problem!
Thanks again for your comments/suggestions,
Lorna
Maybe it's your source code. You have some very long lines. I am
assuming that you are using free form source, even though all your
program statements are indented. [If you are putting TAB characters in
your source file, DON'T do that.] Perhaps some of your source
lines are being truncated? It makes more sense, and it makes it much
easier to read your source if you use line continuations and
have source lines at most 80 characters or less.
Either you can do
real :: foo, bar, baz, ....................... bozo, &
fred, barney, ...... joe_bloggs, &
larry, moe, shemp
or
real :: foo, bar, baz .................... bozo
real :: fred, barney ..........
real :: larry, moe, shemp
Otherwise the chief suspects are
1. argument mismatch between calling routine and called routine
2. array access out of bounds
Allocating an array inside a DO loop is a sign of trouble. If it is a
fixed size, then why not allocate it outside the loop
or just dimension it to the proper size?
Also, most institutions do have academic computing and stats support
people on staff.
--- e
> The subroutine (called mvn_rng_mkl) I've used is similar to ones I
> have used to sample from normal and uniform distributions, which I've
> done without difficulty. The 'vsldeletestream' is not a (sub)routine
> that I've created, but as a previous post mentions, is part of the MKL
> routine (as far as I know!). I've always included it when generating
> random numbers from other distributions, and have not had problems
> with it.
I fail to see what random number generation has to do with the error at
all. As noted, the error message refers to vsldeletestream. Not knowing
anything about that routine, I can't be of much help other than to say
that...
> While I suppose this is a step in the right(!?!) direction, I'm still
> not sure why I'm having such difficulty sampling from a MVN
> distribution, when I am able to sample from other distributions
> without a problem!
No, this doesn't sound like a step in the right direction. I'd say that
your problem likely has little or nothing to do with sampling from
random number distributions. I can't even understand why you bring up
the matter of random number sampling; I was about to say something about
that seeming to be a completely random thing to bring up, but I suppose
that's not the best way for me to expres it.
Since the error message is from vsldetestream, I'd say that the right
direction would be to read whatever documentation you have about it. I
suspect that aimlessly fiddling around with other things is likely to
achieve nothing.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
I would put it more strongly. It is extremely unlikely to have
anything to do with random number generation. From its name, I would
expect it to have to do with I/O.
>No, this doesn't sound like a step in the right direction. I'd say that
>your problem likely has little or nothing to do with sampling from
>random number distributions. I can't even understand why you bring up
>the matter of random number sampling; I was about to say something about
>that seeming to be a completely random thing to bring up, but I suppose
>that's not the best way for me to expres it.
I agree completely.
>Since the error message is from vsldetestream, I'd say that the right
>direction would be to read whatever documentation you have about it. I
>suspect that aimlessly fiddling around with other things is likely to
>achieve nothing.
I wouldn't. I would start looking for array bound errors, argument
errors, scoping errors and other such things that cause generalised
data trashing. My bet is that one of those has occurred and the
actual error message bears very little relationship to the problem.
Regards,
Nick Maclaren.
Arguments are mismatched. RTM.
see http://www.intel.com/software/products/mkl/docs/WebHelp/sf/functn_GaussianMV.html
The mvn_a0a1b_list array is one dimensional in the calling program,
but the documentation says it is two dimensional. I can not tell from
the docs which dimension comes first! Should it be (10,3) or (3,10)? I
don't know. Nor do I know enough stats to guess.
Storing 30 numbers where your array only has room for 10 trashes
something. Perhaps it is stream? I suspect that stream (which is
actually an array) is now invalid, and hence KABOOM.
--- e
>Thanks for your responses. I apologize in advance for my ignorance on
>some of the things you've mentioned - I'm not a programmer, just a
>stats student who was advised to use Fortran for my research :)
>The subroutine (called mvn_rng_mkl) I've used is similar to ones I
>have used to sample from normal and uniform distributions, which I've
>done without difficulty. The 'vsldeletestream' is not a (sub)routine
>that I've created, but as a previous post mentions, is part of the MKL
>routine (as far as I know!). I've always included it when generating
>random numbers from other distributions, and have not had problems
>with it. Interestingly, if I remove the entire line "status =
>vsldeletestream(stream)" from my subroutine mvn_rng_mkl, the MKL error
>messages I was receiving before disappear, and I am left with the
>error:
>forrtl: severe (174): SIGSEGV, segmentation fault occurred
Have you turned on all compiler checks, as suggested.
Have you provided an interface, as I suggested?
>If I include a DEALLOCATE statement in my program, to clear the list
>'mvn_a0a1b_list' at the end of each loop (which I think I need/want to
>do), I get the additional error:
> *** glibc detected *** free(): invalid next size (fast):
>0x600000000001ca90 ***
>*** glibc detected *** free(): invalid next size (fast):
>0x600000000001cae0 ***
>forrtl: severe (174): SIGSEGV, segmentation fault occurred
>While I suppose this is a step in the right(!?!) direction, I'm still
>not sure why I'm having such difficulty sampling from a MVN
>distribution, when I am able to sample from other distributions
>without a problem!
Have you put in an interface? Have you turned on all compiler checks?
I expect her best luck will be w/ the suggestion to find her uni
computing support staff combined w/ spending some time w/ the users
guide for the compiler.
--
Hello,
Thank you for both your comments - they were extremely helpful. I've
corrected the arguments so that the dimensions match, and I've moved
the "allocate" statement outside the DO loop. It looks to be running
smoothly now!
Thanks again,
Lorna
Hello,
Thanks for your suggestions. In terms of an interface for
"vsldeletestream", I'm not sure about this. I've not had to include
an interface for it with any of my other RNG routines, and my basic
understanding of it is that it's included in the MKL. However, I was
able to make some changes based on comments from another post, and
turned on the compiler checks when I tried running it. With the
changes, it's not reporting any errors, and looks to be running o.k.
Thanks again for your help!
Lorna
Hello,
Yes, I'm very unfamiliar with Fortran! I apologize if this or any of
my questions became annoying, but unfortunately I have limited
departmental support for using Fortran. Thanks so much to everyone
for all the comments and suggestions made - things look to be running
smoothly now, and I really appreciate some of the suggestions I can
keep in mind for future problems.
Cheers,
Lorna
I agree with respect to the interface suggestion, but not about compiler
checks. Knowing how to alter compiler diagnostics may not be the most
basic of programming skills, but it's an essential one to doing any
remotely serious work.
> I expect her best luck will be w/ the suggestion to find her uni
> computing support staff combined w/ spending some time w/ the users
> guide for the compiler.
Both are certainly good resources for figuring out how to turn on
compiler checks.