version 1.2

10 files changed, 670 insertions, 1193 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 36f47f6..3f611ee 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 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>.
+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",
@@ -20,8 +20,8 @@ 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
+Version: 1.2
+Date: 2020-07-24
+Packaged: 2020-07-24 07:47:23 UTC; pietro
 Repository: CRAN
-Date/Publication: 2020-07-04 10:50:03 UTC
+Date/Publication: 2020-07-24 14:00:07 UTC
diff --git a/MD5 b/MD5
index f14927a..7aaaa51 100644
--- a/MD5
+++ b/MD5
@@ -1,5 +1,6 @@
-80c587bf8be31b19549b965aa44f1d50 *DESCRIPTION
+7a9c09c1b58001a6ef17502cccb6b9ee *DESCRIPTION
 1f64d30a870f2eae4855e83010aaaa4e *NAMESPACE
+1eae13396706a60fc95f5ff4b751765b *NEWS
 0b87a473118069181446fb007472c042 *R/check_inputs.R
 1df1ecba75175b1665069fc88419cf02 *R/cv_loss.R
 1ada6bc39f15ea58f34a5d46173141fa *R/plot_print_methods.R
@@ -7,12 +8,10 @@
 f26fdd51a2378e500e6ed2eaeda8721a *R/rsc.R
 556a9a10b85bfad178103f3cc5fdc0ff *R/rsc_cv.R
 607039efb30fc06f123290e3dcbce1d4 *R/zzz.R
-ee184b1708ad4ef4735070065fe6737c *inst/CITATION
-fffcd92fdbda643ea59f12298cc04287 *man/plot.cv_rsc.Rd
+02a06053793e54ce413a378a05ccc3ce *inst/CITATION
+1e6af804927098cfa6c235278682578c *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
+8cc5c9373c7aa8835e9818d0de274ece *man/rsc.Rd
+63d3e58e9bdc0da9969feb6ce69d005f *man/rsc_cv.Rd
+f45a14eed865277afe99ebcb1bed66f8 *src/cormaddp.f90
 9ded1c66a2fd469e79d7d1d25d0bd604 *src/init.c
-62f9244cb578a6720b664233da503e32 *src/selectiongeneral.f95
diff --git a/NEWS b/NEWS
new file mode 100755
index 0000000..d7efdcf
--- /dev/null
+++ b/NEWS
@@ -0,0 +1,12 @@
+## 24/07/2020 at 09:44:39 (CEST)

+* v1.2 uploaded

+* fixed compiler issue on Solaris

+

+

+## 17/07/2020 at 09:01:44 (CEST)

+* v1.1  uploaded

+* fixed an issue related to the Fortran compiler on Solaris

+

+

+## 02/07/2020 at 09:00:01 (CEST)

+* v1.0 uploaded 

diff --git a/inst/CITATION b/inst/CITATION
index 8690fa8..2ee61ef 100644
--- a/inst/CITATION
+++ b/inst/CITATION
@@ -14,7 +14,7 @@ citEntry(entry="Manual",
          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",
+                note, "url: https://CRAN.R-project.org/package=RSC",
                 sep='')
         )
 
diff --git a/man/plot.cv_rsc.Rd b/man/plot.cv_rsc.Rd
index 0d7b2c2..030a53f 100755
--- a/man/plot.cv_rsc.Rd
+++ b/man/plot.cv_rsc.Rd
@@ -72,7 +72,7 @@ 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 <- rsc_cv(x = dat, R = 10, K = 10, ncores = 1)
 a
 
 ## plot the cross-validation estimates
diff --git a/man/rsc.Rd b/man/rsc.Rd
index cd252a6..6371af3 100755
--- a/man/rsc.Rd
+++ b/man/rsc.Rd
@@ -84,7 +84,6 @@
 
 
 \examples{
-\donttest{
 ## simulate a random sample from a multivariate Cauchy distribution
 ## note: example in high-dimension are obtained increasing p
 set.seed(1)
@@ -97,7 +96,7 @@ 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 <- rsc_cv(x = dat, R = 10, K = 10, ncores = 1)
 a
 
 ## obtain the RSC matrix with "minimum" flagged solution 
@@ -114,5 +113,4 @@ d
 e <- rsc(cv = a, threshold = 0.5)
 e
 }
-}
 
diff --git a/man/rsc_cv.Rd b/man/rsc_cv.Rd
index 3cd2da6..44459ec 100755
--- a/man/rsc_cv.Rd
+++ b/man/rsc_cv.Rd
@@ -184,7 +184,6 @@
 
 
 \examples{
-\donttest{
 ## simulate a random sample from a multivariate Cauchy distribution
 ## note: example in high-dimension are obtained increasing p
 set.seed(1)
@@ -197,7 +196,7 @@ 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 <- rsc_cv(x = dat, R = 10, K = 10, ncores = 1)
 a
 
    
@@ -224,5 +223,4 @@ d
 e <- rsc(cv = a, threshold = 0.5)
 e
 }
-}
 
diff --git a/src/Makevars b/src/Makevars
deleted file mode 100644
index 8d466ea..0000000
--- a/src/Makevars
+++ /dev/null
@@ -1,9 +0,0 @@
-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
index f9a0f69..b0b6516 100644
--- a/src/cormaddp.f90
+++ b/src/cormaddp.f90
@@ -1,82 +1,642 @@
-   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
+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/selectiongeneral.f95 b/src/selectiongeneral.f95
deleted file mode 100644
index 1078261..0000000
--- a/src/selectiongeneral.f95
+++ /dev/null
@@ -1,1081 +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 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