version 1.0

19 files changed, 2237 insertions, 0 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..36f47f6
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,27 @@
+Package: RSC
+Type: Package
+Title: Robust and Sparse Correlation Matrix
+Description: Performs robust and sparse correlation matrix estimation. Robustness is achieved based on a simple robust pairwise correlation estimator, while sparsity is obtained based on thresholding. The optimal thresholding is tuned via cross-validation. See Serra, Coretto, Fratello, and Tagliaferri (2018) <doi:10.1093/bioinformatics/btx642>.
+Authors@R: c(person("Luca", "Coraggio",
+                    role = c("cre", "aut"), email = "luca.coraggio@unina.it"),
+             person("Pietro", "Coretto",
+	            role = c("aut"), email = "pcoretto@unisa.it"),
+             person("Angela", "Serra",
+	            role = c("aut"), email = "angela.serra@tuni.fi"),
+             person("Roberto", "Tagliaferri",
+	            role = c("ctb"), email = "robtag@unisa.it")) 
+Author: Luca Coraggio [cre, aut],
+  Pietro Coretto [aut],
+  Angela Serra [aut],
+  Roberto Tagliaferri [ctb]
+Maintainer: Luca Coraggio <luca.coraggio@unina.it>
+NeedsCompilation: yes
+Imports: stats, graphics, Matrix, methods, parallel, foreach,
+        doParallel, utils
+License: GPL (>= 2)
+LazyData: TRUE
+Version: 1.0
+Date: 2020-06-29
+Packaged: 2020-07-02 07:30:56 UTC; pietro
+Repository: CRAN
+Date/Publication: 2020-07-04 10:50:03 UTC
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..f14927a
--- /dev/null
+++ b/MD5
@@ -0,0 +1,18 @@
+80c587bf8be31b19549b965aa44f1d50 *DESCRIPTION
+1f64d30a870f2eae4855e83010aaaa4e *NAMESPACE
+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
+ee184b1708ad4ef4735070065fe6737c *inst/CITATION
+fffcd92fdbda643ea59f12298cc04287 *man/plot.cv_rsc.Rd
+57435eb3ed73e1af9435073eeaf57928 *man/rmad.Rd
+18e1373418dce90153d0727d19f5fa9e *man/rsc.Rd
+2ab80b72b7090b2514c9cf4b3fbe006d *man/rsc_cv.Rd
+24cec7b34e7c1931f74d1f25cad05a83 *src/Makevars
+f135feaeda1c3e803c74729a94fcc416 *src/cormaddp.f90
+9ded1c66a2fd469e79d7d1d25d0bd604 *src/init.c
+62f9244cb578a6720b664233da503e32 *src/selectiongeneral.f95
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..ffd8478
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,25 @@
+## import all functions from the following pkgs
+import("stats", "graphics", "Matrix")
+
+## selective imports from pkgs 
+importFrom("methods", "as")
+importFrom("foreach", "foreach", "%dopar%")
+importFrom("parallel","detectCores")
+importFrom("doParallel", "registerDoParallel", "stopImplicitCluster")
+importFrom("utils", "citHeader", "citEntry")
+
+## export objects in /src 
+useDynLib(RSC, .registration = TRUE, .fixes = "F_")
+
+## export objects in /R
+export(rmad)
+export(rsc_cv)
+export(rsc)
+
+
+## export S3 methods
+S3method(plot,  rsc_cv)
+S3method(print, rsc_cv)
+
+## ## exports every object that doesn't start with a dot
+## exportPattern("^[^\\.]")
diff --git a/R/check_inputs.R b/R/check_inputs.R
new file mode 100644
index 0000000..4cbcf0a
--- /dev/null
+++ b/R/check_inputs.R
@@ -0,0 +1,64 @@
+.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")
+        }
+      }
+    }
+    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)
+}
diff --git a/R/cv_loss.R b/R/cv_loss.R
new file mode 100644
index 0000000..ea46007
--- /dev/null
+++ b/R/cv_loss.R
@@ -0,0 +1,17 @@
+.cv_loss <- function(idx, dat, evencorrection, threshold, grid.length, p, nc) {
+  res <- numeric(nc)
+  n1 <- as.integer(sum(idx))
+  C1 <- .Fortran("cormadvecdp", matrix = dat[idx, ], nrow = n1, ncol = p, res = res, 
+    ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res
+  n2 <- as.integer(sum(!idx))
+  C2 <- .Fortran("cormadvecdp", matrix = dat[!idx, ], nrow = n2, ncol = p, res = res, 
+    ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res
+  ans <- rep(0, times = grid.length)
+  for (h in 1:grid.length) {
+    C1[abs(C1) < threshold[h]] <- 0
+    ans[h] <- sum(2 * {
+      C1 - C2
+    }^2)/p
+  }
+  return(ans)
+}
diff --git a/R/plot_print_methods.R b/R/plot_print_methods.R
new file mode 100644
index 0000000..ca163cb
--- /dev/null
+++ b/R/plot_print_methods.R
@@ -0,0 +1,26 @@
+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)
+}
diff --git a/R/rmad.R b/R/rmad.R
new file mode 100644
index 0000000..b8e9d41
--- /dev/null
+++ b/R/rmad.R
@@ -0,0 +1,31 @@
+rmad <- function(x, y = NULL, na.rm = FALSE, even.correction = FALSE) {
+  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)
+  if (even.correction) {
+    evencorrection <- 1L
+  }
+  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)
+  }
+  else {
+    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")
+    if (!is.null(colnames_original)) {
+      dimnames(R)[[1]] <- dimnames(R)[[2]] <- colnames_original
+    }
+    return(R)
+  }
+}
diff --git a/R/rsc.R b/R/rsc.R
new file mode 100644
index 0000000..6cf4e15
--- /dev/null
+++ b/R/rsc.R
@@ -0,0 +1,42 @@
+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)
+}
diff --git a/R/rsc_cv.R b/R/rsc_cv.R
new file mode 100644
index 0000000..4153322
--- /dev/null
+++ b/R/rsc_cv.R
@@ -0,0 +1,183 @@
+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) {
+  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)
+  if ({
+    cv.type != "random"
+  } & {
+    cv.type != "kfold"
+  }) {
+    stop("\"cv.type\" must be either \"random\" (default)  or \"kfold\"")
+  }
+  if (!is.numeric(R)) {
+    stop("\"R\" must be an integer > 1")
+  }
+  else if (R < 1) {
+    stop("\"R\" must be an integer > 1")
+  }
+  if (!is.numeric(K)) {
+    stop("\"K\" must be an integer > 1")
+  }
+  else if (R < 1) {
+    stop("\"K\" must be an integer > 1")
+  }
+  if (length(threshold) == 1) {
+    if (threshold <= 0 | threshold >= 1) {
+      stop("\"threshold\" value does not belong to the interval (0,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).")
+    }
+    grid.length <- length(threshold)
+  }
+  if (even.correction) {
+    evencorrection <- 1L
+  }
+  else {
+    evencorrection <- 0L
+  }
+  if (is.null(ncores)) {
+    DetectedCores <- detectCores()
+    if (DetectedCores <= 2) {
+      ncores <- 1
+    }
+    else {
+      ncores <- {
+        DetectedCores - 1
+      }
+    }
+  }
+  else {
+    ncores <- as.integer(ncores)
+    if (ncores <= 0) {
+      stop("\"ncores\" must be an integer larger or equal to 1.")
+    }
+  }
+  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
+  if (monitor) {
+    t1 <- Sys.time()
+    dt01 <- difftime(t1, t0, units = "auto")
+    message("* RMAD computing time:...... ", round(dt01, 2), " [", attributes(dt01)$units, 
+      "]")
+  }
+  if (cv.type == "random") {
+    if (monitor) {
+      cat("\n")
+      message("Performing cross-validation")
+      t_hat <- round(1.2 * {
+        {
+          dt01 * R * 2
+        }/ncores
+      }, 2)
+      message("* predicted end time (worst case):...... ", Sys.time() + t_hat)
+    }
+    n1 <- n - floor(n/log(n))
+    IDX <- array(FALSE, dim = c(R, n))
+    for (r in 1:R) {
+      IDX[r, ][sample(1:n, size = n1, replace = FALSE)] <- TRUE
+    }
+    registerDoParallel(ncores)
+    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)
+    }
+    stopImplicitCluster()
+    LOSS <- array(0, dim = c(R, grid.length))
+    for (r in 1:R) {
+      LOSS[r, ] <- U[[r]]
+    }
+    avg_loss <- apply(LOSS, 2, mean)
+    se_loss <- apply(LOSS, 2, sd)/sqrt(R)
+  }
+  if (cv.type == "kfold") {
+    if (monitor) {
+      cat("\n")
+      message("Performing cross-validation")
+      t_hat <- round(1.2 * {
+        {
+          dt01 * R * K * 2
+        }/ncores
+      }, 2)
+      message("* predicted end time (worst case):...... ", Sys.time() + t_hat)
+    }
+    idx_fold <- cut(1:n, breaks = K, labels = FALSE)
+    IDX <- array(TRUE, dim = c(R * K, n))
+    row_count <- 1L
+    for (r in 1:R) {
+      idx_fold_shuffle <- sample(idx_fold, size = n, replace = FALSE)
+      for (k in 1:K) {
+        IDX[row_count, ][idx_fold_shuffle == k] <- FALSE
+        row_count <- 1L + row_count
+      }
+    }
+    registerDoParallel(ncores)
+    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)
+    }
+    stopImplicitCluster()
+    if (R == 1) {
+      LOSS <- array(0, dim = c(K, grid.length))
+      for (k in 1:K) {
+        LOSS[k, ] <- U[[k]]
+      }
+    }
+    else {
+      LOSS <- array(0, dim = c(K, grid.length, R))
+      dimnames(LOSS)[[3]] <- paste0("r", 1:R)
+      dimnames(LOSS)[[1]] <- paste0("k", 1:K)
+      row_count <- 1L
+      for (r in 1:R) {
+        for (k in 1:K) {
+          LOSS[k, , r] <- U[[row_count]]
+          row_count <- 1L + row_count
+        }
+      }
+    }
+    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)
+    }
+    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")
+    }
+  }
+  if (monitor) {
+    message("* finished on:.......................... ", Sys.time())
+  }
+  tstar <- which.min(avg_loss)
+  flags <- rep("", grid.length)
+  a <- avg_loss[tstar] - se_loss[tstar]
+  b <- avg_loss[tstar] + se_loss[tstar]
+  flags[{
+    avg_loss >= a
+  } & {
+    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]))
+  class(ans) <- "rsc_cv"
+  return(ans)
+}
diff --git a/R/zzz.R b/R/zzz.R
new file mode 100755
index 0000000..3c33499
--- /dev/null
+++ b/R/zzz.R
@@ -0,0 +1,4 @@
+.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/inst/CITATION b/inst/CITATION
new file mode 100644
index 0000000..8690fa8
--- /dev/null
+++ b/inst/CITATION
@@ -0,0 +1,40 @@
+citHeader("To cite package 'RSC' in publications use:")
+
+year <- sub("-.*", "", meta$Date)
+note <- sprintf("R package version %s", meta$Version)
+
+citEntry(entry="Manual",
+         title   = "RSC: robust and sparse correlation matrix estimation",
+         author  = personList(person(given="Luca",    family="Coraggio"),
+                              person(given="Pietro",  family="Coretto"),
+                              person(given="Angela",  family="Serra"),
+                              person(given="Roberto", family="Tagliaferri")),
+         year    = year,
+         note    = note,
+         textVersion =
+         paste("Coraggio, L., Coretto, P., Serra, A. and Tagliaferri, R. (", year,").
+                RSC: robust and sparse correlation matrix estimation. ",
+                note, "url: https://CRAN.R-project.org/package=rsc",
+                sep='')
+        )
+
+
+
+citEntry(entry="Article",
+         title   = "Robust and sparse correlation matrix estimation for the analysis
+                    of high-dimensional genomics data",
+         author  = personList(person(given="Angela", family="Serra"),
+                              person(given="Pietro", family="Coretto"),
+                              person(given="Michele", family="Fratello"),
+                              person(given="Roberto", family="Tagliaferri")
+                              ),
+         year    = 2018,
+         journal = "Bioinformatics",
+         volume  = 34,
+         number  = 4,
+         pages   = "625--634",
+         doi     = "10.1093/bioinformatics/btx642",
+         textVersion = "Serra, A., Coretto, P., Fratello, M., and Tagliaferri, R. (2018). Robust and sparse correlation matrix estimation for the analysis of high-dimensional genomics data. Bioinformatics, 34(4), 625-634."
+         )
+
+
diff --git a/man/plot.cv_rsc.Rd b/man/plot.cv_rsc.Rd
new file mode 100755
index 0000000..0d7b2c2
--- /dev/null
+++ b/man/plot.cv_rsc.Rd
@@ -0,0 +1,84 @@
+\name{plot.rsc_cv}
+\alias{plot.rsc_cv}
+
+\title{
+  Plot method for rsc_cv objects
+}
+
+\description{
+  Plot the cross-validation estimates of the Frobenius loss.
+}
+
+
+\usage{
+   \method{plot}{rsc_cv}(x, \dots)
+}
+
+
+
+
+\arguments{
+  \item{x}{
+    Output from \code{\link{rsc_cv}}, that is an S3 object of class \code{"rsc_cv"}.
+  }
+  \item{\dots}{
+    additional arguments passed to  \code{\link[graphics]{plot.default}}.
+  }
+}
+
+
+\value{
+   Plot the Frobenius loss estimated via cross-validation (y-axis) vs
+   threshold values (x-axis). The dotted blue line represents the average
+   expected normalized Frobenius loss, while the vertical segments
+   around the average  are  \emph{1-standard-error} error bars
+   (see \emph{Details} in \code{\link{rsc_cv}}.
+
+   The vertical dashed red line identifies the minimum of the average
+   loss, that is the optimal threshold flagged as \code{"minimum"}.  The
+   vertical dashed green line identifies the optimal selection flagged
+   as \code{"minimum1se"} in the output of \code{\link{rsc_cv}} (see
+   \emph{Details} in \code{\link{rsc_cv}}).
+}
+
+
+
+
+\section{References}{
+   Serra, A., Coretto, P., Fratello, M., and Tagliaferri, R. (2018).
+   Robust and sparsecorrelation matrix estimation for the analysis of
+   high-dimensional genomics data. \emph{Bioinformatics}, 34(4),
+   625-634.  doi:10.1093/bioinformatics/btx642
+}
+
+
+\seealso{
+  \code{\link{rsc_cv}}
+}
+
+
+
+\examples{
+\donttest{
+## simulate a random sample from a multivariate Cauchy distribution
+## note: example in high-dimension are obtained increasing p
+set.seed(1)
+n   <- 100  # sample size
+p   <- 10   # dimension
+dat <- matrix(rt(n*p, df = 1), nrow = n, ncol = p)
+colnames(dat) <- paste0("Var", 1:p)
+
+   
+## perform 10-fold cross-validation repeated R=10 times
+## note: for multi-core machines experiment with 'ncores'
+set.seed(2)
+a <- rsc_cv(x = dat, R = 10, K = 10)
+a
+
+## plot the cross-validation estimates
+plot(a)
+   
+## pass additional parameters to graphics::plot
+plot(a , cex = 2)
+}
+}
diff --git a/man/rmad.Rd b/man/rmad.Rd
new file mode 100755
index 0000000..01247cf
--- /dev/null
+++ b/man/rmad.Rd
@@ -0,0 +1,136 @@
+\name{rmad}
+
+\alias{rmad}
+
+\title{RMAD correlation matrix}
+
+\description{
+   Compute the RMAD robust correlation matrix proposed in Serra et
+   al. (2018) based on the robust correlation coefficient proposed in
+   Pasman and Shevlyakov (1987).
+}
+
+
+\usage{
+  rmad(x , y = NULL, na.rm = FALSE , even.correction = FALSE)   
+}
+
+
+\arguments{
+  \item{x}{
+    A numeric vector, a matrix or a data.frame. If \code{x} is a matrix
+    or a data.frame, rows of  \code{x} correspond to sample units
+    and columns correspond to variables. If \code{x} is a numerical
+    vector, and \code{y} is not \code{NULL}, the RMAD correlation
+    coefficient between \code{x} and \code{y} is computed. Categorical
+    variables are not allowed.
+  }
+  \item{y}{
+    A numerical vector if not \code{NULL}. If both \code{x} and \code{y}
+    are numerical vectors, the RMAD correlation coefficient between
+    \code{x} and \code{y} is computed.
+  }
+  \item{na.rm}{
+    A logical value, if \code{TRUE} sample observation
+    containing \code{NA} values are excluded (see \emph{Details}).
+  }
+  \item{even.correction}{
+    A logical value, if \code{TRUE} a correction
+    for the calculation of the medians is applied to reduce the bias
+    when the number of samples even (see \emph{Details}).
+  }
+}
+
+
+\details{
+   The \code{rmad} function computes the correlation matrix based on the
+   pairwise robust correlation coefficient of Pasman and Shevlyakov
+   (1987). This correlation coefficient is based on repeated median
+   calculations for all pairs of variables. This is a computational
+   intensive task when the number of variables (that is \code{ncol(x)})
+   is large.
+   
+   The software is optimized for large dimensional data sets, the median
+   is approximated as the central observation obtained based on the
+   \emph{introselect} sorting algorithm of Musser (1997) implemented in
+   Fortran 95 language. For small samples this may be a crude
+   approximation, however, it makes the computational cost feasible for
+   high-dimensional data sets. With the option \code{even.correction
+   = TRUE} a correction is applied to reduce the bias for data sets with
+   an even number of samples. Although \code{even.correction = TRUE}
+   has a small computational cost for each pair of variables, it is
+   suggested to use the default \code{even.correction = FALSE} for large
+   dimensional data sets.
+
+   The function can handle a data matrix with missing values (\code{NA}
+   records). If \code{na.rm = TRUE} then missing values are handled by
+   casewise deletion (and if there are no complete cases, an error is
+   returned). In practice, if \code{na.rm = TRUE} all rows of
+   \code{x} that contain at least an \code{NA} are removed.
+
+   Since the software is optimized to work with high-dimensional data sets,
+   the output RMAD matrix is packed into a storage efficient format
+   using the \code{"dspMatrix"} S4 class from the \code{\link{Matrix}}
+   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}}.
+}
+
+
+
+\value{
+  \item{If \code{x} is a matrix or a data.frame}{
+    Returns a correlation matrix of class \code{"dspMatrix"} (S4 class object)
+    as defined in the \code{\link{Matrix}} package.
+  }
+  \item{If \code{x} and \code{y} are numerical vectors}{
+    Returns a numerical value, that is the RMAD correlation coefficient
+    between \code{x} and \code{y}.    
+  }
+}
+
+
+
+\section{References}{
+  Musser, D. R. (1997). Introspective sorting and selection algorithms.
+  \emph{Software: Practice and Experience}, 27(8), 983-993.
+  
+  Pasman,V. and Shevlyakov,G. (1987). Robust methods of estimation of
+  correlation coefficient. \emph{Automation Remote Control}, 48, 332-340.
+
+  Serra, A., Coretto, P., Fratello, M., and Tagliaferri, R. (2018).
+  Robust and sparsecorrelation matrix estimation for the analysis of
+  high-dimensional genomics data. \emph{Bioinformatics}, 34(4), 625-634.
+  doi: 10.1093/bioinformatics/btx642 
+}
+
+
+
+\seealso{
+  \code{rsc_cv}, \code{rsc}
+}
+
+
+
+
+
+
+
+\examples{
+## simulate a random sample from a multivariate Cauchy distribution
+set.seed(1)
+n   <- 100    # sample size
+p   <- 7      # dimension
+dat <- matrix(rt(n*p, df = 1), nrow = n, ncol = p)
+colnames(dat) <- paste0("Var", 1:p)
+
+   
+## compute the rmad correlation coefficient between dat[,1] and dat[,2]
+a <- rmad(x = dat[,1], y = dat[,2])
+
+
+## compute the RMAD correlaiton matrix   
+b <- rmad(x = dat)
+b
+}
+
diff --git a/man/rsc.Rd b/man/rsc.Rd
new file mode 100755
index 0000000..cd252a6
--- /dev/null
+++ b/man/rsc.Rd
@@ -0,0 +1,118 @@
+\name{rsc}
+
+\alias{rsc}
+
+\title{Robust and Sparse Correlation Matrix Estimator}
+
+\description{
+  Compute the Robust and Sparse Correlation Matrix (RSC) estimator 
+  proposed in Serra et al. (2018).
+}
+
+
+\usage{
+  rsc(cv, threshold = "minimum") 
+}
+
+
+\arguments{
+  \item{cv}{
+    An S3 object of class \code{"rsc_cv"} (see \code{\link{rsc_cv}}).
+  }
+  \item{threshold}{
+    Threshold parameter to compute the RSC estimate. This
+    is a numeric value taken onto the interval (0,1), or it is
+    equal to  \code{"minimum"} or \code{"minimum1se"} for selecting the
+    optimal threshold according to the selection performed in 
+    \code{\link{rsc_cv}}.
+   }
+}
+
+
+\details{
+   The setting \code{threshold = "minimum"} or \code{threshold =
+   "minimum1se"} applies thresholding according to the criteria
+   discussed in the \emph{Details} section in \code{\link{rsc_cv}}.
+   When \code{cv} is obtained using \code{\link{rsc_cv}} with
+   \code{cv.type = "random"}, the default settings for \code{\link{rsc}}
+   implements exactly the RSC estimator proposed in Serra et al.,
+   (2018).
+  
+   Although \code{threshold = "minimum"} is the default choice, in
+   high-dimensional situations \code{threshold = "minimum1se"} usually
+   provides a more parsimonious representation of the correlation
+   structure. Since the underlying RMAD matrix is passed through the
+   \code{cv} input, any other hand-tuned threshold to the RMAD matrix
+   can be applied without significant additional computational
+   costs. The latter can be done setting \code{threshold} to any value
+   onto the (0,1) interval.
+
+   The software is optimized to handle high-dimensional data sets,
+   therefore, the output RSC matrix is packed into a storage efficient
+   sparse format using the \code{"dsCMatrix"} S4 class from the
+   \code{\link{Matrix}} package. The latter is specifically designed for
+   sparse real symmetric matrices.
+}
+
+
+
+\value{
+  Returns a sparse correlaiton matrix of class \code{"dsCMatrix"}
+  (S4 class object) as defined in the \code{\link{Matrix}} package. 
+}
+
+
+
+
+\section{References}{
+   Serra, A., Coretto, P., Fratello, M., and Tagliaferri, R. (2018).
+   Robust and sparsecorrelation matrix estimation for the analysis of
+   high-dimensional genomics data. \emph{Bioinformatics}, 34(4),
+   625-634.  doi:10.1093/bioinformatics/btx642
+}
+
+
+
+\seealso{
+  \code{\link{rsc_cv}}
+}
+
+
+
+
+
+
+
+\examples{
+\donttest{
+## simulate a random sample from a multivariate Cauchy distribution
+## note: example in high-dimension are obtained increasing p
+set.seed(1)
+n   <- 100  # sample size
+p   <- 10   # dimension
+dat <- matrix(rt(n*p, df = 1), nrow = n, ncol = p)
+colnames(dat) <- paste0("Var", 1:p)
+
+   
+## perform 10-fold cross-validation repeated R=10 times
+## note: for multi-core machines experiment with 'ncores'
+set.seed(2)
+a <- rsc_cv(x = dat, R = 10, K = 10)
+a
+
+## obtain the RSC matrix with "minimum" flagged solution 
+b <- rsc(cv = a, threshold = "minimum")
+b
+   
+## obtain the RSC matrix with "minimum1se" flagged solution 
+d <- rsc(cv = a, threshold = "minimum1se")
+d
+
+## since the object 'a' stores the RMAD underlying estimator, we can
+## apply thresholding at any level without re-estimating the RMAD 
+## matrix
+e <- rsc(cv = a, threshold = 0.5)
+e
+}
+}
+
diff --git a/man/rsc_cv.Rd b/man/rsc_cv.Rd
new file mode 100755
index 0000000..3cd2da6
--- /dev/null
+++ b/man/rsc_cv.Rd
@@ -0,0 +1,228 @@
+\name{rsc_cv}
+
+\alias{rsc_cv}
+
+\title{Optimal threshold selection for the RSC estimator}
+
+\description{
+   Perform cross-validation to select an adaptive optimal threshold for
+   the RSC estimator proposed in Serra et al. (2018).
+}
+
+
+\usage{
+  rsc_cv(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) 
+}
+
+
+\arguments{
+  \item{x}{
+    A matrix or a data.frame. Rows of \code{x} correspond to sample units
+    and columns correspond to variables. Categorical variables are not
+    allowed.
+  }
+  \item{cv.type}{
+    A character string indicating the cross-validation algorithm. Possible
+    values are \code{"kfold"} for repeated K-fold cross-validation, and
+    \code{"random"} for random cross-validation (see \emph{Details}).
+  }
+  \item{R}{
+    An integer corresponding to the number of repeated foldings  when
+    \code{cv.type = "kfold"}. When \code{cv.type = "random"} \code{R}
+    defines the number of random splits (see \emph{Details}).
+  }
+  \item{K}{
+    An integer corresponding to the number of \emph{folds} in K-fold
+    cross-validation. Therefore this argument is not relevant when
+    \code{cv.type = "random"}.
+  }
+  \item{threshold}{
+    A sequence of reals taken onto the interval (0,1) defining the
+    threshold values at which the loss is estimated. 
+  }
+  \item{even.correction}{
+    A logical value. It sets the parameter  \code{even.correction} in
+    each of the underlying RMAD computations (see \emph{Details} in
+    \code{\link{rmad}}). 
+  }
+  \item{na.rm}{
+    A logical value, it defines the treatment of missing values in 
+    each of the underlying RMAD computations (see \emph{Details}).
+  }
+  \item{ncores}{
+    An integer value defining the number of cores used for parallel
+    computing. When \code{ncores=NULL} (default), the number \code{r} of
+    available cores is detected, and \code{(r-1)} of them are used
+    (see \emph{Details}).
+  }
+  \item{monitor}{
+    A logical value. If \code{TRUE} progress  messages are
+    printed on screen.
+  }
+}
+
+
+\details{
+  The \code{rsc_cv} function performs cross-validation to estimate the
+  expected Frobenius loss proposed in Bickel and Levina (2008). The
+  original contribution of Bickel and Levina (2008), and its extension
+  in Serra et al. (2018), is based on a random
+  cross-validation  algorithm where the training/test size depends on
+  the sample size \emph{n}. The latter is implemented selecting
+  \code{cv.type = "ramdom"}, and fixing an appropriate number \code{R} of random
+  train/test splits. \code{R} should be as large as possible, but
+  in practice this impacts the computing time strongly for
+  high-dimensional data sets.  
+
+  Although Serra et al. (2018) showed that the random cross-validation
+  of Bickel and Levina (2008) works well for the RSC estimator,
+  subsequent experiments suggested that repeated K-fold cross-validation
+  on average produces better results. Repeated K-fold cross-validation
+  is implemented with the default \code{cv.type = "kfold"}. In this case
+  \code{K} defines the number of \emph{folds}, while \code{R} defines
+  the number of times that the K-fold cross-validation is repeated with
+  \code{R} independent shuffles of the original data. Selecting
+  \code{R=1} and \code{K=10} one performs the standard 10-fold
+  cross-validation. Ten replicates (\code{R=10}) of the K-fold
+  cross-validation are generally sufficient to obtain reasonable
+  estimates of the underlying loss, but for extremely high-dimensional
+  data \code{R} may be varied to speed up calculations. 
+
+  On multi-core hardware the cross-validation is executed in parallel
+  setting \code{ncores}. The parallelism is implemented on the
+  total number of data splits, that is \code{R} for the random
+  cross-validation, and \code{R*K} for the repeated K-fold
+  cross-validation. The software is optimized so that generally the
+  total computing time scales almost linearly with the number of
+  available computer cores (\code{ncores}). 
+
+  For both the random and the K-fold cross-validation it is computed the
+  normalized version of the expected squared Frobenius loss proposed in
+  Bickel and Levina (2008). The normalization is such
+  that the squared Frobenius norm of the identity matrix equals to 1
+  whatever is its dimension.
+
+  Two optimal threshold selection types are reported with flags (see
+  \emph{Value} section below): \code{"minimum"} and
+  \code{"minimum1se"}. The flag \code{"minimum"} denotes the threshold
+  value that minimizes the average loss. The flag \code{"minimum1se"}
+  implements the so called
+  \emph{1-SE rule}: this is the maximum threshold value such that the
+  corresponding average loss is within \emph{1-standard-error} with
+  respect to the threshold that minimizes the average loss
+  (that is the one corresponding to the \code{"minimum"} flag). 
+
+  Since unbiased standard errors for the K-fold cross-validation are
+  impossible to compute (see Bengio and Grandvalet, 2004), when
+  \code{cv.type="kfold"} the reported standard errors have to be
+  considered as a downward biased approximation.
+}
+
+
+
+\value{
+  An S3 object of class \code{'cv_rsc'} with the following components:
+  \item{rmadvec}{
+    A vector containing the lower triangle of the underlying RMAD
+    matrix.
+  }
+  \item{varnames}{
+    A character vector if variable names are available for the input
+    data set \code{x}. Otherwise this is \code{NULL}.
+  }
+  \item{loss}{
+    A data.frame reporting cross-validation estimates. Columns of
+    \code{loss} are as follows:   \code{loss$Threshold} is the threshold value;
+    \code{loss$Average} is averaged loss;   \code{loss$SE} is the standard error
+    for the average loss; \code{loss$Flag="minimum"} denotes the  threshold
+    achieving the minimum average loss; \code{loss$Flag = "*"} denotes threshold
+    values such that the average loss is within \emph{1-standard-error}
+    with respect to the \code{"minimum"} solution.
+  }
+  \item{minimum}{
+    A numeric value. This is the minimum of the average loss. This
+    corresponds to the flag \code{"minimum"} in the  loss component
+    above (see \emph{Details}).
+  }
+  \item{minimum1se}{
+     A numeric value. This is the largest threshold such that the
+     corresponding \code{flag = "*"}. In practice this selects the
+     optimal threshold based on the \emph{1-SE rule} discussed in the
+     \emph{Details} Section above.
+  }
+}
+
+
+
+
+\section{References}{
+   Bengio, Y., and Grandvalet, Y. (2004). No unbiased estimator of the
+   variance of k-fold cross-validation. \emph{Journal of Machine Learning
+   Research}, 5(Sep), 1089-1105.
+
+   Bickel, P. J., and Levina, E. (2008). Covariance regularization by
+   thresholding. The \emph{Annals of Statistics}, 36(6), 2577-2604.
+   doi:10.1214/08-AOS600
+   
+   Serra, A., Coretto, P., Fratello, M., and Tagliaferri, R. (2018).
+   Robust and sparsecorrelation matrix estimation for the analysis of
+   high-dimensional genomics data. \emph{Bioinformatics}, 34(4),
+   625-634.  doi:10.1093/bioinformatics/btx642
+}
+
+
+
+\seealso{
+   \code{rsc}, \code{plot.rsc_cv}
+}
+
+
+
+
+
+
+
+\examples{
+\donttest{
+## simulate a random sample from a multivariate Cauchy distribution
+## note: example in high-dimension are obtained increasing p
+set.seed(1)
+n   <- 100  # sample size
+p   <- 10   # dimension
+dat <- matrix(rt(n*p, df = 1), nrow = n, ncol = p)
+colnames(dat) <- paste0("Var", 1:p)
+
+   
+## perform 10-fold cross-validation repeated R=10 times
+## note: for multi-core machines experiment with 'ncores'
+set.seed(2)
+a <- rsc_cv(x = dat, R = 10, K = 10)
+a
+
+   
+## threshold selection: note that here, knowing the sampling designs,
+## we would like to threshold any correlation larger than zero in
+## absolute value
+## 
+a$minimum        ## "minimum"    flagged solution 
+a$minimum1se     ## "minimum1se" flagged solution
+
+## plot the cross-validation estimates
+plot(a)
+
+## to obtain the RSC matrix we pass 'a' to the rsc() function
+b <- rsc(cv = a, threshold = "minimum")
+b
+
+d <- rsc(cv = a, threshold = "minimum1se")
+d
+
+## since the object 'a' stores the RMAD underlying estimator, we can
+## apply thresholding at any level without re-estimating the RMAD 
+## matrix
+e <- rsc(cv = a, threshold = 0.5)
+e
+}
+}
+
diff --git a/src/Makevars b/src/Makevars
new file mode 100644
index 0000000..8d466ea
--- /dev/null
+++ b/src/Makevars
@@ -0,0 +1,9 @@
+all: $(SHLIB)
+	$(MAKE) $(SHLIB)
+	rm -f *.o *.mod
+
+cormaddp.o: selectiongeneral.o
+	$(FC) -c -fPIC cormaddp.f90 -o cormaddp.o
+
+selectiongeneral.o: selectiongeneral.f95
+	$(FC) -c -fPIC selectiongeneral.f95 -o selectiongeneral.o
diff --git a/src/cormaddp.f90 b/src/cormaddp.f90
new file mode 100644
index 0000000..f9a0f69
--- /dev/null
+++ b/src/cormaddp.f90
@@ -0,0 +1,82 @@
+   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, fresh
+      real(kind=dp), dimension(:), allocatable :: U, V, 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
new file mode 100644
index 0000000..66209f7
--- /dev/null
+++ b/src/init.c
@@ -0,0 +1,22 @@
+#include <stdlib.h>
+//#include <stdio.h>
+#include <R.h>
+#include <Rinternals.h>
+#include <Rmath.h>
+#include <Rdefines.h>
+#include <R_ext/RS.h>
+#include <R_ext/Rdynload.h>
+
+extern void F77_NAME(cormadvecdp)(void *, void *, void *, void *, void *, void *);
+
+static const R_FortranMethodDef FortMethods[] = {
+   {"cormadvecdp", (DL_FUNC) &F77_NAME(cormadvecdp), 6},
+   {NULL, NULL, 0}
+};
+
+void
+R_init_RSC(DllInfo *dll)
+{
+   R_registerRoutines(dll, NULL, NULL, FortMethods, NULL);
+   R_useDynamicSymbols(dll, FALSE);
+}
diff --git a/src/selectiongeneral.f95 b/src/selectiongeneral.f95
new file mode 100644
index 0000000..1078261
--- /dev/null
+++ b/src/selectiongeneral.f95
@@ -0,0 +1,1081 @@
+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 quickselectsp, quickselectscalarsp,quickselectdp, quickselectscalardp
+   end interface qselect
+
+   interface iselect
+      procedure introselectsp, introselectscalarsp, introselectdp, introselectscalardp
+   end interface iselect
+
+   contains
+
+      !!!!!!! SINGLE PRECISION
+   !from here there are the selectionalgo and selectionalgoimpl condensed in one file
+      recursive subroutine quickselectrecursivesp(invector,tovector,k,output)
+         real(kind=sp), dimension(:), allocatable, target, intent(inout) :: invector
+         real(kind=sp),  intent(inout), pointer, contiguous :: tovector(:)
+         integer, intent(inout) :: k
+         real(kind=sp), intent(out) :: output
+         real(kind=sp) :: swapper, pvt
+         integer :: i, r, subst
+         real(kind=sp), 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 quickselectrecursivesp(invector,tovector,k, output)
+            case(:-1)
+               k=k-subst
+               tovector => tovector(subst+1:)
+               call quickselectrecursivesp(invector,tovector,k, output)
+            end select
+         end select
+      end subroutine quickselectrecursivesp
+
+      recursive subroutine introselectrecursivesp(invector,tovector,k,output)
+         real(kind=sp), dimension(:), allocatable, target, intent(inout) :: invector
+         real(kind=sp),  intent(inout), pointer, contiguous :: tovector(:)
+         integer, intent(inout) :: k
+         real(kind=sp), intent(out) :: output
+         real(kind=sp) :: swapper, pvt
+         integer, save :: switch
+         integer :: i, r, subst, r_min, switch_after
+         real(kind=sp), 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 medofmedsp
+                     switch=0 
+                     call medofmedsp(invector,tovector,k, output)
+                  else
+                     call introselectrecursivesp(invector,tovector,k, output)
+                  end if
+               else
+                  switch=0
+                  call introselectrecursivesp(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 medofmedsp
+                     switch=0 
+                     call medofmedsp(invector,tovector,k, output)
+                  else
+                     call introselectrecursivesp(invector,tovector,k, output)
+                  end if
+               else
+                  switch=0
+                  call introselectrecursivesp(invector,tovector,k, output)
+               end if
+            end select
+         end select
+      end subroutine introselectrecursivesp
+
+      recursive subroutine medofmedsp(invector,tovector,k,output)
+         real(kind=sp), dimension(:), allocatable, target, intent(inout) :: invector
+         real(kind=sp),  intent(inout), pointer, contiguous :: tovector(:)
+         real(kind=sp), intent(out) :: output
+         integer, intent(inout) :: k
+
+         real(kind=sp), pointer, contiguous :: subsec(:)
+         real(kind=sp),  dimension(:),allocatable :: meds
+         integer :: i,subst,r
+         real(kind=sp) ::  pvt, swapper
+         real(kind=sp), 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(i/5+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(i/5+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(i/5+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(i/5+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(i/5+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 introselectrecursivesp(invector,tovector,k, output)
+      end subroutine medofmedsp
+
+      function introselectsp(invector,ord,evencorrection)
+         real(kind=sp), intent(in), dimension(:) :: invector
+         logical, optional, intent(in) :: evencorrection
+         integer, intent(in), optional :: ord
+         real(kind=sp) :: introselectsp
+
+         real(kind=sp), dimension(:), allocatable, target :: vector
+         real(kind=sp), 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 introselectrecursivesp(vector,tovector,k, introselectsp)
+
+         if (k1/=0) then
+            block
+               real(kind=sp) :: temp
+               temp=introselectsp
+               tovector=>vector
+               call introselectrecursivesp(vector,tovector,k1, introselectsp)
+               introselectsp=(introselectsp+temp)/2.0
+            end block
+         end if
+      end function introselectsp
+
+      function quickselectsp(invector,ord,evencorrection)
+         real(kind=sp), intent(in), dimension(:) :: invector
+         logical, optional, intent(in) :: evencorrection
+         integer, intent(in), optional :: ord
+         real(kind=sp) :: quickselectsp
+
+         real(kind=sp), dimension(:), allocatable, target :: vector
+         real(kind=sp), 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 quickselectrecursivesp(vector,tovector,k, quickselectsp)
+
+         if (k1/=0) then
+            block
+               real(kind=sp) :: temp
+               temp=quickselectsp
+               tovector=>vector(1:k1+1)
+               call quickselectrecursivesp(vector,tovector,k1, quickselectsp)
+               quickselectsp=(quickselectsp+temp)/2.0
+            end block
+         end if
+      end function quickselectsp
+
+      !deal with scalar input
+      function introselectscalarsp(invector,ord, evencorrection) 
+         real(kind=sp), intent(in) :: invector
+         integer, intent(in), optional :: ord
+         logical, optional, intent(in) :: evencorrection
+         real(kind=sp) :: introselectscalarsp
+
+         introselectscalarsp=invector
+      end function introselectscalarsp
+
+      function quickselectscalarsp(invector,ord, evencorrection) 
+         real(kind=sp), intent(in) :: invector
+         integer, intent(in), optional :: ord
+         logical, optional, intent(in) :: evencorrection
+         real(kind=sp) :: quickselectscalarsp
+
+         quickselectscalarsp=invector
+      end function quickselectscalarsp
+
+      !!!!!!! DOUBLE PRECISION
+   !from here there are the selectionalgo and selectionalgoimpl condensed in one file
+      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, 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(i/5+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(i/5+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(i/5+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(i/5+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(i/5+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
+
+         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
+
+         quickselectscalardp=invector
+      end function quickselectscalardp
+
+
+end module selectionalgo