> dmvnorm <- function (x, mean, sigma, log = FALSE) {
+     if (is.vector(x)) {
+         x <- matrix(x, ncol = length(x))
+     }
+     if (missing(mean)) {
+         mean <- rep(0, length = ncol(x))
+     }
+     if (missing(sigma)) {
+         sigma <- diag(ncol(x))
+     }
+     if (NCOL(x) != NCOL(sigma)) {
+         stop("x and sigma have non-conforming size")
+     }
+     if (!isSymmetric(sigma, tol = sqrt(.Machine$double.eps),
+         check.attributes = FALSE)) {
+         stop("sigma must be a symmetric matrix")
+     }
+     if (length(mean) != NROW(sigma)) {
+         stop("mean and sigma have non-conforming size")
+     }
+     distval <- mahalanobis(x, center = mean, cov = sigma)
+     logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
+     logretval <- -(ncol(x) * log(2 * pi) + logdet + distval)/2
+     if (log)
+         return(logretval)
+     exp(logretval)
+ }
> 
> m <- 3
> sigma <- diag(3)
> sigma[2,1] <- sigma[1, 2] <- 3/5 ; sigma[3,1] <- sigma[1, 3] <- 1/3
> sigma[3,2] <- sigma[2, 3] <- 11/15
> adaptIntegrate(dmvnorm, lower=rep(-0.5, m), upper=c(1,4,2),
+                         mean=rep(0, m), sigma=sigma, log=FALSE,
+                maxEval=10000)
$integral
[1] 0.3341125

$error
[1] 4.185435e-06

$functionEvaluations
[1] 10065

$returnCode
[1] 0