version 1.2

1 files changed, 0 insertions, 1081 deletions

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