diff --git a/DESCRIPTION b/DESCRIPTION
index 5aef20a..652c579 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: TruncatedNormal
 Type: Package
 Title: Truncated Multivariate Normal and Student Distributions
-Version: 2.0
+Version: 2.0.600
 Date: 2019-06-07
 Authors@R: c(person(given="Zdravko", family="Botev", role = "aut", email = "botev@unsw.edu.au", comment = c(ORCID = "0000-0001-9054-3452")), person(given="Leo", family="Belzile", role = c("aut", "cre"), email = "belzilel@gmail.com", comment = c(ORCID = "0000-0002-9135-014X")))
 Description: A collection of functions to deal with the truncated univariate and multivariate normal and Student distributions, described in Botev (2017) <doi:10.1111/rssb.12162> and Botev and L'Ecuyer (2015) <doi:10.1109/WSC.2015.7408180>.
diff --git a/NEWS b/NEWS
index 64ffb64..d3f5095 100644
--- a/NEWS
+++ b/NEWS
@@ -1,5 +1,9 @@
 TruncatedNormal changelog
 
+version 2.1
+
+* bug fix: Cholesky decomposition with permutation now checks the arguments to ensure lb less than ub, to avoid segfault errors. In case of degenerate bounds, the conditional distribution is used with fixed components for the degenerate variables (issue #2)
+
 version 2.0
 
 * Can specify mean vector directly in the functions
diff --git a/R/RcppExports.R b/R/RcppExports.R
index b040ed2..2dfa4db 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -50,7 +50,7 @@
 #' @useDynLib TruncatedNormal, .registration = TRUE
 #' @importFrom Rcpp evalCpp
 #' @export
-lnNpr <- function(a, b, check = FALSE) {
+lnNpr <- function(a, b, check = TRUE) {
     .Call(`_TruncatedNormal_lnNpr`, a, b, check)
 }
 
diff --git a/R/tmvnorm.R b/R/tmvnorm.R
index 17558ed..74e61ba 100644
--- a/R/tmvnorm.R
+++ b/R/tmvnorm.R
@@ -163,10 +163,35 @@ rtmvnorm <- function(n, mu, sigma, lb, ub){
   if(missing(ub)){
     ub <- rep(Inf, d) 
   }
-  if(n == 1){
-    as.vector(mvrandn(l = lb, u = ub, Sig = sigma, n = n, mu = mu)) 
+  stopifnot(length(lb) == length(ub), length(lb) == d, lb <= ub)
+  if(!any((ub - lb) < 1e-10)){
+    if(n == 1){
+      as.vector(mvrandn(l = lb, u = ub, Sig = sigma, n = n, mu = mu)) 
+    } else{
+      t(mvrandn(l = lb, u = ub, Sig = sigma, n = n, mu = mu))
+    }
   } else{
-    t(mvrandn(l = lb, u = ub, Sig = sigma, n = n, mu = mu))
+    warning("Some variables have a degenerate distribution.")
+    ind <- which((ub - lb) >= 1e-10)
+    # check covariance matrix
+    stopifnot(isSymmetric(sigma), all(eigen(sigma, only.values = TRUE)$value > 0))
+    # compute conditional Gaussian
+      schurcomp <- function(sigma, ind) {
+        stopifnot(c(length(ind) > 0, ncol(sigma) - length(ind) > 0))
+        sigma[ind, ind, drop = FALSE] - sigma[ind, -ind, drop = FALSE] %*%
+          solve(sigma[-ind, -ind, drop = FALSE]) %*% sigma[-ind, ind, drop = FALSE]
+      }  
+      sigmap <- schurcomp(sigma, ind)
+      mup <- c(mu[ind] + sigma[ind, -ind, drop = FALSE] %*% solve(sigma[-ind, -ind, drop = FALSE]) %*% (lb[-ind] - mu[-ind]))
+      #
+      res <- matrix(0, nrow = n, ncol = d)
+      res[, -ind] <- rep(lb[-ind], each = n)
+    if(n == 1){
+      res[, ind] <- as.vector(mvrandn(l = lb[ind], u = ub[ind], Sig = sigmap, n = n, mu = mup)) 
+      return(as.vector(res))
+    } else{
+      res[, ind] <- t(mvrandn(l = lb[ind], u = ub[ind], Sig = sigmap, n = n, mu = mup))
+    }
   }
 }
 
diff --git a/R/tmvt.R b/R/tmvt.R
index 7410609..8616822 100644
--- a/R/tmvt.R
+++ b/R/tmvt.R
@@ -187,10 +187,39 @@ rtmvt <- function(n, mu, sigma, df, lb, ub){
   if(missing(ub)){
     ub <- rep(Inf, d) 
   }
-  if(n == 1){
-    as.vector(mvrandt(l = lb, u = ub, Sig = sigma, df = df, n = n, mu = mu))
+  stopifnot(length(lb) == length(ub), length(lb) == d, lb <= ub)
+  if(!any((ub - lb) < 1e-10)){
+    
+    if(n == 1){
+      as.vector(mvrandt(l = lb, u = ub, Sig = sigma, df = df, n = n, mu = mu))
+    } else{
+      t(mvrandt(l = lb, u = ub, Sig = sigma, df = df, n = n, mu = mu))
+    }
   } else{
-    t(mvrandt(l = lb, u = ub, Sig = sigma, df = df, n = n, mu = mu))
+    warning("Some variables have a degenerate distribution.")
+    ind <- which((ub - lb) >= 1e-10)
+    # check covariance matrix
+    stopifnot(isSymmetric(sigma), all(eigen(sigma, only.values = TRUE)$value > 0))
+    # compute conditional Gaussian
+    schurcomp <- function(sigma, ind) {
+      stopifnot(c(length(ind) > 0, ncol(sigma) - length(ind) > 0))
+      sigma[ind, ind, drop = FALSE] - sigma[ind, -ind, drop = FALSE] %*%
+        solve(sigma[-ind, -ind, drop = FALSE]) %*% sigma[-ind, ind, drop = FALSE]
+    }  
+    sigmap <- c(df + t(lb[-ind] - mu[-ind]) %*% solve(sigma[-ind, -ind, drop = FALSE]) %*% (lb[-ind] - mu[-ind])) /
+      (df + d - length(ind)) * schurcomp(sigma, ind)
+    mup <- c(mu[ind] + sigma[ind, -ind, drop = FALSE] %*% solve(sigma[-ind, -ind, drop = FALSE]) %*% (lb[-ind] - mu[-ind]))
+    # matrix to store results
+    res <- matrix(0, nrow = n, ncol = d)
+    res[, -ind] <- rep(lb[-ind], each = n)
+    if(n == 1){
+      res[, ind] <-  as.vector(mvrandt(l = lb[ind], u = ub[ind], Sig = sigmap, 
+                                       df = df + d - length(ind), n = n, mu = mup))
+      res <- as.vector(res)
+    } else{
+      res[, ind] <- t(mvrandt(l = lb[ind], u = ub[ind], Sig = sigmap, df = df + d - length(ind), n = n, mu = mup))
+    }
+  return(res)
   }
 }
 
diff --git a/src/TruncatedNormal.so b/src/TruncatedNormal.so
index 307ce18..df47a98 100755
Binary files a/src/TruncatedNormal.so and b/src/TruncatedNormal.so differ
diff --git a/src/lnNpr_cholperm_Phinv.o b/src/lnNpr_cholperm_Phinv.o
index 52924e2..7156c7c 100644
Binary files a/src/lnNpr_cholperm_Phinv.o and b/src/lnNpr_cholperm_Phinv.o differ