diff options
Diffstat (limited to 'src/cormaddp.f90')
| -rw-r--r-- | src/cormaddp.f90 | 724 |
1 files changed, 642 insertions, 82 deletions
diff --git a/src/cormaddp.f90 b/src/cormaddp.f90 index f9a0f69..b0b6516 100644 --- a/src/cormaddp.f90 +++ b/src/cormaddp.f90 @@ -1,82 +1,642 @@ - subroutine cormadvecdp(matrix,nrow,ncol,res,ressize,correcteven) - - use selectionalgo - - implicit none - - !note kind dp is defined in selectionalgo - integer, intent(in) :: nrow, ncol,ressize, correcteven !.Fortran R function will not deal with deferred size arrays - real(kind=dp), dimension(nrow,ncol), intent(inout) :: matrix - real(kind=dp), dimension(ressize), intent(out) :: res - - integer :: i, j, n, p - real(kind=dp) :: med,mad, fresh - real(kind=dp), dimension(:), allocatable :: U, V, A, B - - - p=ubound(matrix,2) - n=ubound(matrix,1) - allocate(U(n),V(n)) - - if (correcteven==1) then - do i=1,p-3,3 - U=matrix(:,i) - med=iselect(U,evencorrection=.true.) - matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.)) - U=matrix(:,i+1) - med=iselect(U,evencorrection=.true.) - matrix(:,i+1)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.)) - U=matrix(:,i+2) - med=iselect(U,evencorrection=.true.) - matrix(:,i+2)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.)) - end do - do i=i,p - U=matrix(:,i) - med=iselect(U,evencorrection=.true.) - matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.)) - end do - !unrolled loops - n=1 - do i=1,p-1 - do j=i+1,p - U=matrix(:,i)+matrix(:,j) - V=-matrix(:,i)+matrix(:,j) - mad=(cost*iselect(abs(U-iselect(U,evencorrection=.true.)),evencorrection=.true.))**2 - med=(cost*iselect(abs(V-iselect(V,evencorrection=.true.)),evencorrection=.true.))**2 - res(n)=(mad-med)/(mad+med) - n=n+1 - end do - end do - else - do i=1,p-3,3 - U=matrix(:,i) - med=iselect(U) - matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med))) - U=matrix(:,i+1) - med=iselect(U) - matrix(:,i+1)=(U-med)/(sqrt2*cost*iselect(abs(U-med))) - U=matrix(:,i+2) - med=iselect(U) - matrix(:,i+2)=(U-med)/(sqrt2*cost*iselect(abs(U-med))) - end do - do i=i,p - U=matrix(:,i) - med=iselect(U) - matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med))) - end do - !unrolled loops - n=1 - do i=1,p-1 - do j=i+1,p - U=matrix(:,i)+matrix(:,j) - V=-matrix(:,i)+matrix(:,j) - mad=(cost*iselect(abs(U-iselect(U))))**2 - med=(cost*iselect(abs(V-iselect(V))))**2 - res(n)=(mad-med)/(mad+med) - n=n+1 - end do - end do - end if - - - end subroutine cormadvecdp +module selectionalgo + use, intrinsic :: iso_fortran_env + + implicit none + + real, parameter :: cost=1.4826, sqrt2=sqrt(2.) + integer, parameter :: sp = real32, dp=real64 + + interface qselect + procedure quickselectdp, quickselectscalardp + end interface qselect + + interface iselect + procedure introselectdp, introselectscalardp + end interface iselect + + contains + + !!!!!!! DOUBLE PRECISION + recursive subroutine quickselectrecursivedp(invector,tovector,k,output) + real(kind=dp), dimension(:), allocatable, target, intent(inout) :: invector + real(kind=dp), intent(inout), pointer, contiguous :: tovector(:) + integer, intent(inout) :: k + real(kind=dp), intent(out) :: output + real(kind=dp) :: swapper, pvt + integer :: i, r, subst + real(kind=dp), pointer :: a,b,c + + r=size(tovector) + + select case(r) + case(1) + output=tovector(1) + return + case(2) + select case(k) + case(1) + output=minval(tovector) + case(2) + output=maxval(tovector) + end select + return + case(3) + select case(k) + case(1) + output=minval(tovector) + case(3) + output=maxval(tovector) + case(2) + i=r/2+1 + output=sum(tovector)-minval(tovector)-maxval(tovector) + end select + return + case default + i=r/2+1 + + !pivoting section (pivot of 3) + a=>tovector(1) + b=>tovector(i) + c=>tovector(r) + pvt=a+b+c + swapper=max(a,b,c) + b=swapper + pvt=pvt-swapper + swapper=min(a,b,c) + a=swapper + pvt=pvt-swapper + c=pvt + + !one pass section + subst=1 + b=>tovector(subst) + do i=1,(r-1-3),3 + a => tovector(i) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + a => tovector(i+1) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + a => tovector(i+2) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + end do + do i=i,(r-1) + a => tovector(i) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + end do + tovector(r)=tovector(subst) + tovector(subst)=pvt + + !decide next step + select case(subst-k) + case(0) + output=tovector(k) + case(1:) + tovector => tovector(1:subst-1) + call quickselectrecursivedp(invector,tovector,k, output) + case(:-1) + k=k-subst + tovector => tovector(subst+1:) + call quickselectrecursivedp(invector,tovector,k, output) + end select + end select + end subroutine quickselectrecursivedp + + recursive subroutine introselectrecursivedp(invector,tovector,k,output) + real(kind=dp), dimension(:), allocatable, target, intent(inout) :: invector + real(kind=dp), intent(inout), pointer, contiguous :: tovector(:) + integer, intent(inout) :: k + real(kind=dp), intent(out) :: output + real(kind=dp) :: swapper, pvt + integer, save :: switch + integer :: i, r, subst, r_min, switch_after + real(kind=dp), pointer :: a,b,c + + r=ubound(tovector,1) + r_min = 3000 ! at least this bigger to switch to median of medians + switch_after = 5 ! at least fails this amount of time to reorder + ! one third of the vector to consider med of med + + select case(r) + case(1) + output=tovector(1) + switch=0 + return + case(2) + select case(k) + case(1) + output=minval(tovector) + case(2) + output=maxval(tovector) + end select + switch=0 + return + case(3) + select case(k) + case(1) + output=minval(tovector) + case(3) + output=maxval(tovector) + case(2) + i=r/2+1 + output=sum(tovector)-minval(tovector)-maxval(tovector) + end select + switch=0 + return + case default + i=r/2+1 + + !pivoting section (pivot of 3) + a=>tovector(1) + b=>tovector(i) + c=>tovector(r) + pvt=a+b+c + swapper=max(a,b,c) + b=swapper + pvt=pvt-swapper + swapper=min(a,b,c) + a=swapper + pvt=pvt-swapper + c=pvt + + !one pass section + subst=1 + b=>tovector(subst) + do i=1,(r-1-3),3 + a => tovector(i) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + a => tovector(i+1) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + a => tovector(i+2) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + end do + do i=i,(r-1) + a => tovector(i) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + end do + tovector(r)=tovector(subst) + tovector(subst)=pvt + + !decide next step + select case(subst-k) + case(0) + output=tovector(k) + switch=0 + case(1:) + tovector => tovector(1:subst-1) + if ((subst-1)>r*2/3) then + switch=switch+1 + if (switch==switch_after .and. r>r_min) then + !reset switch and call medofmeddp + switch=0 + call medofmeddp(invector,tovector,k, output) + else + call introselectrecursivedp(invector,tovector,k, output) + end if + else + switch=0 + call introselectrecursivedp(invector,tovector,k, output) + end if + case(:-1) + k=k-subst + tovector => tovector(subst+1:) + if ((r-subst)>r*2/3) then + switch=switch+1 + if (switch==switch_after .and. r>r_min) then + !reset switch and call medofmeddp + switch=0 + call medofmeddp(invector,tovector,k, output) + else + call introselectrecursivedp(invector,tovector,k, output) + end if + else + switch=0 + call introselectrecursivedp(invector,tovector,k, output) + end if + end select + end select + end subroutine introselectrecursivedp + + recursive subroutine medofmeddp(invector,tovector,k,output) + real(kind=dp), dimension(:), allocatable, target, intent(inout) :: invector + real(kind=dp), intent(inout), pointer, contiguous :: tovector(:) + real(kind=dp), intent(out) :: output + integer, intent(inout) :: k + + real(kind=dp), pointer, contiguous :: subsec(:) + real(kind=dp), dimension(:),allocatable :: meds + integer :: i,subst,r + real(kind=dp) :: pvt + !real(kind=dp) :: swapper + !real(kind=dp), pointer :: a,b + + r=ubound(tovector,1) + select case(mod(r,5)) + case(0) + allocate(meds(r/5)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(int(i/5.0)+1)=pvt + end do + case(1) + allocate(meds(r/5+1)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(int(i/5.0)+1)=pvt + end do + meds(r/5+1)=tovector(i) + case(2) + allocate(meds(r/5+1)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(int(i/5.0)+1)=pvt + end do + subsec=>tovector(i:) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + meds(r/5+1)=subsec(2) + case(3) + allocate(meds(r/5+1)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(int(i/5.0)+1)=pvt + end do + subsec=>tovector(i:) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + meds(r/5+1)=pvt + case(4) + allocate(meds(r/5+1)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(int(i/5.0)+1)=pvt + end do + subsec=>tovector(i:) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(r/5+1)=pvt + end select + + !set pivot + pvt=qselect(meds) + do i=1,r-3,3 + if (tovector(i)==pvt) then + tovector(i)=tovector(r) + tovector(r)=pvt + exit + else if (tovector(i+1)==pvt) then + tovector(i+1)=tovector(r) + tovector(r)=pvt + exit + else if (tovector(i+2)==pvt) then + tovector(i+2)=tovector(r) + tovector(r)=pvt + exit + else + cycle + end if + end do + do i=i,r + if (tovector(i)==pvt) then + tovector(i)=tovector(r) + tovector(r)=pvt + exit + else + cycle + end if + end do + + call introselectrecursivedp(invector,tovector,k, output) + end subroutine medofmeddp + + function introselectdp(invector,ord,evencorrection) + real(kind=dp), intent(in), dimension(:) :: invector + logical, optional, intent(in) :: evencorrection + integer, intent(in), optional :: ord + real(kind=dp) :: introselectdp + + real(kind=dp), dimension(:), allocatable, target :: vector + real(kind=dp), pointer, contiguous :: tovector(:) => null() + integer :: k1, k + + k1=size(invector) + allocate(vector(k1)) + vector=invector + tovector => vector + + k1=ubound(invector,1) + if(present(ord)) then + k=ord + else + k=k1/2+1 + end if + + if(present(evencorrection)) then + if (k/=k1/2+1 .or. mod(k1,2)/=0) then + !write(*,'(a)') 'Warning: evencorrection argument ignored.' + k1=0 + else + k1=0 + if (evencorrection) k1=k-1 + end if + else + k1=0 + end if + + call introselectrecursivedp(vector,tovector,k, introselectdp) + + if (k1/=0) then + block + real(kind=dp) :: temp + temp=introselectdp + !tovector=>vector(1:k1+1) + tovector=>vector + call introselectrecursivedp(vector,tovector,k1, introselectdp) + introselectdp=(introselectdp+temp)/2.0 + end block + end if + end function introselectdp + + function quickselectdp(invector,ord,evencorrection) + real(kind=dp), intent(in), dimension(:) :: invector + logical, optional, intent(in) :: evencorrection + integer, intent(in), optional :: ord + real(kind=dp) :: quickselectdp + + real(kind=dp), dimension(:), allocatable, target :: vector + real(kind=dp), pointer, contiguous :: tovector(:) => null() + integer :: k1, k + + k1=size(invector) + allocate(vector(k1)) + vector=invector + tovector => vector + + k1=ubound(invector,1) + if(present(ord)) then + k=ord + else + k=k1/2+1 + end if + + if(present(evencorrection)) then + if (k/=k1/2+1 .or. mod(k1,2)/=0) then + !write(*,'(a)') 'Warning: evencorrection argument ignored.' + k1=0 + else + k1=0 + if (evencorrection) k1=k-1 + end if + else + k1=0 + end if + + call quickselectrecursivedp(vector,tovector,k, quickselectdp) + + if (k1/=0) then + block + real(kind=dp) :: temp + temp=quickselectdp + tovector=>vector(1:k1+1) + call quickselectrecursivedp(vector,tovector,k1, quickselectdp) + quickselectdp=(quickselectdp+temp)/2.0 + end block + end if + end function quickselectdp + + !deal with scalar input + function introselectscalardp(invector,ord, evencorrection) + real(kind=dp), intent(in) :: invector + integer, intent(in), optional :: ord + logical, optional, intent(in) :: evencorrection + real(kind=dp) :: introselectscalardp + + if (present(evencorrection) .or. present(ord)) then + continue + end if + + introselectscalardp=invector + end function introselectscalardp + + function quickselectscalardp(invector,ord, evencorrection) + real(kind=dp), intent(in) :: invector + integer, intent(in), optional :: ord + logical, optional, intent(in) :: evencorrection + real(kind=dp) :: quickselectscalardp + + if (present(evencorrection) .or. present(ord)) then + continue + end if + + quickselectscalardp=invector + end function quickselectscalardp + + +end module selectionalgo + +subroutine cormadvecdp(matrix,nrow,ncol,res,ressize,correcteven) + + use selectionalgo + + implicit none + + !note kind dp is defined in selectionalgo + integer, intent(in) :: nrow, ncol,ressize, correcteven !.Fortran R function will not deal with deferred size arrays + real(kind=dp), dimension(nrow,ncol), intent(inout) :: matrix + real(kind=dp), dimension(ressize), intent(out) :: res + + integer :: i, j, n, p + real(kind=dp) :: med,mad + !real(kind=dp) :: fresh + real(kind=dp), dimension(:), allocatable :: U, V + !real(kind=dp), dimension(:), allocatable :: A, B + + + p=ubound(matrix,2) + n=ubound(matrix,1) + allocate(U(n),V(n)) + + if (correcteven==1) then + do i=1,p-3,3 + U=matrix(:,i) + med=iselect(U,evencorrection=.true.) + matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.)) + U=matrix(:,i+1) + med=iselect(U,evencorrection=.true.) + matrix(:,i+1)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.)) + U=matrix(:,i+2) + med=iselect(U,evencorrection=.true.) + matrix(:,i+2)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.)) + end do + do i=i,p + U=matrix(:,i) + med=iselect(U,evencorrection=.true.) + matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.)) + end do + !unrolled loops + n=1 + do i=1,p-1 + do j=i+1,p + U=matrix(:,i)+matrix(:,j) + V=-matrix(:,i)+matrix(:,j) + mad=(cost*iselect(abs(U-iselect(U,evencorrection=.true.)),evencorrection=.true.))**2 + med=(cost*iselect(abs(V-iselect(V,evencorrection=.true.)),evencorrection=.true.))**2 + res(n)=(mad-med)/(mad+med) + n=n+1 + end do + end do + else + do i=1,p-3,3 + U=matrix(:,i) + med=iselect(U) + matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med))) + U=matrix(:,i+1) + med=iselect(U) + matrix(:,i+1)=(U-med)/(sqrt2*cost*iselect(abs(U-med))) + U=matrix(:,i+2) + med=iselect(U) + matrix(:,i+2)=(U-med)/(sqrt2*cost*iselect(abs(U-med))) + end do + do i=i,p + U=matrix(:,i) + med=iselect(U) + matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med))) + end do + !unrolled loops + n=1 + do i=1,p-1 + do j=i+1,p + U=matrix(:,i)+matrix(:,j) + V=-matrix(:,i)+matrix(:,j) + mad=(cost*iselect(abs(U-iselect(U))))**2 + med=(cost*iselect(abs(V-iselect(V))))**2 + res(n)=(mad-med)/(mad+med) + n=n+1 + end do + end do + end if +end subroutine cormadvecdp |