diff --git a/DESCRIPTION b/DESCRIPTION
index f110339..9a10202 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: PrestoGP
 Type: Package
 Title: Penalized Regression for Spatio-Temporal Outcomes via Gaussian Processes 
-Version: 0.2.0.9050
+Version: 0.2.0.9053
 Authors@R: c(
     person(given = "Eric",
            family = "Bair",
@@ -21,7 +21,7 @@ Authors@R: c(
 	   role = "aut"))
 Description: Simultaneous variable seletion and estimation of LUR models with spatiotemporally correlated errors that is scalable for big data.
 Depends:
-    R (>= 3.5.0)
+    R (>= 4.1.0)
 LinkingTo:
     Rcpp, RcppArmadillo, BH
 Imports:
@@ -67,3 +67,4 @@ Suggests:
     testthat (>= 3.0.0)
 Config/testthat/edition: 3
 URL: https://niehs.github.io/PrestoGP/, https://github.com/NIEHS/PrestoGP
+Config/roxygen2/version: 8.0.0
diff --git a/NEWS.md b/NEWS.md
index 8114166..9ba59ad 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,28 @@
+# PrestoGP 0.2.0.9053 (2026-5-26)
+
+## BUG FIXES
+
+* Fixed a bug where specifying the beta.hat parameter to `prestogp_fit`
+  caused an error due to a faulty length check.
+
+# PrestoGP 0.2.0.9052 (2026-1-20)
+
+## BUG FIXES
+
+* Created local versions of several GPvecchia utility functions (e.g.,
+  `vecchia_likelihood`, `vecchia_prediction`) that are robust to numerically
+  singular U/V matrices.
+
+* Modified `vecchia_prediction` and `vecchia_Mprediction` to avoid crashes
+  when the estimated V matrix is numerically singular.
+
+# PrestoGP 0.2.0.9051 (2026-1-14)
+
+## NEW FEATURES
+
+* Added a parameter `omp.cores` to `prestogp_fit` to specify the number of
+  cores used by OMP when computing the U matrix.
+
 # PrestoGP 0.2.0.9050 (2025-11-13)
 
 ## BUG FIXES
diff --git a/R/Log_Likelihood.R b/R/Log_Likelihood.R
index b9692ed..be8175e 100644
--- a/R/Log_Likelihood.R
+++ b/R/Log_Likelihood.R
@@ -59,6 +59,20 @@ negloglik.full <- function(logparams, d, y, param.seq) {
   -1 * mvtnorm::dmvnorm(y, rep(0, N), cov.mat, log = TRUE)
 }
 
+vecchia_likelihood <- function(z, vecchia.approx, covparms, nuggets,
+  covmodel = "matern") {
+  if (vecchia.approx$cond.yz == "zy")
+    warning("cond.yz='zy' will produce a poor likelihood approximation. Use 'SGV' instead.")
+  removeNAs <- getFromNamespace("removeNAs", "GPvecchia")
+  removeNAs()
+  U.obj <- createU(vecchia.approx, covparms, nuggets, covmodel)
+  res <- try(junk <- vecchia_likelihood_U(z, U.obj), silent = TRUE)
+  if (inherits(res, "try-error")) {
+    junk <- -1 * Inf
+  }
+  junk
+}
+
 #' Evaluation of the multivariate Vecchia likelihood
 #'
 #' This function is used to evaluate the multivariate Vecchia likelihood.
@@ -118,6 +132,30 @@ vecchia_Mlikelihood <- function(z, vecchia.approx, covparams) {
 }
 
 
+vecchia_likelihood_U <- function(z, U.obj) {
+  U <- U.obj$U
+  latent <- U.obj$latent
+  zord <- z[U.obj$ord.z]
+  const <- sum(!latent) * log(2 * pi)
+  z1 <- Matrix::crossprod(U[!latent, ], zord)
+  quadform.num <- sum(z1^2)
+  logdet.num <- -2 * sum(log(Matrix::diag(U)))
+  if (sum(latent) == 0) {
+    logdet.denom <- quadform.denom <- 0
+  } else {
+    U.y <- U[latent, ]
+    z2 <- as.numeric(U.y %*% z1)
+    V.ord <- U2V(U.obj)
+    z3 <- Matrix::solve(V.ord, rev(z2), system = "L")
+    quadform.denom <- sum(z3^2)
+    logdet.denom <- -2 * sum(log(Matrix::diag(V.ord)))
+  }
+  neg2loglik <- logdet.num - logdet.denom + quadform.num - quadform.denom +
+    const
+  loglik <- -neg2loglik / 2
+  loglik
+}
+
 ##############################################################################
 ### Flexible Multivariate Matern Negative Loglikelihood Function ###########
 
diff --git a/R/PrestoGP_CreateU_Multivariate.R b/R/PrestoGP_CreateU_Multivariate.R
index 3e4e35f..4fb7fc6 100644
--- a/R/PrestoGP_CreateU_Multivariate.R
+++ b/R/PrestoGP_CreateU_Multivariate.R
@@ -103,7 +103,7 @@ create_param_sequence <- function(P, ns = 1) {
   }
 
   param.sequence.begin <- c(1, P + 1, seq(P * (ns + 1) + 1, length = 3, by = P))
-  param.sequence.end <- c(P, ns * P, P, P, nk) %>% cumsum()
+  param.sequence.end <- c(P, ns * P, P, P, nk) |> cumsum()
   param.sequence <- cbind(param.sequence.begin, param.sequence.end)
 
   param.sequence
@@ -470,7 +470,7 @@ vecchia_Mspecify <- function(locs.list, m, locs.list.pred = NULL,
     ord.pred = ordering.pred, cond.yz = "SGV", conditioning = "NN",
     P = P, ondx = ondx, dist.func = dist.func,
     dist.func.code = dist.func.code, q.list = q.list,
-    n.neighbors = m
+    n.neighbors = m, n.cores = -1
   )
 }
 
@@ -609,7 +609,8 @@ createUMultivariate <- function(vec.approx, params, cov_func = NULL) {
     cur.qys <- do.call(cbind, q.list$q.y)
     cur.qzs <- do.call(cbind, q.list$q.z)
     # browser()
-    U <- createU_helper_mat(olocs, ondx, cur.qys, cur.qzs, vijs, aijs, full_const, nugget, sig2, uvec)
+    U <- createU_helper_mat(olocs, ondx, cur.qys, cur.qzs, vijs, aijs,
+      full_const, nugget, sig2, uvec, vec.approx$n.cores)
     # U <- sparseMatrix(i=U1[1,], j=U1[2,], x = U1[3,], triangular = TRUE)
   }
 
diff --git a/R/PrestoGP_Model.R b/R/PrestoGP_Model.R
index 4495999..17e7cb9 100644
--- a/R/PrestoGP_Model.R
+++ b/R/PrestoGP_Model.R
@@ -49,6 +49,8 @@
 #' @slot nscale The number of scale parameters in the model.
 #' @slot common_scale Do all columns of locs have the same scale parameter?
 #' @slot param_sequence Records the indices of the various Matern parameters.
+#' @slot omp_cores Number of cores used by OMP to compute the U matrix.
+#' @slot logparams Transformed version of the Matern parameters.
 #' See \code{\link{create_param_sequence}}.
 #'
 #' @seealso \code{\link{VecchiaModel-class}}, \code{\link{FullModel-class}},
@@ -97,6 +99,7 @@ PrestoGPModel <- setClass("PrestoGPModel",
     nscale = "numeric", # the number of scale parameters
     common_scale = "logical", # is there a common scale parameter?
     param_sequence = "matrix", # maps the indices of the various Matern parameters
+    omp_cores = "numeric", # Number of cores used by OMP to compute the U matrix
     logparams = "numeric" # transformed version of the Matern parameters
   )
 )
@@ -171,7 +174,7 @@ setGeneric(
     optim.control = list(trace = 0, reltol = 1e-3, maxit = 5000),
     penalty = c("lasso", "relaxed", "MCP", "SCAD"), alpha = 1,
     family = c("gaussian", "binomial"), nfolds = 10, foldid = NULL,
-    parallel = FALSE, cluster = NULL, adaptive = FALSE) {
+    omp.cores = -1, parallel = FALSE, cluster = NULL, adaptive = FALSE) {
     standardGeneric("prestogp_fit")
   }
 )
@@ -827,6 +830,9 @@ setMethod(
 #' what fold each observation should be assigned to in the cv.glmnet
 #' cross-validation. See \code{\link[glmnet]{cv.glmnet}} and
 #' \code{\link[ncvreg]{cv.ncvreg}}.
+#' @param omp.cores Number of parallel cores to be used by OMP when computing
+#' the U matrix in the Vecchia approximation. Defaults to -1 (which causes
+#' all available cores to be used).
 #' @param parallel Should parallel "foreach" be used to speed up the model
 #' fitting procedure where possible? Defaults to FALSE. Specifically,
 #' parallelization will be used for imputation and fitting the cv.glmnet
@@ -834,6 +840,9 @@ setMethod(
 #' only applies to glmnet models (where penalty="lasso" or
 #' penalty="relaxed"). Models using ncvreg (where penalty="MCP" or
 #' penalty="SCAD") require a cluster argument for parallelization (see below).
+#' Note that parallelization will be used to compute the Vecchia
+#' approximation even if this parameter is set to FALSE (assuming that OMP is
+#' available on the system). To disable this, set omp.cores=1.
 #' @param cluster A cluster for running cv.ncvreg in parallel. See
 #' \code{\link[ncvreg]{cv.ncvreg}} and \code{\link[parallel]{makeCluster}}.
 #' This must be specified to run cv.ncvreg in parallel. It is ignored for
@@ -947,9 +956,19 @@ setMethod(
     quiet = FALSE, verbose = FALSE, optim.method = "Nelder-Mead",
     optim.control = list(trace = 0, reltol = 1e-3, maxit = 5000),
     penalty = c("lasso", "relaxed", "MCP", "SCAD"), alpha = 1,
-    family = c("gaussian", "binomial"),
-    nfolds = 10, foldid = NULL, parallel = FALSE, cluster = NULL,
+    family = c("gaussian", "binomial"), nfolds = 10, foldid = NULL,
+    omp.cores = -1, parallel = FALSE, cluster = NULL,
     adaptive = FALSE) {
+    if (!is.numeric(omp.cores)) {
+      stop("omp.cores must be a positive integer >= 1")
+    }
+    if (length(omp.cores) != 1) {
+      stop("omp.cores must be a positive integer >= 1")
+    }
+    if (omp.cores < -1 | omp.cores == 0 | omp.cores %% 1 != 0) {
+      stop("omp.cores must be a positive integer >= 1")
+    }
+    model@omp_cores <- omp.cores
     penalty <- match.arg(penalty)
     model@penalty <- penalty
     family <- match.arg(family)
@@ -994,8 +1013,7 @@ setMethod(
       if (!is.vector(beta.hat) | !is.numeric(beta.hat)) {
         stop("beta.hat parameter must be a numeric vector")
       }
-      if (length(beta.hat) != (ncol(model@X_train) +
-            length(model@locs_train))) {
+      if (length(beta.hat) != (ncol(model@X_train) + 1)) {
         stop("Length of beta.hat must match the number of predictors")
       }
       beta.hat <- as.matrix(beta.hat)
diff --git a/R/PrestoGP_Multivariate_Vecchia.R b/R/PrestoGP_Multivariate_Vecchia.R
index 9e60fdc..39c662c 100644
--- a/R/PrestoGP_Multivariate_Vecchia.R
+++ b/R/PrestoGP_Multivariate_Vecchia.R
@@ -542,6 +542,7 @@ setMethod("specify", "MultivariateVecchiaModel", function(model) {
   locs <- model@locs_train
   locs.scaled <- scale_locs(model, locs)
   model@vecchia_approx <- vecchia_Mspecify(locs.scaled, model@n_neighbors)
+  model@vecchia_approx$n.cores <- model@omp_cores
   if (!model@common_scale) {
     olocs.scaled <- model@vecchia_approx$locsord
     for (i in seq_along(locs)) {
diff --git a/R/PrestoGP_Util_Functions.R b/R/PrestoGP_Util_Functions.R
index 1481cb8..6137485 100644
--- a/R/PrestoGP_Util_Functions.R
+++ b/R/PrestoGP_Util_Functions.R
@@ -226,6 +226,44 @@ revMat <- function(mat) {
   mat.out
 }
 
+vecchia_prediction <- function(z, vecchia.approx, covparms, nuggets,
+  var.exact, covmodel = "matern", return.values = "all") {
+  removeNAs <- getFromNamespace("removeNAs", "GPvecchia")
+  removeNAs()
+  U.obj <- createU(vecchia.approx, covparms, nuggets, covmodel)
+  V.ord <- U2V(U.obj)
+  if (length(U.obj$zero.nugg) > 0)
+    warning("Rows/cols of V have been removed for data with zero noise")
+  vecchia_mean <- getFromNamespace("vecchia_mean", "GPvecchia")
+  V.singular <- FALSE
+  res <- try(vecchia.mean <- vecchia_mean(z, U.obj, V.ord), silent = TRUE)
+  if (inherits(res, "try-error")) {
+    warning("V is numerically singular. Predicted means are unreliable.")
+    V.ord.pd <- Matrix::nearPD(V.ord)$mat
+    V.singular <- TRUE
+    vecchia.mean <- vecchia_mean(z, U.obj, V.ord.pd)
+  }
+  return.list <- list(mu.pred = vecchia.mean$mu.pred,
+    mu.obs = vecchia.mean$mu.obs, var.pred = NULL, var.obs = NULL,
+    V.ord = NULL, U.obj = NULL)
+  if (return.values == "meanmat" || return.values == "all") {
+    return.list$V.ord <- V.ord
+    return.list$U.obj <- U.obj
+  }
+  if (return.values == "meanvar" || return.values == "all") {
+    if (V.singular) {
+      stop("V is numerically singular. Prediction variance cannot be computed.")
+    }
+    if (missing(var.exact))
+      var.exact <- (sum(!vecchia.approx$obs) < 4 * 10000)
+    vecchia_var <- getFromNamespace("vecchia_var", "GPvecchia")
+    vars.vecchia <- vecchia_var(U.obj, V.ord, exact = var.exact)
+    return.list$var.pred <- vars.vecchia$vars.pred
+    return.list$var.obs <- vars.vecchia$vars.obs
+  }
+  return.list
+}
+
 #' Multivariate Vecchia prediction
 #'
 #' This function is used to make predictions based on multivariate Vecchia
@@ -296,7 +334,14 @@ vecchia_Mprediction <- function(z, vecchia.approx, covparms, var.exact = NULL, r
   #    if (length(U.obj$zero.nugg) > 0)
   #        warning("Rows/cols of V have been removed for data with zero noise")
   vecchia_mean <- getFromNamespace("vecchia_mean", "GPvecchia")
-  vecchia.mean <- vecchia_mean(z, U.obj, V.ord)
+  V.singular <- FALSE
+  res <- try(vecchia.mean <- vecchia_mean(z, U.obj, V.ord), silent = TRUE)
+  if (inherits(res, "try-error")) {
+    warning("V is numerically singular. Predicted means are unreliable.")
+    V.ord.pd <- Matrix::nearPD(V.ord)$mat
+    V.singular <- TRUE
+    vecchia.mean <- vecchia_mean(z, U.obj, V.ord.pd)
+  }
   return.list <- list(
     mu.pred = vecchia.mean$mu.pred, mu.obs = vecchia.mean$mu.obs,
     var.pred = NULL, var.obs = NULL, V.ord = NULL, U.obj = NULL
@@ -306,6 +351,9 @@ vecchia_Mprediction <- function(z, vecchia.approx, covparms, var.exact = NULL, r
     return.list$U.obj <- U.obj
   }
   if (return.values == "meanvar" || return.values == "all") {
+    if (V.singular) {
+      stop("V is numerically singular. Prediction variance cannot be computed.")
+    }
     if (is.null(var.exact)) {
       var.exact <- (sum(!vecchia.approx$obs) < 4 * 10000)
     }
diff --git a/R/PrestoGP_Vecchia.R b/R/PrestoGP_Vecchia.R
index bf289d7..9f1a719 100644
--- a/R/PrestoGP_Vecchia.R
+++ b/R/PrestoGP_Vecchia.R
@@ -383,6 +383,9 @@ setMethod("impute_y_lod", "VecchiaModel", function(model, lodu, lodl, n.mi = 10,
 setMethod("specify", "VecchiaModel", function(model) {
   locs.scaled <- scale_locs(model, model@locs_train)
   model@vecchia_approx <- vecchia_specify(locs.scaled[[1]], model@n_neighbors)
+  if (model@omp_cores > -1) {
+    model@vecchia_approx$U.prep$n.cores <- model@omp_cores
+  }
   if (!model@common_scale) {
     olocs.scaled <- model@vecchia_approx$locsord
     for (j in 1:model@nscale) {
diff --git a/R/RcppExports.R b/R/RcppExports.R
index dd16aaf..03aeb6a 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -5,7 +5,7 @@ na_omit_c <- function(x) {
     .Call('_PrestoGP_na_omit_c', PACKAGE = 'PrestoGP', x)
 }
 
-createU_helper_mat <- function(olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning) {
-    .Call('_PrestoGP_createU_helper_mat', PACKAGE = 'PrestoGP', olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning)
+createU_helper_mat <- function(olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning, n_cores) {
+    .Call('_PrestoGP_createU_helper_mat', PACKAGE = 'PrestoGP', olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning, n_cores)
 }
 
diff --git a/man/PrestoGP-package.Rd b/man/PrestoGP-package.Rd
index 84e2edb..7b9c43a 100644
--- a/man/PrestoGP-package.Rd
+++ b/man/PrestoGP-package.Rd
@@ -23,6 +23,7 @@ Useful links:
 
 Authors:
 \itemize{
+  \item Eric Bair \email{eric.bair@sciome.com}
   \item Brian Kidd
   \item Eric Wimberley
   \item Deepak Mav
diff --git a/man/PrestoGPModel-class.Rd b/man/PrestoGPModel-class.Rd
index 70965db..48ca2dd 100644
--- a/man/PrestoGPModel-class.Rd
+++ b/man/PrestoGPModel-class.Rd
@@ -81,7 +81,11 @@ approximation. Ignored for full models.}
 
 \item{\code{common_scale}}{Do all columns of locs have the same scale parameter?}
 
-\item{\code{param_sequence}}{Records the indices of the various Matern parameters.
+\item{\code{param_sequence}}{Records the indices of the various Matern parameters.}
+
+\item{\code{omp_cores}}{Number of cores used by OMP to compute the U matrix.}
+
+\item{\code{logparams}}{Transformed version of the Matern parameters.
 See \code{\link{create_param_sequence}}.}
 }}
 
diff --git a/man/prestogp_fit-PrestoGPModel-method.Rd b/man/prestogp_fit-PrestoGPModel-method.Rd
index 8e3c2ea..71ac9cb 100644
--- a/man/prestogp_fit-PrestoGPModel-method.Rd
+++ b/man/prestogp_fit-PrestoGPModel-method.Rd
@@ -34,6 +34,7 @@
   family = c("gaussian", "binomial"),
   nfolds = 10,
   foldid = NULL,
+  omp.cores = -1,
   parallel = FALSE,
   cluster = NULL,
   adaptive = FALSE
@@ -145,13 +146,20 @@ what fold each observation should be assigned to in the cv.glmnet
 cross-validation. See \code{\link[glmnet]{cv.glmnet}} and
 \code{\link[ncvreg]{cv.ncvreg}}.}
 
+\item{omp.cores}{Number of parallel cores to be used by OMP when computing
+the U matrix in the Vecchia approximation. Defaults to -1 (which causes
+all available cores to be used).}
+
 \item{parallel}{Should parallel "foreach" be used to speed up the model
 fitting procedure where possible? Defaults to FALSE. Specifically,
 parallelization will be used for imputation and fitting the cv.glmnet
 object. See \code{\link[glmnet]{cv.glmnet}}. Note that this
 only applies to glmnet models (where penalty="lasso" or
 penalty="relaxed"). Models using ncvreg (where penalty="MCP" or
-penalty="SCAD") require a cluster argument for parallelization (see below).}
+penalty="SCAD") require a cluster argument for parallelization (see below).
+Note that parallelization will be used to compute the Vecchia
+approximation even if this parameter is set to FALSE (assuming that OMP is
+available on the system). To disable this, set omp.cores=1.}
 
 \item{cluster}{A cluster for running cv.ncvreg in parallel. See
 \code{\link[ncvreg]{cv.ncvreg}} and \code{\link[parallel]{makeCluster}}.
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 3f445f7..1f3f18b 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -23,8 +23,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // createU_helper_mat
-arma::sp_mat createU_helper_mat(const arma::mat& olocs, const arma::vec& ondx, const arma::mat& curqys, const arma::mat& curqzs, const arma::mat& vijs, const arma::mat& aijs, const arma::mat& full_const, const arma::vec& nugget, const arma::vec& sig2, const arma::vec& U_beginning);
-RcppExport SEXP _PrestoGP_createU_helper_mat(SEXP olocsSEXP, SEXP ondxSEXP, SEXP curqysSEXP, SEXP curqzsSEXP, SEXP vijsSEXP, SEXP aijsSEXP, SEXP full_constSEXP, SEXP nuggetSEXP, SEXP sig2SEXP, SEXP U_beginningSEXP) {
+arma::sp_mat createU_helper_mat(const arma::mat& olocs, const arma::vec& ondx, const arma::mat& curqys, const arma::mat& curqzs, const arma::mat& vijs, const arma::mat& aijs, const arma::mat& full_const, const arma::vec& nugget, const arma::vec& sig2, const arma::vec& U_beginning, const int n_cores);
+RcppExport SEXP _PrestoGP_createU_helper_mat(SEXP olocsSEXP, SEXP ondxSEXP, SEXP curqysSEXP, SEXP curqzsSEXP, SEXP vijsSEXP, SEXP aijsSEXP, SEXP full_constSEXP, SEXP nuggetSEXP, SEXP sig2SEXP, SEXP U_beginningSEXP, SEXP n_coresSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -38,14 +38,15 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< const arma::vec& >::type nugget(nuggetSEXP);
     Rcpp::traits::input_parameter< const arma::vec& >::type sig2(sig2SEXP);
     Rcpp::traits::input_parameter< const arma::vec& >::type U_beginning(U_beginningSEXP);
-    rcpp_result_gen = Rcpp::wrap(createU_helper_mat(olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning));
+    Rcpp::traits::input_parameter< const int >::type n_cores(n_coresSEXP);
+    rcpp_result_gen = Rcpp::wrap(createU_helper_mat(olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning, n_cores));
     return rcpp_result_gen;
 END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
     {"_PrestoGP_na_omit_c", (DL_FUNC) &_PrestoGP_na_omit_c, 1},
-    {"_PrestoGP_createU_helper_mat", (DL_FUNC) &_PrestoGP_createU_helper_mat, 10},
+    {"_PrestoGP_createU_helper_mat", (DL_FUNC) &_PrestoGP_createU_helper_mat, 11},
     {NULL, NULL, 0}
 };
 
diff --git a/src/createUMC.cpp b/src/createUMC.cpp
index 151acfd..d063ddb 100644
--- a/src/createUMC.cpp
+++ b/src/createUMC.cpp
@@ -51,7 +51,8 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx,
                                 const arma::mat &aijs,
                                 const arma::mat &full_const,
                                 const arma::vec &nugget, const arma::vec &sig2,
-                                const arma::vec &U_beginning) {
+                                const arma::vec &U_beginning,
+				const int n_cores) {
   int n = arma::as_scalar(ondx.n_elem);
   int m = arma::as_scalar(curqys.n_rows);
   int n_inds = 2 * n * (m + 3);
@@ -69,6 +70,9 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx,
   // feeder(2, arma::span(0,6)) = U_beginning.t();
   int ind = 7;
   // int ind = 0;
+  if (n_cores > -1) {
+    omp_set_num_threads(n_cores);
+  }
 #ifdef _OPENMP
 #pragma omp parallel for
 #endif
diff --git a/tests/testthat/test-PrestoGP_Model.R b/tests/testthat/test-PrestoGP_Model.R
index 28cb315..8b5ce2f 100644
--- a/tests/testthat/test-PrestoGP_Model.R
+++ b/tests/testthat/test-PrestoGP_Model.R
@@ -1,3 +1,48 @@
+test_that("omp.cores not numeric", {
+  source("sim_vecchia_small.R")
+  pgp.model1 <- new("VecchiaModel")
+  expect_error(
+    prestogp_fit(pgp.model1, y, X, locs, omp.cores = "foo"),
+    "omp.cores must be a positive integer >= 1"
+  )
+})
+
+test_that("omp.cores not a scalar", {
+  source("sim_vecchia_small.R")
+  pgp.model1 <- new("VecchiaModel")
+  expect_error(
+    prestogp_fit(pgp.model1, y, X, locs, omp.cores = c(1, 2)),
+    "omp.cores must be a positive integer >= 1"
+  )
+})
+
+test_that("omp.cores negative", {
+  source("sim_vecchia_small.R")
+  pgp.model1 <- new("VecchiaModel")
+  expect_error(
+    prestogp_fit(pgp.model1, y, X, locs, omp.cores = -2),
+    "omp.cores must be a positive integer >= 1"
+  )
+})
+
+test_that("omp.cores zero", {
+  source("sim_vecchia_small.R")
+  pgp.model1 <- new("VecchiaModel")
+  expect_error(
+    prestogp_fit(pgp.model1, y, X, locs, omp.cores = 0),
+    "omp.cores must be a positive integer >= 1"
+  )
+})
+
+test_that("omp.cores not an integer", {
+  source("sim_vecchia_small.R")
+  pgp.model1 <- new("VecchiaModel")
+  expect_error(
+    prestogp_fit(pgp.model1, y, X, locs, omp.cores = 2.5),
+    "omp.cores must be a positive integer >= 1"
+  )
+})
+
 test_that("beta.hat not numeric", {
   source("sim_vecchia_small.R")
   pgp.model1 <- new("VecchiaModel")