diff options
Diffstat (limited to 'src/selectiongeneral.f95')
| -rw-r--r-- | src/selectiongeneral.f95 | 1081 |
1 files changed, 1081 insertions, 0 deletions
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 |