From 511e3ca9e5235e018f772693907d9ec10002b02a Mon Sep 17 00:00:00 2001
From: Luca Coraggio <luca.coraggio@unina.it>
Date: Sat, 4 Jul 2020 09:50:03 +0000
Subject: version 1.0

---
 src/Makevars             |    9 +
 src/cormaddp.f90         |   82 ++++
 src/init.c               |   22 +
 src/selectiongeneral.f95 | 1081 ++++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 1194 insertions(+)
 create mode 100644 src/Makevars
 create mode 100644 src/cormaddp.f90
 create mode 100644 src/init.c
 create mode 100644 src/selectiongeneral.f95

(limited to 'src')

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
-- 
cgit v1.2.3