version 2.0

18 files changed, 849 insertions, 871 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 3f611ee..dcafcfb 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -19,9 +19,8 @@ NeedsCompilation: yes
 Imports: stats, graphics, Matrix, methods, parallel, foreach,
         doParallel, utils
 License: GPL (>= 2)
-LazyData: TRUE
-Version: 1.2
-Date: 2020-07-24
-Packaged: 2020-07-24 07:47:23 UTC; pietro
+Version: 2.0
+Date: 2021-10-14
+Packaged: 2021-10-14 14:44:39 UTC; luco
 Repository: CRAN
-Date/Publication: 2020-07-24 14:00:07 UTC
+Date/Publication: 2021-10-17 20:00:08 UTC
diff --git a/MD5 b/MD5
index 7aaaa51..c43967d 100644
--- a/MD5
+++ b/MD5
@@ -1,17 +1,20 @@
-7a9c09c1b58001a6ef17502cccb6b9ee *DESCRIPTION
-1f64d30a870f2eae4855e83010aaaa4e *NAMESPACE
-1eae13396706a60fc95f5ff4b751765b *NEWS
-0b87a473118069181446fb007472c042 *R/check_inputs.R
-1df1ecba75175b1665069fc88419cf02 *R/cv_loss.R
-1ada6bc39f15ea58f34a5d46173141fa *R/plot_print_methods.R
-6dcfdaa8b26686b70794ba2c9d6147ca *R/rmad.R
-f26fdd51a2378e500e6ed2eaeda8721a *R/rsc.R
-556a9a10b85bfad178103f3cc5fdc0ff *R/rsc_cv.R
-607039efb30fc06f123290e3dcbce1d4 *R/zzz.R
+f8d0f0d74876f46c0818b142ffa180aa *DESCRIPTION
+9e6c9e63967c10a2dad469679534e3b2 *NAMESPACE
+12aa4cb77fcbd1096f7c868a332debac *NEWS
+852c9f1717e4dbd26ed8c549edc3c588 *R/check_inputs.R
+228fb1b443a4c9eaf4fddcdf9c8a31a9 *R/cv_loss.R
+8c8c70764c8a9ab72f23ba031e0f330d *R/plot_print_methods.R
+720ba075c0b402a232b5528774fa1230 *R/rmad.R
+3e98806cb071a0b84028b7e5bb022cc5 *R/rsc.R
+98992e4ed337da761f641cb2960f7579 *R/rsc_cv.R
+cafa1636fd3af0dd260cbfcb27df210e *R/zzz.R
 02a06053793e54ce413a378a05ccc3ce *inst/CITATION
 1e6af804927098cfa6c235278682578c *man/plot.cv_rsc.Rd
-57435eb3ed73e1af9435073eeaf57928 *man/rmad.Rd
+a762f06952b9b9f89cab9826674be146 *man/rmad.Rd
 8cc5c9373c7aa8835e9818d0de274ece *man/rsc.Rd
 63d3e58e9bdc0da9969feb6ce69d005f *man/rsc_cv.Rd
-f45a14eed865277afe99ebcb1bed66f8 *src/cormaddp.f90
-9ded1c66a2fd469e79d7d1d25d0bd604 *src/init.c
+95e3011e37d9dde0d75f3a3819b2acd3 *src/Makevars
+4aad47542a9a6ce7b9f153dccd1690aa *src/RSCdefines.h
+7ef7aff199844c0d42d07684e58e63d1 *src/cormad_ptr.c
+b518af5dd0a97f37ab37e4d915dc5303 *src/init.c
+cdb3d966f6db89bf54d3fac4b9314bc6 *src/selection_algos_ptr.c
diff --git a/NAMESPACE b/NAMESPACE
index ffd8478..3d5ba56 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -9,7 +9,7 @@ importFrom("doParallel", "registerDoParallel", "stopImplicitCluster")
 importFrom("utils", "citHeader", "citEntry")
 
 ## export objects in /src 
-useDynLib(RSC, .registration = TRUE, .fixes = "F_")
+useDynLib(RSC, .registration = TRUE, .fixes = "C_")
 
 ## export objects in /R
 export(rmad)
diff --git a/NEWS b/NEWS
index d7efdcf..f79bd3a 100755
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,13 @@
+## 14/10/2021 at 11:51:10 (CEST)

+* New version: v2.0.

+   + Implemented the low-level routines in C.

+   + Moved from introselect to quickselect (median calculation);

+      faster on benchmarks with intended input.

+   + Added support for parallel execution.

+

+## 07/10/2021 at 11:51:10 (CEST)

+* v1.3 does not compile on Solaris (syntax issues)

+

 ## 24/07/2020 at 09:44:39 (CEST)

 * v1.2 uploaded

 * fixed compiler issue on Solaris

diff --git a/R/check_inputs.R b/R/check_inputs.R
index 4cbcf0a..ccc54d4 100644
--- a/R/check_inputs.R
+++ b/R/check_inputs.R
@@ -1,64 +1,81 @@
-.check_input_data_matrix <- function(x, y, na.rm) {
-  if (is.null(y)) {
-    if (!is.array(x) & !is.data.frame(x) & !is.matrix(x)) {
-      stop("\"x\" must be a numeric matrix or any other array-type object that can be converted or a \"matrix\" object with with ncol>=2.\n\n")
-    }
-    if (is.vector(x)) {
-      stop("\"x\" must be a numeric matrix or any other array-type object that can be converted or a \"matrix\" object.  with with ncol>=2.\n\n")
-    }
-    if (!is.matrix(x)) {
-      x <- data.matrix(x)
-    }
-    if (!is.numeric(x)) {
-      stop("\"x\" must be numeric.")
-    }
-    if (nrow(x) < 2 | ncol(x) < 2) {
-      stop("nrow(x)>=2 and ncol(xa)>=2 are required\n\n")
-    }
-    is_na_data <- is.na(x)
-    if (any(is_na_data)) {
-      if (na.rm == FALSE) {
-        stop("\"x\" contains NA records. You may want to filter NAs by setting \"na.rm=TRUE\" (see documentation for more details).\n\n")
-      }
-      else {
-        idx_na <- which(rowSums(is_na_data) >= 1)
-        x <- x[-idx_na, , drop = FALSE]
-        if (nrow(x) < 2) {
-          stop("nrow(x)<2 after NA removal.\n\n")
+## Check input data and return a valid matrix object
+##
+.check_input_data_matrix <- function(x , y , na.rm){
+
+    if(is.null(y)){
+        if(!is.array(x) &  !is.data.frame(x) & !is.matrix(x)){
+            stop('"x" must be a numeric matrix or any other array-type object that can be converted or a "matrix" object with with ncol>=2.\n\n')  
         }
-      }
-    }
-    if (any(!is.finite(x))) {
-      stop("\"x\" contains Inf values\n\n")
-    }
-  }
-  else {
-    if (!is.vector(x) | !is.vector(y)) {
-      stop("If \"y\" is given, \"x\" and  \"y\" must be both numeric.\n\n")
-    }
-    if (!is.numeric(x) | !is.numeric(y)) {
-      stop("\"x\" and \"y\" must be numeric.\n\n")
-    }
-    if (length(x) != length(y)) {
-      stop("\"x\" and \"y\" have different length.\n\n")
-    }
-    x <- cbind(x, y, deparse.level = 0)
-    is_na_data <- is.na(x)
-    if (any(is_na_data)) {
-      if (na.rm == FALSE) {
-        stop("\"x\" or \"y\" contains NA records. You may want to filter NAs by setting \"na.rm=TRUE\" (see documentation for more details).\n\n")
-      }
-      else {
-        idx_na <- which(rowSums(is_na_data) >= 1)
-        x <- x[-idx_na, , drop = FALSE]
-        if (nrow(x) < 2) {
-          stop("length(x)<2 and/or length(y)<2 after NA removal.\n\n")
+
+        if(is.vector(x)){
+            stop('"x" must be a numeric matrix or any other array-type object that can be converted or a "matrix" object.  with with ncol>=2.\n\n')  
         }
-      }
-    }
-    if (any(!is.finite(x))) {
-      stop("\"x\" and/or \"y\"  contains Inf values\n\n")
+
+        if(!is.matrix(x)){
+            x <- data.matrix(x)
+        }
+        
+        if( !is.numeric(x) ){
+            stop('"x" must be numeric.')
+        }
+        
+        if(nrow(x)<2 | ncol(x)<2){
+            stop('nrow(x)>=2 and ncol(xa)>=2 are required\n\n')  
+        }
+
+        is_na_data <- is.na(x)
+        if(any(is_na_data)){
+            if(na.rm==FALSE){
+                stop('"x" contains NA records. You may want to filter NAs by setting "na.rm=TRUE" (see documentation for more details).\n\n')
+            }else{
+                idx_na <- which( rowSums(is_na_data) >=1 )
+                x   <- x[-idx_na , , drop = FALSE]
+                if(nrow(x)<2){
+                    stop('nrow(x)<2 after NA removal.\n\n')  
+                }   
+            }
+        }
+
+        if(any(!is.finite(x))){
+            stop('"x" contains Inf values\n\n')  
+        }
+        
+    }else{
+
+        if(!is.vector(x) | !is.vector(y)){
+            stop('If "y" is given, "x" and  "y" must be both numeric.\n\n')
+        }
+        
+        if(!is.numeric(x) | !is.numeric(y)){
+            stop('"x" and "y" must be numeric.\n\n')
+        }
+
+        if(length(x) != length(y)){
+            stop('"x" and "y" have different length.\n\n')
+        }
+
+        x <- cbind( x , y , deparse.level = 0) 
+        
+        is_na_data <- is.na(x)
+        if(any(is_na_data)){
+            if(na.rm==FALSE){
+                stop('"x" or "y" contains NA records. You may want to filter NAs by setting "na.rm=TRUE" (see documentation for more details).\n\n')
+            }else{
+                idx_na <- which( rowSums(is_na_data) >=1 )
+                x   <- x[-idx_na , , drop = FALSE]
+                if(nrow(x)<2){
+                    stop('length(x)<2 and/or length(y)<2 after NA removal.\n\n')  
+                }   
+            }
+        }
+
+        if(any(!is.finite(x))){
+            stop('"x" and/or "y"  contains Inf values\n\n')  
+        }   
     }
-  }
-  return(x)
-}
+
+    
+    return(x)
+    
+} ### END function 
+
diff --git a/R/cv_loss.R b/R/cv_loss.R
index ea46007..d855362 100644
--- a/R/cv_loss.R
+++ b/R/cv_loss.R
@@ -1,17 +1,29 @@
-.cv_loss <- function(idx, dat, evencorrection, threshold, grid.length, p, nc) {
-  res <- numeric(nc)
+## Normalized squared Frobenius loss for all threshold values at a given train/test
+## set, where
+##
+##   * idx = TRUE for sample points into the train set
+##   * train set = dat[  idx , ]
+##   * test  set = dat[ !idx , ]
+##
+.cv_loss <- function(idx, dat, evencorrection, threshold, grid.length, p) {
+
+
+  ## compute RMAD  on train set
   n1 <- as.integer(sum(idx))
-  C1 <- .Fortran("cormadvecdp", matrix = dat[idx, ], nrow = n1, ncol = p, res = res, 
-    ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res
+  C1 <- .Call(C_cormad_C, dat[idx, ], n1, p, evencorrection, num.threads = 1)
+
+  ## compute RMAD  on test set
   n2 <- as.integer(sum(!idx))
-  C2 <- .Fortran("cormadvecdp", matrix = dat[!idx, ], nrow = n2, ncol = p, res = res, 
-    ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res
+  C2 <- .Call(C_cormad_C, dat[!idx, ], n2, p, correcteven = evencorrection, num.threads = 1)
+
+  ## apply thresholds
   ans <- rep(0, times = grid.length)
   for (h in 1:grid.length) {
-    C1[abs(C1) < threshold[h]] <- 0
+    C1[abs(C1) < threshold[h]] <- 0 ## fit on train set
     ans[h] <- sum(2 * {
       C1 - C2
-    }^2)/p
+    }^2) / p
   }
+
   return(ans)
 }
diff --git a/R/plot_print_methods.R b/R/plot_print_methods.R
index ca163cb..746179a 100644
--- a/R/plot_print_methods.R
+++ b/R/plot_print_methods.R
@@ -1,26 +1,41 @@
-print.rsc_cv <- function(x, ...) {
-  cat("\n")
-  cat("================================================\n")
-  cat("   Expected Normalized Squared Frobenius Loss   \n")
-  cat("================================================\n")
-  print(x$loss)
-  cat("================================================\n")
-  cat("\n")
+
+## print method for the crsc_cv class
+print.rsc_cv <-function(x, ...){
+    cat("\n")
+    cat("================================================\n")
+    cat("   Expected Normalized Squared Frobenius Loss   \n")
+    cat("================================================\n")
+    print(x$loss)
+    cat("================================================\n")
+    cat("\n")
 }
-plot.rsc_cv <- function(x, ...) {
-  tstar <- which(x$loss$Flag == "minimum")
-  hstar <- x$loss$Threshold[tstar]
-  inf_loss <- x$loss$Average - x$loss$SE
-  sup_loss <- x$loss$Average + x$loss$SE
-  a <- inf_loss[tstar]
-  b <- sup_loss[tstar]
-  hstar1se <- max(x$loss$Threshold[which(x$loss$Flag == "*")])
-  Ylim <- range(c(inf_loss, sup_loss))
-  plot(x$loss$Threshold, x$loss$Average, t = "b", ylim = Ylim, pch = 20, col = "#0052A5", 
-    lwd = 2, main = "RSC Optimal Threshold Selection", xlab = "Threshold", ylab = "Average loss", 
-    ...)
-  arrows(x$loss$Threshold, inf_loss, x$loss$Threshold, sup_loss, length = 0.05, 
-    angle = 90, code = 3, col = "#0052A5")
-  abline(v = hstar1se, col = "#31A853", lty = 2, lwd = 2)
-  abline(v = hstar, col = "#E0162B", lty = 2, lwd = 2)
+
+
+
+## ## ## plot method for the crsc_cv class
+plot.rsc_cv <- function(x, ...){
+
+    ## add check object
+    
+    tstar <- which(x$loss$Flag == "minimum")
+    hstar <- x$loss$Threshold[tstar]
+    
+    inf_loss <- x$loss$Average - x$loss$SE
+    sup_loss <- x$loss$Average + x$loss$SE
+    a <- inf_loss[tstar]
+    b <- sup_loss[tstar]
+    hstar1se <- max(x$loss$Threshold[which(x$loss$Flag == "*")])
+
+    Ylim     <- range(c(inf_loss, sup_loss))
+    plot(x$loss$Threshold, x$loss$Average, t='b', ylim = Ylim ,
+         pch=20 , col= "#0052A5", lwd = 2,
+         main = "RSC Optimal Threshold Selection",
+         xlab = "Threshold",
+         ylab = "Average loss", ...)
+    arrows(x$loss$Threshold, inf_loss, x$loss$Threshold, sup_loss ,
+           length=.05, angle=90, code=3, col="#0052A5")
+    abline(v = hstar1se,     col = "#31A853", lty=2, lwd=2)
+    abline(v = hstar,        col = "#E0162B", lty=2, lwd=2)
 }
+
+
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
diff --git a/R/rsc.R b/R/rsc.R
index 6cf4e15..8182d1c 100644
--- a/R/rsc.R
+++ b/R/rsc.R
@@ -1,42 +1,49 @@
-rsc <- function(cv, threshold = "minimum") {
-  if (class(cv) == "rsc_cv") {
-    if (is.numeric(threshold)) {
-      if (length(threshold) > 1) {
-        stop("if a specific value for 'threshold' is chosen, this must be a single numeric value in (0,1)")
-      }
-      else if (threshold <= 0 | threshold >= 1) {
-        stop("if a specific value for 'threshold' is chosen, this must be a single numeric value in (0,1)")
-      }
-    }
-    else {
-      if ({
-        threshold != "minimum"
-      } & {
-        threshold != "minimum1se"
-      }) {
-        stop("'threshold' must be one of the following: 'minimum', 'minimum1se', a numeric value in (0,1).")
-      }
-      if (threshold == "minimum") {
-        threshold <- cv$minimum
-      }
-      else if (threshold == "minimum1se") {
-        threshold <- cv$minimum1se
-      }
-    }
-    cv$rmadvec[abs(cv$rmadvec) < threshold] <- 0
-    nc <- length(cv$rmadvec)
-    p <- {
-      1 + sqrt(1 + 8 * nc)
-    }/2
-    R <- Matrix(1, nrow = p, ncol = p, sparse = TRUE)
-    R[lower.tri(R, diag = FALSE)] <- cv$rmadvec
-    R <- forceSymmetric(R, uplo = "L")
-    if (!is.null(cv$varnames)) {
-      dimnames(R)[[1]] <- dimnames(R)[[2]] <- cv$varnames
-    }
-  }
-  else {
-    stop("'cv' must be a an object of class 'rsc_cv' obtained from 'rsc::rsc_cv'")
-  }
-  return(R)
+rsc <- function(cv, threshold = "minimum"){
+
+    ## inputs
+    ## cv         = u   ## a class cv_rsc or any other correlation matrix
+    ## threshold  = "minimum" ## "minimum", "minimum1se" or numeric in (0,1)
+
+    if(class(cv) == "rsc_cv"){
+
+        ## check threshold
+        if(is.numeric(threshold)){
+            if(length(threshold)>1){
+                stop("if a specific value for 'threshold' is chosen, this must be a single numeric value in (0,1)")
+            }else if(threshold <=0 | threshold >=1){
+                stop("if a specific value for 'threshold' is chosen, this must be a single numeric value in (0,1)")
+            }
+        }else{
+            if({threshold != "minimum"} & {threshold != "minimum1se"}){
+                stop("'threshold' must be one of the following: 'minimum', 'minimum1se', a numeric value in (0,1).")
+            }
+
+            if(threshold == "minimum"){
+                threshold <- cv$minimum
+            }else if(threshold == "minimum1se"){
+                threshold <- cv$minimum1se
+            }
+        }
+
+
+        ## threshold the rmadvec 
+        cv$rmadvec[ abs(cv$rmadvec) < threshold ] <- 0
+
+        nc <- length(cv$rmadvec)
+        p  <- {1 + sqrt( 1 + 8 * nc ) } / 2
+        R  <- Matrix(1, nrow = p, ncol = p, sparse = TRUE)
+
+        R[lower.tri(R , diag = FALSE)]  <- cv$rmadvec
+        R <- forceSymmetric(R , uplo="L")
+
+        ## attach dimnames if needed 
+        if(!is.null(cv$varnames)){
+            dimnames(R)[[1]] <- dimnames(R)[[2]] <- cv$varnames
+        }
+    }else{
+        stop("'cv' must be a an object of class 'rsc_cv' obtained from 'rsc::rsc_cv'")
+        }
+
+
+    return(R)
 }
diff --git a/R/rsc_cv.R b/R/rsc_cv.R
index 4153322..660695e 100644
--- a/R/rsc_cv.R
+++ b/R/rsc_cv.R
@@ -1,79 +1,140 @@
-rsc_cv <- function(x, cv.type = "kfold", R = 10, K = 10, threshold = seq(0.05, 0.95, 
-  by = 0.025), even.correction = FALSE, na.rm = FALSE, ncores = NULL, monitor = TRUE) {
+rsc_cv <- function(x,
+                   cv.type = "kfold",
+                   R = 10,
+                   K = 10,
+                   threshold = seq(0.05, 0.95, by = 0.025),
+                   even.correction = FALSE,
+                   na.rm = FALSE,
+                   ncores = NULL,
+                   monitor = TRUE) {
+
+  ## ## *****************************************************************************
+  ## ## RSC inputs
+  ## ## *****************************************************************************
+  ## x               =  X
+  ## cv.type         =  "random" ## "kfold"   "random"
+  ## R               =  10      ## replicate (for K-fold) or splits for random
+  ## K               =  10      ## folds in kfcv
+  ## threshold       =  seq(0.025, 0.975, by = 0.025)
+  ## opt             = "min" , ## "min" "min1se"
+  ## even.correction =  FALSE
+  ## na.rm           =  TRUE
+  ## ncores          =  6
+  ## monitor         =  TRUE
+  ## ## *****************************************************************************
+
+
+  ## check input data
   dat <- .check_input_data_matrix(x = x, y = NULL, 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)
+
+  ## check cv.type
   if ({
     cv.type != "random"
   } & {
     cv.type != "kfold"
   }) {
-    stop("\"cv.type\" must be either \"random\" (default)  or \"kfold\"")
+    stop('"cv.type" must be either "random" (default)  or "kfold"')
   }
+
+
+  ## check R
   if (!is.numeric(R)) {
-    stop("\"R\" must be an integer > 1")
-  }
-  else if (R < 1) {
-    stop("\"R\" must be an integer > 1")
+    stop('"R" must be an integer > 1')
+  } else if (R < 1) {
+    stop('"R" must be an integer > 1')
   }
+
+
+  ## check K
   if (!is.numeric(K)) {
-    stop("\"K\" must be an integer > 1")
-  }
-  else if (R < 1) {
-    stop("\"K\" must be an integer > 1")
+    stop('"K" must be an integer > 1')
+  } else if (R < 1) {
+    stop('"K" must be an integer > 1')
   }
+
+
+  ## check threshold
   if (length(threshold) == 1) {
     if (threshold <= 0 | threshold >= 1) {
-      stop("\"threshold\" value does not belong to the interval (0,1).")
+      stop('"threshold" value does not belong to the interval (0,1).')
     }
-  }
-  else if (length(threshold) > 1) {
+  } else if (length(threshold) > 1) {
+    ##
     if (any(threshold < 0) | any(threshold >= 1)) {
-      stop("Some of the \"threshold\" values do not belong to the interval (0,1).")
+      stop('Some of the "threshold" values do not belong to the interval (0,1).')
     }
     grid.length <- length(threshold)
   }
+
+
+
+
+
+  ## set even correction
   if (even.correction) {
     evencorrection <- 1L
-  }
-  else {
+  } else {
     evencorrection <- 0L
   }
+
+
+
+  ## set and check ncores
   if (is.null(ncores)) {
     DetectedCores <- detectCores()
     if (DetectedCores <= 2) {
       ncores <- 1
-    }
-    else {
+    } else {
       ncores <- {
         DetectedCores - 1
       }
     }
-  }
-  else {
+  } else {
     ncores <- as.integer(ncores)
     if (ncores <= 0) {
-      stop("\"ncores\" must be an integer larger or equal to 1.")
+      stop('"ncores" must be an integer larger or equal to 1.')
     }
   }
+
+
+
+
+
+
+
+
+
+
+  ## compute the RMAD vec form
   if (monitor) {
     cat("\n")
     message("Computing the RMAD matrix")
     t0 <- Sys.time()
   }
-  rmad_vec <- .Fortran("cormadvecdp", matrix = dat, nrow = n, ncol = p, res = numeric(nc), 
-    ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res
+
+  rmad_vec <- .Call(C_cormad_C, dat, n, p, evencorrection, num.threads = 1)
+
   if (monitor) {
     t1 <- Sys.time()
     dt01 <- difftime(t1, t0, units = "auto")
-    message("* RMAD computing time:...... ", round(dt01, 2), " [", attributes(dt01)$units, 
-      "]")
+    message(
+      "* RMAD computing time:...... ", round(dt01, 2),
+      " [", attributes(dt01)$units, "]"
+    )
   }
+
+
+
+
+
+
+
+
+
   if (cv.type == "random") {
     if (monitor) {
       cat("\n")
@@ -81,28 +142,55 @@ rsc_cv <- function(x, cv.type = "kfold", R = 10, K = 10, threshold = seq(0.05, 0
       t_hat <- round(1.2 * {
         {
           dt01 * R * 2
-        }/ncores
+        } / ncores
       }, 2)
       message("* predicted end time (worst case):...... ", Sys.time() + t_hat)
     }
-    n1 <- n - floor(n/log(n))
+
+
+    ## IDX[k,i] = TRUE means dat[i, ] is in train set at the k-th split
+    n1 <- n - floor(n / log(n)) ## train set
     IDX <- array(FALSE, dim = c(R, n))
     for (r in 1:R) {
       IDX[r, ][sample(1:n, size = n1, replace = FALSE)] <- TRUE
     }
+
+    ## register parallel backend
     registerDoParallel(ncores)
+
+    ## parallel computation of losses over splits
     U <- foreach(r = 1:R) %dopar% {
-      .cv_loss(idx = IDX[r, ], dat = dat, evencorrection = evencorrection, 
-        threshold = threshold, grid.length = grid.length, p = p, nc = nc)
+      .cv_loss(
+        idx = IDX[r, ], dat = dat, evencorrection = evencorrection,
+        threshold = threshold, grid.length = grid.length,
+        p = p
+      )
     }
+
+    ## stop parallel backend
     stopImplicitCluster()
+
+    ## LOSS[split , threshold]
+    ## array with normalized squared Frobenius loss
     LOSS <- array(0, dim = c(R, grid.length))
     for (r in 1:R) {
       LOSS[r, ] <- U[[r]]
     }
+
+    ## Estimate average loss with std errors
     avg_loss <- apply(LOSS, 2, mean)
-    se_loss <- apply(LOSS, 2, sd)/sqrt(R)
+    se_loss <- apply(LOSS, 2, sd) / sqrt(R)
   }
+
+
+
+
+
+
+
+
+
+
   if (cv.type == "kfold") {
     if (monitor) {
       cat("\n")
@@ -110,35 +198,59 @@ rsc_cv <- function(x, cv.type = "kfold", R = 10, K = 10, threshold = seq(0.05, 0
       t_hat <- round(1.2 * {
         {
           dt01 * R * K * 2
-        }/ncores
+        } / ncores
       }, 2)
       message("* predicted end time (worst case):...... ", Sys.time() + t_hat)
     }
+
+    ## create K deterministic folds
     idx_fold <- cut(1:n, breaks = K, labels = FALSE)
+
+    ## IDX :: indexes of all train sets
+    ## IDX[k , i] = TRUE means dat[i, ] is in the train set at some k-th  split
+    ## rows {{r-1}*K + 1 }:{r*K} of IDX correspond to the r-th replica
     IDX <- array(TRUE, dim = c(R * K, n))
+
+    ## set initial IDX row counter
     row_count <- 1L
     for (r in 1:R) {
+      ## at each replicate r shuffle the original fold indexes
       idx_fold_shuffle <- sample(idx_fold, size = n, replace = FALSE)
+
+      ## for each fold shuffle make up the K-fold
       for (k in 1:K) {
         IDX[row_count, ][idx_fold_shuffle == k] <- FALSE
         row_count <- 1L + row_count
       }
     }
+
+
+    ## register parallel backend
     registerDoParallel(ncores)
+
+    ## parallel computation of losses over splits
     U <- foreach(r = 1:{
       R * K
     }) %dopar% {
-      .cv_loss(idx = IDX[r, ], dat = dat, evencorrection = evencorrection, 
-        threshold = threshold, grid.length = grid.length, p = p, nc = nc)
+      .cv_loss(
+        idx = IDX[r, ], dat = dat, evencorrection = evencorrection,
+        threshold = threshold, grid.length = grid.length,
+        p = p
+      )
     }
+
+    ## stop parallel backend
     stopImplicitCluster()
+
+
+    ## LOSS[fold , threshold, replica]
     if (R == 1) {
       LOSS <- array(0, dim = c(K, grid.length))
+      ## names ??
       for (k in 1:K) {
         LOSS[k, ] <- U[[k]]
       }
-    }
-    else {
+    } else {
       LOSS <- array(0, dim = c(K, grid.length, R))
       dimnames(LOSS)[[3]] <- paste0("r", 1:R)
       dimnames(LOSS)[[1]] <- paste0("k", 1:K)
@@ -150,21 +262,33 @@ rsc_cv <- function(x, cv.type = "kfold", R = 10, K = 10, threshold = seq(0.05, 0
         }
       }
     }
+
+
+    ## compute average losses and standard errors
     avg_loss_r <- sd_loss_r <- matrix(0, nrow = R, ncol = grid.length)
     for (r in 1:R) {
       avg_loss_r[r, ] <- apply(LOSS[, , r], 2, mean)
-      sd_loss_r[r, ] <- apply(LOSS[, , r], 2, sd)/sqrt(K)
+      sd_loss_r[r, ] <- apply(LOSS[, , r], 2, sd) / sqrt(K)
     }
     avg_loss <- apply(avg_loss_r, 2, mean)
     se_loss <- apply(sd_loss_r, 2, mean)
+
+
+
     if (monitor) {
       t2 <- Sys.time()
       dt02 <- difftime(t2, t0, units = "auto")
     }
-  }
+  } ## END if(cv.type == "kfold"){
+
+
+
   if (monitor) {
     message("* finished on:.......................... ", Sys.time())
   }
+
+
+  ## Organize cross-validation results
   tstar <- which.min(avg_loss)
   flags <- rep("", grid.length)
   a <- avg_loss[tstar] - se_loss[tstar]
@@ -175,9 +299,23 @@ rsc_cv <- function(x, cv.type = "kfold", R = 10, K = 10, threshold = seq(0.05, 0
     avg_loss <= b
   }] <- "*"
   flags[tstar] <- "minimum"
-  res <- data.frame(Threshold = threshold, Average = avg_loss, SE = se_loss, Flag = flags)
-  ans <- list(rmadvec = rmad_vec, varnames = colnames_original, loss = res, minimum = threshold[tstar], 
-    minimum1se = max(threshold[avg_loss >= a & avg_loss <= b]))
+  res <- data.frame(
+    Threshold = threshold,
+    Average = avg_loss,
+    SE = se_loss,
+    Flag = flags
+  )
+
+  ## output list
+  ans <- list(
+    rmadvec = rmad_vec,
+    varnames = colnames_original,
+    loss = res,
+    minimum = threshold[tstar],
+    minimum1se = max(threshold[avg_loss >= a & avg_loss <= b])
+  )
+
+
   class(ans) <- "rsc_cv"
   return(ans)
 }
diff --git a/R/zzz.R b/R/zzz.R
index 3c33499..ca8c78e 100755
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -1,4 +1,7 @@
-.onAttach <- function(lib, pkg) {
-  packageStartupMessage("\nRSC: robust and sparse correlation matrix estimation\n     Type 'citation(\"RSC\")' for citing this package\n")
-  invisible()
+.onAttach <- function(lib, pkg){
+   packageStartupMessage("\nRSC: robust and sparse correlation matrix estimation\n     Type 'citation(\"RSC\")' for citing this package\n")
+   invisible()
 }
+
+
+
diff --git a/man/rmad.Rd b/man/rmad.Rd
index 01247cf..26417a0 100755
--- a/man/rmad.Rd
+++ b/man/rmad.Rd
@@ -12,7 +12,7 @@
 
 
 \usage{
-  rmad(x , y = NULL, na.rm = FALSE , even.correction = FALSE)   
+  rmad(x , y = NULL, na.rm = FALSE , even.correction = FALSE, num.threads = "half-max")   
 }
 
 
@@ -39,6 +39,9 @@
     for the calculation of the medians is applied to reduce the bias
     when the number of samples even (see \emph{Details}).
   }
+  \item{num.threads}{
+     An integer value or the string \code{"half-max"} (default), specifying the number of threads for parallel execution (see \emph{Details}).
+  }
 }
 
 
@@ -74,6 +77,12 @@
    package. The latter is  specifically designed for dense real symmetric
    matrices. A sparse correlation matrix can be obtained applying
    thresholding using the \code{\link{rsc_cv}} and \code{\link{rsc}}.
+
+   \code{rmad} function supports parallel execution.
+   This is provided via \emph{openmp} (http://www.openmp.org), which must be available; otherwise, falls back to single-core execution.
+   If \code{num.threads > 0}, function is executed using \code{min(num.threads, max.threads)} threads, where \code{max.threads} is the maximum number of available threads. That is, if positive use the specified number of threads (up to the maximum available).
+   If \code{num.threads < 0}, function is executed using \code{max(max.threads - num.threads, 1)} threads, i.e. when negative \code{num.threads} indicates the number of threads not to use.
+   If \code{num.threads == 0} or \code{num.threads == "half-max"}, function is executed using half of the available threads (\code{max(max.threads/2, 1)}). This is the default.
 }
 
 
diff --git a/src/Makevars b/src/Makevars
new file mode 100644
index 0000000..0ee3615
--- /dev/null
+++ b/src/Makevars
@@ -0,0 +1,2 @@
+PKG_CFLAGS = $(SHLIB_OPENMP_CFLAGS)
+PKG_LIBS = $(SHLIB_OPENMP_CFLAGS)
diff --git a/src/RSCdefines.h b/src/RSCdefines.h
new file mode 100644
index 0000000..75203ac
--- /dev/null
+++ b/src/RSCdefines.h
@@ -0,0 +1,11 @@
+/* Defines quickselect and cormad (double precision) */
+
+double quickselect_recursive(double *vector_extract_k, int vec_size, int extract_this_element);
+
+void cormad(double *matrix, int n_row, int n_col, double *output, int evencorrect);
+
+#ifdef _OPENMP
+void cormad_parallel(double *matrix, int n_row, int n_col, double *output,
+                      int evencorrect, int num_threads);
+#endif
+
diff --git a/src/cormad_ptr.c b/src/cormad_ptr.c
new file mode 100644
index 0000000..670a71a
--- /dev/null
+++ b/src/cormad_ptr.c
@@ -0,0 +1,180 @@
+#include "RSCdefines.h"
+#include <math.h>
+#include <stdio.h>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#define SQRT2 (sqrt(2))
+#define CONST 1.4826
+
+void cormad(double *matrix, int n_row, int n_col, double *output,
+            int evencorrect) {
+
+  int k = n_row / 2; // position of the median
+  int output_size = (n_col - 1) * n_col / 2;
+  double med, mad;           // store variables for medians
+  double U[n_row], V[n_row]; // help vectors
+
+  /* Transform matrix columns (CORMAD part 1)
+   * each column has n entries;
+   * matrix is assumed to be streamed in onevector
+   */
+  int i = 0;
+  for (int l = 0; l < n_col * n_row; l++) {
+    if (i == (n_row - 1)) {
+      U[i] = *matrix;
+      med = quickselect_recursive(U, n_row, k);
+      if (evencorrect == 1) { // handle even correction
+        med = (med + quickselect_recursive(U, n_row, k - 1)) / 2;
+      }
+      for (int j = 0; j < n_row; j++) {
+        U[j] = *(matrix - (n_row - 1) + j) - med;
+        V[j] = fabs(U[j]);
+      }
+      med = quickselect_recursive(V, n_row, k);
+      if (evencorrect == 1) { // handle even correction
+        med = (med + quickselect_recursive(V, n_row, k - 1)) / 2;
+      }
+      for (int j = 0; j < n_row; j++) // reassign
+        *(matrix - (n_row - 1) + j) = U[j] / (SQRT2 * CONST * med);
+      // prepare next iter
+      matrix++;
+      i = 0;
+    } else {
+      U[i] = *matrix;
+      // prepare next iter
+      matrix++;
+      i++;
+    }
+  }
+
+  matrix = matrix - (n_col * n_row); /* reset pointer */
+
+  /* Operate on columns pairs (CORMAD part 2) */
+  double *matrix_2 = matrix; // used to point at second column
+  int first_col = 0;         /* Running first column */
+  int second_col = 0;        /* Running second column */
+  for (int l = 0; l < output_size; l++) {
+    if (second_col == n_col - 1) {
+      first_col++;
+      second_col = first_col;
+      /* set pointers to columns */
+      matrix += n_row;
+      matrix_2 = matrix;
+    }
+    second_col++;
+    matrix_2 += n_row;
+
+    for (int i = 0; i < n_row; i++) { // auxiliary vectors from matrix
+      U[i] = *(matrix + i) + *(matrix_2 + i);
+      V[i] = -*(matrix + i) + *(matrix_2 + i);
+    }
+    mad = quickselect_recursive(U, n_row, k);
+    med = quickselect_recursive(V, n_row, k);
+    if (evencorrect == 1) {
+      mad = (mad + quickselect_recursive(U, n_row, k - 1)) / 2;
+      med = (med + quickselect_recursive(V, n_row, k - 1)) / 2;
+    }
+    for (int i = 0; i < n_row; i++) { // reassign for new medians
+      U[i] = fabs(U[i] - mad);
+      V[i] = fabs(V[i] - med);
+    }
+    mad = quickselect_recursive(U, n_row, k);
+    med = quickselect_recursive(V, n_row, k);
+    if (evencorrect == 1) {
+      mad = (mad + quickselect_recursive(U, n_row, k - 1)) / 2;
+      med = (med + quickselect_recursive(V, n_row, k - 1)) / 2;
+    }
+    mad = pow(CONST * mad, 2);
+    med = pow(CONST * med, 2);
+
+    // Assign output
+    *output = (mad - med) / (mad + med);
+    output++;
+  }
+}
+
+#ifdef _OPENMP
+void cormad_parallel(double *matrix, int n_row, int n_col, double *output,
+                     int evencorrect, int num_threads) {
+  int k = n_row / 2; // position of the median
+  int output_size = (n_col - 1) * n_col / 2;
+  double med, mad;           // store variables for medians
+  double U[n_row], V[n_row]; // help vectors
+
+  /* Transform matrix columns (CORMAD part 1) */
+  double *help_matrix = matrix; /*  help pointer for matrix */
+#pragma omp parallel for private(med, mad, U, V, help_matrix)                  \
+    num_threads(num_threads)
+  for (int j = 0; j < n_col; j++) {   // iterate on cols
+    help_matrix = matrix + n_row * j; /* set pointer at beg of col */
+    for (int i = 0; i < n_row; i++)
+      U[i] = *(help_matrix + i);
+    med = quickselect_recursive(U, n_row, k);
+    if (evencorrect == 1) { // handle even correction
+      med = (med + quickselect_recursive(U, n_row, k - 1)) / 2;
+    }
+    for (int i = 0; i < n_row; i++) {
+      U[i] = *(help_matrix + i) - med;
+      V[i] = fabs(U[i]);
+    }
+    med = quickselect_recursive(V, n_row, k);
+    if (evencorrect == 1) { // handle even correction
+      med = (med + quickselect_recursive(V, n_row, k - 1)) / 2;
+    }
+    for (int i = 0; i < n_row; i++) // reassign
+      *(help_matrix + i) = U[i] / (SQRT2 * CONST * med);
+  }
+
+  // int l = 0; // used to iterate over output
+  /* Operate on columns pairs (CORMAD part 2) */
+  help_matrix = matrix; /* help pointers for matrix */
+  double *help_matrix_2 = matrix;
+#pragma omp parallel for num_threads(num_threads) private(                     \
+    med, mad, U, V, help_matrix, help_matrix_2)
+  for (int l = 0; l < output_size; l++) {
+    /*  Detrmine columns pairs */
+    int col1 = 0, col2 = 0;
+    int copy_l = l + 1;
+    for (int last_elem = n_col - 1; last_elem > 0; last_elem--) {
+      copy_l -= last_elem;
+      if (copy_l <= 0) {
+        col2 = (n_col - 1) + copy_l;
+        break;
+      } else {
+        col1++;
+      }
+    }
+    /* Set pointers to columns */
+    help_matrix = matrix + n_row * col1;   // col1
+    help_matrix_2 = matrix + n_row * col2; // col2
+
+    for (int i = 0; i < n_row; i++) { // auxiliary vectors from matrix
+      U[i] = *(help_matrix + i) + *(help_matrix_2 + i);
+      V[i] = -*(help_matrix + i) + *(help_matrix_2 + i);
+    }
+    mad = quickselect_recursive(U, n_row, k);
+    med = quickselect_recursive(V, n_row, k);
+    if (evencorrect == 1) {
+      mad = (mad + quickselect_recursive(U, n_row, k - 1)) / 2;
+      med = (med + quickselect_recursive(V, n_row, k - 1)) / 2;
+    }
+    for (int i = 0; i < n_row; i++) { // reassign for new medians
+      U[i] = fabs(U[i] - mad);
+      V[i] = fabs(V[i] - med);
+    }
+    mad = quickselect_recursive(U, n_row, k);
+    med = quickselect_recursive(V, n_row, k);
+    if (evencorrect == 1) {
+      mad = (mad + quickselect_recursive(U, n_row, k - 1)) / 2;
+      med = (med + quickselect_recursive(V, n_row, k - 1)) / 2;
+    }
+    mad = pow(CONST * mad, 2);
+    med = pow(CONST * med, 2);
+
+    // Assign output
+    *(output + l) = (mad - med) / (mad + med);
+  }
+}
+#endif
diff --git a/src/cormaddp.f90 b/src/cormaddp.f90
deleted file mode 100644
index b0b6516..0000000
--- a/src/cormaddp.f90
+++ /dev/null
@@ -1,642 +0,0 @@
-module selectionalgo 
-   use, intrinsic :: iso_fortran_env
-
-   implicit none
-
-   real, parameter :: cost=1.4826, sqrt2=sqrt(2.) 
-   integer, parameter :: sp = real32, dp=real64
-
-   interface qselect
-      procedure quickselectdp, quickselectscalardp
-   end interface qselect
-
-   interface iselect
-      procedure introselectdp, introselectscalardp
-   end interface iselect
-
-   contains
-
-      !!!!!!! DOUBLE PRECISION
-      recursive subroutine quickselectrecursivedp(invector,tovector,k,output)
-         real(kind=dp), dimension(:), allocatable, target, intent(inout) :: invector
-         real(kind=dp),  intent(inout), pointer, contiguous :: tovector(:)
-         integer, intent(inout) :: k
-         real(kind=dp), intent(out) :: output
-         real(kind=dp) :: swapper, pvt
-         integer :: i, r, subst
-         real(kind=dp), pointer :: a,b,c
-
-         r=size(tovector)
-
-         select case(r)
-         case(1)
-            output=tovector(1)
-            return
-         case(2)
-            select case(k)
-               case(1)
-                  output=minval(tovector)
-               case(2)
-                  output=maxval(tovector)
-               end select
-            return
-         case(3)
-            select case(k)
-               case(1)
-                  output=minval(tovector)
-               case(3)
-                  output=maxval(tovector)
-               case(2)
-               i=r/2+1
-               output=sum(tovector)-minval(tovector)-maxval(tovector)
-            end select
-            return
-         case default
-            i=r/2+1
-
-            !pivoting section (pivot of 3)
-            a=>tovector(1)
-            b=>tovector(i)
-            c=>tovector(r)
-            pvt=a+b+c
-            swapper=max(a,b,c)
-            b=swapper
-            pvt=pvt-swapper
-            swapper=min(a,b,c)
-            a=swapper
-            pvt=pvt-swapper
-            c=pvt
-
-            !one pass section
-            subst=1
-            b=>tovector(subst)
-            do i=1,(r-1-3),3
-               a => tovector(i)
-               if (a<pvt) then
-                     swapper=b
-                     b=a
-                     a=swapper
-                  subst=subst+1
-                  b=>tovector(subst)
-               end if
-               a => tovector(i+1)
-               if (a<pvt) then
-                     swapper=b
-                     b=a
-                     a=swapper
-                  subst=subst+1
-                  b=>tovector(subst)
-               end if
-               a => tovector(i+2)
-               if (a<pvt) then
-                     swapper=b
-                     b=a
-                     a=swapper
-                  subst=subst+1
-                  b=>tovector(subst)
-               end if
-            end do
-            do i=i,(r-1)
-               a => tovector(i)
-               if (a<pvt) then
-                     swapper=b
-                     b=a
-                     a=swapper
-                  subst=subst+1
-                  b=>tovector(subst)
-               end if
-            end do
-            tovector(r)=tovector(subst)
-            tovector(subst)=pvt
-
-            !decide next step
-            select case(subst-k)
-            case(0)
-               output=tovector(k)
-            case(1:)
-               tovector => tovector(1:subst-1)
-               call quickselectrecursivedp(invector,tovector,k, output)
-            case(:-1)
-               k=k-subst
-               tovector => tovector(subst+1:)
-               call quickselectrecursivedp(invector,tovector,k, output)
-            end select
-         end select
-      end subroutine quickselectrecursivedp
-
-      recursive subroutine introselectrecursivedp(invector,tovector,k,output)
-         real(kind=dp), dimension(:), allocatable, target, intent(inout) :: invector
-         real(kind=dp),  intent(inout), pointer, contiguous :: tovector(:)
-         integer, intent(inout) :: k
-         real(kind=dp), intent(out) :: output
-         real(kind=dp) :: swapper, pvt
-         integer, save :: switch
-         integer :: i, r, subst, r_min, switch_after
-         real(kind=dp), pointer :: a,b,c
-
-         r=ubound(tovector,1)
-         r_min = 3000 ! at least this bigger to switch to median of medians
-         switch_after = 5 ! at least fails this amount of time to reorder
-                          ! one third of the vector to consider med of med
-
-         select case(r)
-         case(1)
-            output=tovector(1)
-            switch=0
-            return
-         case(2)
-            select case(k)
-               case(1)
-                  output=minval(tovector)
-               case(2)
-                  output=maxval(tovector)
-            end select
-            switch=0
-            return
-         case(3)
-            select case(k)
-               case(1)
-                  output=minval(tovector)
-               case(3)
-                  output=maxval(tovector)
-               case(2)
-               i=r/2+1
-               output=sum(tovector)-minval(tovector)-maxval(tovector)
-            end select
-            switch=0
-            return
-         case default
-            i=r/2+1
-
-            !pivoting section (pivot of 3)
-            a=>tovector(1)
-            b=>tovector(i)
-            c=>tovector(r)
-            pvt=a+b+c
-            swapper=max(a,b,c)
-            b=swapper
-            pvt=pvt-swapper
-            swapper=min(a,b,c)
-            a=swapper
-            pvt=pvt-swapper
-            c=pvt
-
-            !one pass section
-            subst=1
-            b=>tovector(subst)
-            do i=1,(r-1-3),3
-               a => tovector(i)
-               if (a<pvt) then
-                     swapper=b
-                     b=a
-                     a=swapper
-                  subst=subst+1
-                  b=>tovector(subst)
-               end if
-               a => tovector(i+1)
-               if (a<pvt) then
-                     swapper=b
-                     b=a
-                     a=swapper
-                  subst=subst+1
-                  b=>tovector(subst)
-               end if
-               a => tovector(i+2)
-               if (a<pvt) then
-                     swapper=b
-                     b=a
-                     a=swapper
-                  subst=subst+1
-                  b=>tovector(subst)
-               end if
-            end do
-            do i=i,(r-1)
-               a => tovector(i)
-               if (a<pvt) then
-                     swapper=b
-                     b=a
-                     a=swapper
-                  subst=subst+1
-                  b=>tovector(subst)
-               end if
-            end do
-            tovector(r)=tovector(subst)
-            tovector(subst)=pvt
-
-            !decide next step
-            select case(subst-k)
-            case(0)
-               output=tovector(k)
-               switch=0
-            case(1:)
-               tovector => tovector(1:subst-1)
-               if ((subst-1)>r*2/3) then
-                  switch=switch+1
-                  if (switch==switch_after .and. r>r_min) then
-                     !reset switch and call medofmeddp
-                     switch=0 
-                     call medofmeddp(invector,tovector,k, output)
-                  else
-                     call introselectrecursivedp(invector,tovector,k, output)
-                  end if
-               else
-                  switch=0
-                  call introselectrecursivedp(invector,tovector,k, output)
-               end if
-            case(:-1)
-               k=k-subst
-               tovector => tovector(subst+1:)
-               if ((r-subst)>r*2/3) then
-                  switch=switch+1
-                  if (switch==switch_after .and. r>r_min) then
-                     !reset switch and call medofmeddp
-                     switch=0 
-                     call medofmeddp(invector,tovector,k, output)
-                  else
-                     call introselectrecursivedp(invector,tovector,k, output)
-                  end if
-               else
-                  switch=0
-                  call introselectrecursivedp(invector,tovector,k, output)
-               end if
-            end select
-         end select
-      end subroutine introselectrecursivedp
-
-      recursive subroutine medofmeddp(invector,tovector,k,output)
-         real(kind=dp), dimension(:), allocatable, target, intent(inout) :: invector
-         real(kind=dp),  intent(inout), pointer, contiguous :: tovector(:)
-         real(kind=dp), intent(out) :: output
-         integer, intent(inout) :: k
-
-         real(kind=dp), pointer, contiguous :: subsec(:)
-         real(kind=dp),  dimension(:),allocatable :: meds
-         integer :: i,subst,r
-         real(kind=dp) :: pvt
-         !real(kind=dp) :: swapper
-         !real(kind=dp), pointer :: a,b
-
-         r=ubound(tovector,1)
-         select case(mod(r,5))
-            case(0)
-               allocate(meds(r/5))
-               do i=1,r-4,5
-                  subsec=>tovector(i:i+4)
-                  subst=minloc(subsec,1)
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(1)
-                  subsec(1)=pvt
-                  subst=minloc(subsec(2:),1)+1
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(2)
-                  subsec(2)=pvt
-                  subst=minloc(subsec(3:),1)+2
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(3)
-                  subsec(3)=pvt
-                  meds(int(i/5.0)+1)=pvt
-               end do
-            case(1)
-               allocate(meds(r/5+1))
-               do i=1,r-4,5
-                  subsec=>tovector(i:i+4)
-                  subst=minloc(subsec,1)
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(1)
-                  subsec(1)=pvt
-                  subst=minloc(subsec(2:),1)+1
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(2)
-                  subsec(2)=pvt
-                  subst=minloc(subsec(3:),1)+2
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(3)
-                  subsec(3)=pvt
-                  meds(int(i/5.0)+1)=pvt
-               end do
-               meds(r/5+1)=tovector(i)
-            case(2)
-               allocate(meds(r/5+1))
-               do i=1,r-4,5
-                  subsec=>tovector(i:i+4)
-                  subst=minloc(subsec,1)
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(1)
-                  subsec(1)=pvt
-                  subst=minloc(subsec(2:),1)+1
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(2)
-                  subsec(2)=pvt
-                  subst=minloc(subsec(3:),1)+2
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(3)
-                  subsec(3)=pvt
-                  meds(int(i/5.0)+1)=pvt
-               end do
-               subsec=>tovector(i:)
-               subst=minloc(subsec,1)
-               pvt=subsec(subst)
-               subsec(subst)=subsec(1)
-               subsec(1)=pvt
-               meds(r/5+1)=subsec(2)
-            case(3)
-               allocate(meds(r/5+1))
-               do i=1,r-4,5
-                  subsec=>tovector(i:i+4)
-                  subst=minloc(subsec,1)
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(1)
-                  subsec(1)=pvt
-                  subst=minloc(subsec(2:),1)+1
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(2)
-                  subsec(2)=pvt
-                  subst=minloc(subsec(3:),1)+2
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(3)
-                  subsec(3)=pvt
-                  meds(int(i/5.0)+1)=pvt
-               end do
-               subsec=>tovector(i:)
-               subst=minloc(subsec,1)
-               pvt=subsec(subst)
-               subsec(subst)=subsec(1)
-               subsec(1)=pvt
-               subst=minloc(subsec(2:),1)+1
-               pvt=subsec(subst)
-               subsec(subst)=subsec(2)
-               subsec(2)=pvt
-               meds(r/5+1)=pvt
-            case(4)
-               allocate(meds(r/5+1))
-               do i=1,r-4,5
-                  subsec=>tovector(i:i+4)
-                  subst=minloc(subsec,1)
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(1)
-                  subsec(1)=pvt
-                  subst=minloc(subsec(2:),1)+1
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(2)
-                  subsec(2)=pvt
-                  subst=minloc(subsec(3:),1)+2
-                  pvt=subsec(subst)
-                  subsec(subst)=subsec(3)
-                  subsec(3)=pvt
-                  meds(int(i/5.0)+1)=pvt
-               end do
-               subsec=>tovector(i:)
-               subst=minloc(subsec,1)
-               pvt=subsec(subst)
-               subsec(subst)=subsec(1)
-               subsec(1)=pvt
-               subst=minloc(subsec(2:),1)+1
-               pvt=subsec(subst)
-               subsec(subst)=subsec(2)
-               subsec(2)=pvt
-               subst=minloc(subsec(3:),1)+2
-               pvt=subsec(subst)
-               subsec(subst)=subsec(3)
-               subsec(3)=pvt
-               meds(r/5+1)=pvt
-            end select
-
-            !set pivot
-            pvt=qselect(meds)
-            do i=1,r-3,3
-               if (tovector(i)==pvt) then
-                  tovector(i)=tovector(r)
-                  tovector(r)=pvt
-                  exit
-               else if (tovector(i+1)==pvt) then
-                  tovector(i+1)=tovector(r)
-                  tovector(r)=pvt
-                  exit
-               else if (tovector(i+2)==pvt) then
-                  tovector(i+2)=tovector(r)
-                  tovector(r)=pvt
-                  exit
-               else 
-                  cycle
-               end if
-            end do
-            do i=i,r
-               if (tovector(i)==pvt) then
-                  tovector(i)=tovector(r)
-                  tovector(r)=pvt
-                  exit
-               else 
-                  cycle
-               end if
-            end do
-
-            call introselectrecursivedp(invector,tovector,k, output)
-      end subroutine medofmeddp
-
-      function introselectdp(invector,ord,evencorrection)
-         real(kind=dp), intent(in), dimension(:) :: invector
-         logical, optional, intent(in) :: evencorrection
-         integer, intent(in), optional :: ord
-         real(kind=dp) :: introselectdp
-
-         real(kind=dp), dimension(:), allocatable, target :: vector
-         real(kind=dp), pointer, contiguous :: tovector(:) => null()
-         integer ::  k1, k
-
-         k1=size(invector)
-         allocate(vector(k1))
-         vector=invector
-         tovector => vector
-
-         k1=ubound(invector,1)
-         if(present(ord)) then
-            k=ord
-         else
-            k=k1/2+1
-         end if
-
-         if(present(evencorrection)) then
-            if (k/=k1/2+1 .or. mod(k1,2)/=0) then
-               !write(*,'(a)') 'Warning: evencorrection argument ignored.'
-               k1=0
-            else 
-               k1=0
-               if (evencorrection) k1=k-1
-            end if
-         else
-            k1=0
-         end if
-
-         call introselectrecursivedp(vector,tovector,k, introselectdp)
-
-         if (k1/=0) then
-            block
-               real(kind=dp) :: temp
-               temp=introselectdp
-               !tovector=>vector(1:k1+1)
-               tovector=>vector
-               call introselectrecursivedp(vector,tovector,k1, introselectdp)
-               introselectdp=(introselectdp+temp)/2.0
-            end block
-         end if
-      end function introselectdp
-
-      function quickselectdp(invector,ord,evencorrection)
-         real(kind=dp), intent(in), dimension(:) :: invector
-         logical, optional, intent(in) :: evencorrection
-         integer, intent(in), optional :: ord
-         real(kind=dp) :: quickselectdp
-
-         real(kind=dp), dimension(:), allocatable, target :: vector
-         real(kind=dp), pointer, contiguous :: tovector(:) => null()
-         integer ::  k1, k
-
-         k1=size(invector)
-         allocate(vector(k1))
-         vector=invector
-         tovector => vector
-
-         k1=ubound(invector,1)
-         if(present(ord)) then
-            k=ord
-         else
-            k=k1/2+1
-         end if
-
-         if(present(evencorrection)) then
-            if (k/=k1/2+1 .or. mod(k1,2)/=0) then
-               !write(*,'(a)') 'Warning: evencorrection argument ignored.'
-               k1=0
-            else 
-               k1=0
-               if (evencorrection) k1=k-1
-            end if
-         else
-            k1=0
-         end if
-
-         call quickselectrecursivedp(vector,tovector,k, quickselectdp)
-
-         if (k1/=0) then
-            block
-               real(kind=dp) :: temp
-               temp=quickselectdp
-               tovector=>vector(1:k1+1)
-               call quickselectrecursivedp(vector,tovector,k1, quickselectdp)
-               quickselectdp=(quickselectdp+temp)/2.0
-            end block
-         end if
-      end function quickselectdp
-
-      !deal with scalar input
-      function introselectscalardp(invector,ord, evencorrection) 
-         real(kind=dp), intent(in) :: invector
-         integer, intent(in), optional :: ord
-         logical, optional, intent(in) :: evencorrection
-         real(kind=dp) :: introselectscalardp
-
-         if (present(evencorrection) .or. present(ord)) then
-            continue
-         end if
-
-         introselectscalardp=invector
-      end function introselectscalardp
-
-      function quickselectscalardp(invector,ord, evencorrection) 
-         real(kind=dp), intent(in) :: invector
-         integer, intent(in), optional :: ord
-         logical, optional, intent(in) :: evencorrection
-         real(kind=dp) :: quickselectscalardp
-
-         if (present(evencorrection) .or. present(ord)) then
-            continue
-         end if
-
-         quickselectscalardp=invector
-      end function quickselectscalardp
-
-
-end module selectionalgo
-
-subroutine cormadvecdp(matrix,nrow,ncol,res,ressize,correcteven)
-
-   use selectionalgo
-
-   implicit none
-   
-   !note kind dp is defined in selectionalgo
-   integer, intent(in) :: nrow, ncol,ressize, correcteven     !.Fortran R function will not deal with deferred size arrays
-   real(kind=dp), dimension(nrow,ncol), intent(inout) :: matrix
-   real(kind=dp), dimension(ressize), intent(out) :: res
-
-   integer :: i, j, n, p
-   real(kind=dp)  :: med,mad
-   !real(kind=dp) :: fresh
-   real(kind=dp), dimension(:), allocatable :: U, V
-   !real(kind=dp), dimension(:), allocatable :: A, B
-
-
-   p=ubound(matrix,2)
-   n=ubound(matrix,1)
-   allocate(U(n),V(n))
-
-   if (correcteven==1) then
-      do i=1,p-3,3
-         U=matrix(:,i)
-         med=iselect(U,evencorrection=.true.)
-         matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.))
-         U=matrix(:,i+1)
-         med=iselect(U,evencorrection=.true.)
-         matrix(:,i+1)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.))
-         U=matrix(:,i+2)
-         med=iselect(U,evencorrection=.true.)
-         matrix(:,i+2)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.))
-      end do
-      do i=i,p
-         U=matrix(:,i)
-         med=iselect(U,evencorrection=.true.)
-         matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.))
-      end do
-      !unrolled loops
-      n=1
-      do i=1,p-1
-         do j=i+1,p
-            U=matrix(:,i)+matrix(:,j)
-            V=-matrix(:,i)+matrix(:,j)
-            mad=(cost*iselect(abs(U-iselect(U,evencorrection=.true.)),evencorrection=.true.))**2
-            med=(cost*iselect(abs(V-iselect(V,evencorrection=.true.)),evencorrection=.true.))**2
-            res(n)=(mad-med)/(mad+med)
-            n=n+1
-         end do 
-      end do
-   else
-      do i=1,p-3,3
-         U=matrix(:,i)
-         med=iselect(U)
-         matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med)))
-         U=matrix(:,i+1)
-         med=iselect(U)
-         matrix(:,i+1)=(U-med)/(sqrt2*cost*iselect(abs(U-med)))
-         U=matrix(:,i+2)
-         med=iselect(U)
-         matrix(:,i+2)=(U-med)/(sqrt2*cost*iselect(abs(U-med)))
-      end do
-      do i=i,p
-         U=matrix(:,i)
-         med=iselect(U)
-         matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med)))
-      end do
-      !unrolled loops
-      n=1
-      do i=1,p-1
-         do j=i+1,p
-            U=matrix(:,i)+matrix(:,j)
-            V=-matrix(:,i)+matrix(:,j)
-            mad=(cost*iselect(abs(U-iselect(U))))**2
-            med=(cost*iselect(abs(V-iselect(V))))**2
-            res(n)=(mad-med)/(mad+med)
-            n=n+1
-         end do 
-      end do
-   end if
-end subroutine cormadvecdp
diff --git a/src/init.c b/src/init.c
index 66209f7..32fd0f1 100644
--- a/src/init.c
+++ b/src/init.c
@@ -1,22 +1,66 @@
-#include <stdlib.h>
-//#include <stdio.h>
+#include "RSCdefines.h"
 #include <R.h>
-#include <Rinternals.h>
-#include <Rmath.h>
-#include <Rdefines.h>
 #include <R_ext/RS.h>
 #include <R_ext/Rdynload.h>
+#include <Rdefines.h>
+#include <Rinternals.h>
+#include <Rmath.h>
+#include <stdio.h>
+#include <stdlib.h>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+// Wrapper for C cormad (with pointers)
+SEXP cormad_C(SEXP R_matrix, SEXP len_rows, SEXP len_cols, SEXP R_evencorrect,
+              SEXP R_num_threads) {
+  // prepare input arguments
+  int n_row = asInteger(len_rows);
+  int n_col = asInteger(len_cols);
+  int output_size = (n_col - 1) * n_col / 2;
+  int evencorrect =
+      (n_row % 2 == 1) ? 0 : asInteger(R_evencorrect); // ingore for odd n
+  int num_threads = asInteger(R_num_threads);
+  SEXP output = PROTECT(allocVector(REALSXP, output_size)); // define output
 
-extern void F77_NAME(cormadvecdp)(void *, void *, void *, void *, void *, void *);
+  // create pointers to (duplicate) matrix column
+  double *dupmat =
+      REAL(PROTECT(duplicate(R_matrix))); // copy input (to avoid modifications)
+
+  // Call cormad
+#ifdef _OPENMP
+  int max_threads = omp_get_max_threads();
+  if (num_threads == 0) /*  use max number of threads / 2 */
+    num_threads = (max_threads / 2 == 0) ? 1 : max_threads / 2;
+  else if (num_threads > max_threads)
+    num_threads = max_threads;
+  else if (num_threads < 0)
+    num_threads =
+        (max_threads + num_threads > 0) ? max_threads + num_threads : 1;
+
+  if (num_threads > 1) {
+    cormad_parallel(dupmat, n_row, n_col, REAL(output), evencorrect,
+                    num_threads);
+  } else {
+    cormad(dupmat, n_row, n_col, REAL(output), evencorrect);
+  }
+#else
+  if (num_threads != 1)
+    Rprintf("\nSpecified num_threads=%d, but OPENMP support not found; "
+            "switching to single core.\n\n",
+            num_threads);
+  cormad(dupmat, n_row, n_col, REAL(output), evencorrect);
+#endif
+  // Remove protect to allow R grabage collect
+  UNPROTECT(2);
+  return output;
+}
 
-static const R_FortranMethodDef FortMethods[] = {
-   {"cormadvecdp", (DL_FUNC) &F77_NAME(cormadvecdp), 6},
-   {NULL, NULL, 0}
-};
+static const R_CallMethodDef callMethods[] = {
+    {"cormad_C", (DL_FUNC)&cormad_C, 5}, {NULL, NULL, 0}};
 
-void
-R_init_RSC(DllInfo *dll)
-{
-   R_registerRoutines(dll, NULL, NULL, FortMethods, NULL);
-   R_useDynamicSymbols(dll, FALSE);
+void R_init_RSC(DllInfo *dll) {
+  R_registerRoutines(dll, NULL, callMethods, NULL, NULL);
+  R_useDynamicSymbols(dll, FALSE);
+  R_forceSymbols(dll, TRUE);
 }
diff --git a/src/selection_algos_ptr.c b/src/selection_algos_ptr.c
new file mode 100644
index 0000000..2abecac
--- /dev/null
+++ b/src/selection_algos_ptr.c
@@ -0,0 +1,154 @@
+#include <stdio.h>
+/*  Implements functions quickselect and introselect */
+
+static double select_corner_cases(double *vector_shorter_3, int size_of_vector,
+                                  int extract_this_element_Ccorrected);
+
+static double pivot_of_3(double *vector_to_be_pivoted, int size);
+
+static int quickselect_onepass(double *vector, int size);
+
+/* Functions */
+
+double quickselect_recursive(double *vector, int size, int k) {
+  /* Arguments:
+   * *vector: pointer to vector
+   * idN: size of vector
+   * k: k-th element to be extracted from the vector output: the
+   * extracted element of order k */
+
+  /* Decide if we need to go with quickselect_recursive or return output
+   * (quickselect_recursive works only for vectors with at least 3 elems) */
+  if (size < 3) { // select the returning value
+    return select_corner_cases(vector, size, k);
+  }
+
+  int subst = quickselect_onepass(vector, size);
+
+  /* Determine where to go next (left of subst or right?) */
+  if (subst == k) {
+    return *(vector + subst);
+  } else if (subst > k) { // go left
+    size = subst;
+    return quickselect_recursive(vector, size, k);
+  } else {          // go right
+                    // readjust k
+    k -= subst + 1; // take it back to orginal value (for full vector)
+    vector += subst + 1;
+    size -= subst + 1;
+    return quickselect_recursive(vector, size, k);
+  }
+}
+
+static double select_corner_cases(double *vector, int size, int k) {
+  /* Used to handel corner cases not handeled by recursive strategy.
+   * Recursive strategy needs vectors of at least 3 elements
+
+   *  Arguments:
+   * *vector: pointer to vector of doubles to sort
+   * size: size of the vector
+   * k: k-th element to be extracted from the vector output: the
+   * extracted element of order k
+  */
+
+  // total elements in vector are idN-id0+1
+  double ret = -111;
+  switch (size) {
+  case 1:
+    return *vector;
+  case 2:
+    switch (k) {
+    case 0: // return the smallest of the two
+      if (*vector < *(vector + 1))
+        return *vector;
+      else
+        return *(vector + 1);
+    case 1: // return the biggest of the two
+      if (*vector > *(vector + 1))
+        return *vector;
+      else
+        return *(vector + 1);
+    }
+  }
+  return ret;
+}
+
+static double pivot_of_3(double *vector, int size) {
+  /* perform pivot of 3 returning the median element and
+   * rearranging the vector so to have:
+   * id0, i(see below) , idN -> min, max, med
+
+   * Arguments
+   * vector: vector where to make substitutions
+   * id0, i, idN: positions of first, median and last element to consider
+  */
+
+  int idN = size - 1; // index of last element of the vector
+  int i = idN / 2;    // index of middle element of the vector
+  double a, b, c;     // auxiliary variables (avoid dereferncing too much)
+  a = *vector;
+  b = *(vector + i);
+  c = *(vector + idN);
+
+  /* Nota; vogliamo che l'elemento mediano si trovi alla fine */
+  double swapper;
+  if ((a > b) ^ (a > c)) { // id0 is median element
+    swapper = c;
+    c = a;
+    if (swapper > b) {
+      a = b;
+      b = swapper;
+    } else {
+      a = swapper;
+    }
+  } else if ((b > a) ^ (b > c)) { // i is median element
+    swapper = c;
+    c = b;
+    if (swapper > a) {
+      b = swapper;
+    } else {
+      b = a;
+      a = swapper;
+    }
+  } else { // idN is median element
+    if (a > b) {
+      swapper = a;
+      a = b;
+      b = swapper;
+    }
+  }
+  *vector = a;
+  *(vector + i) = b;
+  *(vector + idN) = c;
+  return c;
+}
+
+static int quickselect_onepass(double *vector, int size) {
+  /* perform one pass for quickselect and returns the element where the pivot
+   * was substituted: this is used to decide whether to go left or right */
+
+  /* Pivoting section (pivot of 3) */
+  double pivot = pivot_of_3(vector, size);
+
+  /* One pass on vector */
+  double swapper;
+  double *subst_el = vector; // address of element to substitute (copy not to
+                             // alter main pointer)
+  int subst = 0;
+  for (int i = 0; i < size - 1; i++) {
+    if (*subst_el < pivot) {
+      if (subst == i) {
+        subst++;
+      } else {
+        swapper = *subst_el;
+        *subst_el = *(vector + subst);
+        *(vector + subst) = swapper;
+        subst++;
+      }
+    }
+    subst_el++;
+  }
+  *subst_el = *(vector + subst);
+  *(vector + subst) = pivot;
+  return subst;
+}