diff --git a/radixalgorithm/mergesort.f90 b/radixalgorithm/mergesort.f90 index f8935074e69d289abe20766c78bec944495f10b9..0bbe950f3fb862f23673fb9f53988a962ed07f28 100644 --- a/radixalgorithm/mergesort.f90 +++ b/radixalgorithm/mergesort.f90 @@ -1,66 +1,67 @@ - subroutine merge_sort(n, a, indices) -! This routine takes an input array 'a' of dimension (n) -! and an input array 'indices' of dimension (n) -! and conducts a merge sort on array 'a' while saving the original order -! to array 'indices' - -! The first pass sorts each pair of consecutive elements [a(i) vs. a(i+1)]. -! The second pass sorts adjacent sequences of length 2 [a(i), a(i+1) vs. a(i+2), a(i+3)]; the sequences are already ordered, due to step 1. -! The third pass sorts adjacent sequences of length 4 [a(i),...,a(i+3) vs. a(i+4),...,a(i+7)]; the sequences are already ordered, from step 2. -! This continues for L1 passes, where 2**(L1-1) < n < 2**L1. Note that n does not need to be a power of 2. - - implicit none - integer :: i, imax, i0, j, jmax, k, L, L1, m, n - real, intent(inout) :: a(n) - integer, intent(out) :: indices(n) - real, allocatable :: b(:) - allocate (b(n)) - L1 = 1 - m = 1 - do while (m < n) ! Determine L1 so that 2**(L1-1) < n < 2**L1 - m = m + m - L1 = L1 + 1 - end do - L1 = L1 - 1 - m = 1 - do L=1, L1 - k=1 - do i0=1, n-m+1, m+m - i=i0 - j=i+m - imax=j - jmax=j+m - if(imax > n) imax = n + 1 - if(jmax > n) jmax = n + 1 - do while(i < imax .and. j < jmax) - if(a(i) < a(j)) then - b(k) = a(i) - indices(k) = i - i = i + 1 - else - b(k) = a(j) - indices(k) = j - j = j + 1 - end if - k = k + 1 - end do - do while(i < imax) - b(k) = a(i) - indices(k) = i - i = i + 1 - k = k + 1 - end do - do while(j < jmax) - b(k) = a(j) - indices(k) = j - j = j + 1 - k = k + 1 - end do - end do - m = m + m - a = b - end do - if(allocated(b)) deallocate(b) - return - end subroutine merge_sort - + subroutine merge_sort(n, a, indices) +! This routine takes an input array 'a' of dimension (n) +! and an input array 'indices' of dimension (n) +! and conducts a merge sort on array 'a' while saving the original order +! to array 'indices' + +! The first pass sorts each pair of consecutive elements [a(i) vs. a(i+1)]. +! The second pass sorts adjacent sequences of length 2 [a(i), a(i+1) vs. a(i+2), a(i+3)]; the sequences are already ordered, due to step 1. +! The third pass sorts adjacent sequences of length 4 [a(i),...,a(i+3) vs. a(i+4),...,a(i+7)]; the sequences are already ordered, from step 2. +! This continues for L1 passes, where 2**(L1-1) < n < 2**L1. Note that n does not need to be a power of 2. + + implicit none + !DEC\$ ATTRIBUTES DLLEXPORT::merge_sort + integer :: i, imax, i0, j, jmax, k, L, L1, m, n + real, intent(inout) :: a(n) + integer, intent(out) :: indices(n) + real, allocatable :: b(:) + allocate (b(n)) + L1 = 1 + m = 1 + do while (m < n) ! Determine L1 so that 2**(L1-1) < n < 2**L1 + m = m + m + L1 = L1 + 1 + end do + L1 = L1 - 1 + m = 1 + do L=1, L1 + k=1 + do i0=1, n-m+1, m+m + i=i0 + j=i+m + imax=j + jmax=j+m + if(imax > n) imax = n + 1 + if(jmax > n) jmax = n + 1 + do while(i < imax .and. j < jmax) + if(a(i) < a(j)) then + b(k) = a(i) + indices(k) = i + i = i + 1 + else + b(k) = a(j) + indices(k) = j + j = j + 1 + end if + k = k + 1 + end do + do while(i < imax) + b(k) = a(i) + indices(k) = i + i = i + 1 + k = k + 1 + end do + do while(j < jmax) + b(k) = a(j) + indices(k) = j + j = j + 1 + k = k + 1 + end do + end do + m = m + m + a = b + end do + if(allocated(b)) deallocate(b) + return + end subroutine merge_sort +