r - Draw from a conditional multivariate normal distribution in fortran -

i trying write fortran subroutine draw subsample multivariate normal distribution conditional on state of other subspace. basically:

(x1, x2)' ~ n( (mu1, mu2)', sigma)

where covariance matrix sigma can partitioned in 4 submatrices

sigma=( s11, s12; s21, s22)

textbooks & wikipedia tell me conditional distribution of x1 on x2=a is:

x1|x1=a ~ n( mu, sigma*)

where

mu = mu1 + s12 * s22^-1 * (a - mu2)

sigma* = s11 - s12 * s22^-1 * s21

when writing in r works charm. in fortran not much.

subroutine dcondmvnorm ( didx, ndraw, sigma, nsigma, mu, tmcurr)  implicit none  integer            :: i, nsigma, ndraw, info integer            :: didx(ndraw), nidx(nsigma-ndraw), allidx(nsigma) logical            :: idxmask(nsigma) double precision   :: sigma11(ndraw, ndraw), sigma22(nsigma-ndraw,nsigma-ndraw) double precision   :: sigma(nsigma,nsigma) double precision   :: sigma12minv22(nsigma-ndraw,ndraw), isigma22(nsigma-ndraw,nsigma-ndraw) double precision   :: randnums(ndraw), dummy1(ndraw), meandiff(nsigma-ndraw ) double precision   :: tmcurr(nsigma), mu(nsigma)  ! create indeces _not_ draw (nidx) idxmask = .false. idxmask(didx) = .true. allidx = (/ (i, i=1, nsigma) /) nidx = pack( allidx, .not. idxmask)  sigma11 = sigma( didx, didx) sigma22 = sigma( nidx, nidx) isigma22 =0.0d0 = 1, nsigma-ndraw   isigma22(i,i) = 1.0d0 end call dposv( 'u', nsigma-ndraw,nsigma-ndraw, sigma22, nsigma-ndraw, isigma22, nsigma-ndraw, info) call dgemm ( 'n', 'n', ndraw, nsigma-ndraw, nsigma-ndraw, 1.0d0, sigma(didx,nidx), ndraw, &    isigma22, nsigma-ndraw, 0.0d0, sigma12minv22, ndraw )  call dgemm ( 'n', 'n', ndraw, ndraw, nsigma-ndraw, -1.0d0, sigma12minv22, ndraw, &    sigma(nidx,didx), nsigma-ndraw, +1.0d0, sigma11, ndraw) call dpotrf( 'u', ndraw, sigma11, ndraw, info) = 1, ndraw-1   sigma11(i+1:ndraw,i) = 0.0d0 end ! sigma11 cholesky decomposition of matrix sigma* meandiff = tmcurr(nidx) - mu(nidx) call dgemv( 'n', ndraw, nsigma-ndraw, 1.0d0, sigma12minv22, ndraw, meandiff, 1, 0.0d0, dummy1(1), 1)  ! sorry, 1 self written , returns ndraw random numbers i.i.d. n(0,1) using marsaglia's algorithm call getzig(randnums, ndraw) call dgemv( 'n', ndraw, ndraw, 1.0d0, sigma11, ndraw, randnums(1), 1, 1.0d0, dummy1(1), 1)  tmcurr(didx) = dummy1 end subroutine dcondmvnorm 

so build (it part of larger module working on) call r using

covmat <- diag(4) covmat[1:3,2:4] <- covmat[1:3,2:4] + diag(3)*.5 covmat[2:4,1:3] <- covmat[2:4,1:3] + diag(3)*.5 covmat[3:4,1:2] <- covmat[3:4,1:2] + diag(2)*.2 covmat[1:2,3:4] <- covmat[1:2,3:4] + diag(2)*.2 library(mass) dyn.load("tm_updater.so") testmat2 <- matrix(na,0,4) (a in seq(500) ){   y <- mvrnorm(1,rep(0,2), covmat[3:4,3:4])   x <- .fortran("dcondmvnorm", as.integer(c(1,2)),as.integer(2), covmat, as.integer(4), c(0.0,0.0,0.0,0.0), c(0.0,0.0,y))[[6]]   testmat2 <- rbind(testmat2, c(x[1:2],y) ) } dyn.unload("tm_updater.so") cov(testmat2) 

and returns

> cov(testmat2)         [,1]      [,2]      [,3]        [,4] [1,] 1.179618573 0.4183372 0.1978489 0.002156081 [2,] 0.418337156 0.8317497 0.4891746 0.204091537 [3,] 0.197848928 0.4891746 0.9649001 0.498660858 [4,] 0.002156081 0.2040915 0.4986609 1.032272666 

clearly, covariance of [1,1] high , way no matter how (or how long) run it. missing? covariance matrix calculated fortran matches 1 calculated hand, means... issues different accuracies?

plus there's weirdness dgemv need give exact starting address (see last call dgemv) of vector y (as called in documentary) in order get

y := alpha *x + beta * y, beta != 0 

any appreciated!

i feel embarrassed, future reference blunder shall remain available see.

the problem converting vector of i.i.d. n(0,1) random numbers target multivariate normal. checking textbooks need cholesky decomposition of covariance matrix s

s = aa'

note lower triangular matrix interested in, not upper calculated.

solution: in last call dgemv change 'n' 't' or calculate 'l' triangle in call dposv , rewrite zeroing out of upper triangle in following lines.