1 files changed, 29 insertions, 13 deletions
diff --git a/R/rmad.R b/R/rmad.R
index b8e9d41..46e9b35 100644
--- a/R/rmad.R
+++ b/R/rmad.R
@@ -1,31 +1,47 @@
-rmad <- function(x, y = NULL, na.rm = FALSE, even.correction = FALSE) {
+rmad <- function(x, y = NULL, na.rm = FALSE, even.correction = FALSE, num.threads = "half-max") {
+
+  ## check input data
   dat <- .check_input_data_matrix(x = x, y = y, na.rm = na.rm)
   colnames_original <- colnames(dat)
   storage.mode(dat) <- "double"
   n <- as.integer(nrow(dat))
   p <- as.integer(ncol(dat))
-  nc <- as.integer({
-    p^2 - p
-  }/2)
+
+
+  ## set even correction
   if (even.correction) {
     evencorrection <- 1L
-  }
-  else {
+  } else {
     evencorrection <- 0L
   }
-  u <- .Fortran("cormadvecdp", matrix = dat, nrow = n, ncol = p, res = numeric(nc), 
-    ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res
-  if (!is.null(y)) {
-    return(u)
+
+  ## set number of threads
+  if (num.threads == "half-max") {
+    num.threads <- 0L
+  } else {
+    storage.mode(num.threads) <- "integer"
   }
-  else {
+
+
+  ## Call C code
+  u <- .Call(C_cormad_C, dat, n, p, evencorrection, num.threads)
+
+
+  if (!is.null(y)) { ## 2-dimensional
+    return(u)
+  } else { ## p-dimensional
+
+    ## assemble the matrix using the lower triangle
     R <- Matrix(1, nrow = p, ncol = p, sparse = FALSE)
     R[lower.tri(R, diag = FALSE)] <- u
     R <- forceSymmetric(R, uplo = "L")
     R <- as(R, "dspMatrix")
+
+    ## attach dimnames if needed
     if (!is.null(colnames_original)) {
       dimnames(R)[[1]] <- dimnames(R)[[2]] <- colnames_original
     }
+
     return(R)
-  }
-}
+  } ## END if(!is.null(y)){ ## 2-dimensional
+} ## END function