From 4925226beec2aeb5c7d58f38a4a19ab284e9c67b Mon Sep 17 00:00:00 2001 From: "0.382" <18322825326@163.com> Date: Sun, 7 May 2023 16:47:56 +0800 Subject: [PATCH 1/7] Add radix_sort --- example/sorting/CMakeLists.txt | 1 + example/sorting/example_radix_sort.f90 | 63 ++++ src/CMakeLists.txt | 1 + src/stdlib_sorting.fypp | 91 +++++ src/stdlib_sorting_radix_sort.f90 | 454 +++++++++++++++++++++++++ test/sorting/test_sorting.f90 | 97 ++++++ 6 files changed, 707 insertions(+) create mode 100644 example/sorting/example_radix_sort.f90 create mode 100644 src/stdlib_sorting_radix_sort.f90 diff --git a/example/sorting/CMakeLists.txt b/example/sorting/CMakeLists.txt index 8a08998d5..4800cc376 100644 --- a/example/sorting/CMakeLists.txt +++ b/example/sorting/CMakeLists.txt @@ -1,2 +1,3 @@ ADD_EXAMPLE(ord_sort) ADD_EXAMPLE(sort) +ADD_EXAMPLE(radix_sort) \ No newline at end of file diff --git a/example/sorting/example_radix_sort.f90 b/example/sorting/example_radix_sort.f90 new file mode 100644 index 000000000..ccdae2c18 --- /dev/null +++ b/example/sorting/example_radix_sort.f90 @@ -0,0 +1,63 @@ +program example_radix_sort + use iso_fortran_env + use iso_c_binding + use ieee_arithmetic + use stdlib_sorting, only: radix_sort, ord_sort + implicit none + integer(int8), allocatable :: arri8(:) + integer(int16), allocatable :: arri16(:) + real(real64), allocatable, target :: arrf64(:), x + + real(real32), dimension(:, :), allocatable :: arr1, arr2 + real(real32), dimension(:), allocatable :: work + integer :: i, test_size, repeat + real :: start, stop + + arri8 = [-128, 127, 0, -1, 1] + call radix_sort(arri8) + print *, arri8 + + arri16 = [-32767, 32767, 0, 0, -3, 2, -3] + call radix_sort(arri16, reverse=.true.) + print *, arri16 + + allocate (arrf64(10)) + x = 0.0 ! divide zero will arise compile error + arrf64(1) = 1.0/x + arrf64(2) = 0.0 + arrf64(3) = 0.0/x + arrf64(4) = -1.0/x + arrf64(5) = -0.0 + arrf64(6) = 1.0 + arrf64(7) = -1.0 + arrf64(8) = 3.45 + arrf64(9) = -3.14 + arrf64(10) = 3.44 + call radix_sort(arrf64) + print *, arrf64 + ! In my computer, it gives + ! nan, -inf, -3.14, -1.0, -0.0, 0.0, 1.0, 3.44, 3.45, inf + ! but the position of nan is undefined, the position of `±inf`, `±0.0` is as expected + + test_size = 100000 + repeat = 100 + allocate (arr1(test_size, repeat)) + allocate (arr2(test_size, repeat)) + call random_number(arr1) + arr1 = arr1 - 0.5 + arr2(:, :) = arr1(:, :) + allocate (work(test_size)) + call cpu_time(start) + do i = 1, repeat + call ord_sort(arr1(:, i), work) + end do + call cpu_time(stop) + print *, "ord_sort time = ", (stop - start) + call cpu_time(start) + do i = 1, repeat + call radix_sort(arr2(:, i), work) + end do + call cpu_time(stop) + print *, "radix_sort time = ", (stop - start) + print *, all(arr1 == arr2) +end program example_radix_sort diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ceb1bd2b9..8a6fe66cc 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -76,6 +76,7 @@ set(SRC stdlib_hashmap_chaining.f90 stdlib_hashmap_open.f90 stdlib_logger.f90 + stdlib_sorting_radix_sort.f90 stdlib_system.F90 stdlib_specialfunctions.f90 stdlib_specialfunctions_legendre.f90 diff --git a/src/stdlib_sorting.fypp b/src/stdlib_sorting.fypp index 3ef5ca9ae..9d9316fc2 100644 --- a/src/stdlib_sorting.fypp +++ b/src/stdlib_sorting.fypp @@ -380,6 +380,97 @@ module stdlib_sorting end interface ord_sort + public radix_sort +!! Version: experimental +!! +!! The generic subroutine implementing the LSD radix sort algorithm to return +!! an input array with its elements sorted in order of (non-)decreasing +!! value. Its use has the syntax: +!! +!! call radix_sort( array[, work, reverse] ) +!! +!! with the arguments: +!! +!! * array: the rank 1 array to be sorted. It is an `intent(inout)` +!! argument of any of the types `integer(int8)`, `integer(int16)`, +!! `integer(int32)`, `integer(int64)`, `real(real32)`, `real(real64)`. +!! If both the type of `array` is real and at least one of the +!! elements is a `NaN`, then the ordering of the result is undefined. +!! Otherwise it is defined to be the original elements in +!! non-decreasing order. Especially, -0.0 is lesser than 0.0. +!! +!! * work (optional): shall be a rank 1 array of the same type as +!! `array`, and shall have at least `size(array)` elements. It is an +!! `intent(inout)` argument to be used as buffer. Its value on return is +!! undefined. If it is not present, `radix_sort` will allocate a +!! buffer for use, and deallocate it bufore return. If you do several +!! similar `radix_sort`, reuse the `work` array is a good parctice. +!! This argument is not present for `int8_radix_sort` because it use +!! counting sort, so no buffer is needed. +!! +!! * `reverse` (optional): shall be a scalar of type default logical. It +!! is an `intent(in)` argument. If present with a value of `.true.` then +!! `array` will be sorted in order of non-increasing values in stable +!! order. Otherwise index will sort `array` in order of non-decreasing +!! values in stable order. +!! +!!#### Example +!! +!!```fortran +!! ... +!! ! Read random data from a file +!! call read_file( 'dummy_file', array ) +!! ! Sort the random data +!! call radix_sort( array ) +!! ! Process the sorted data +!! call array_search( array, values ) +!! ... +!!``` + + interface radix_sort +!! Version: experimental +!! +!! The generic subroutine interface implementing the LSD radix sort algorithm, +!! see https://en.wikipedia.org/wiki/Radix_sort for more details. +!! It is always O(N) in sorting random data, but need a O(N) buffer. +!! + + pure module subroutine int8_radix_sort(array, reverse) + integer(kind=int8), dimension(:), intent(inout) :: array + logical, intent(in), optional :: reverse + end subroutine int8_radix_sort + + pure module subroutine int16_radix_sort(array, work, reverse) + integer(kind=int16), dimension(:), intent(inout) :: array + integer(kind=int16), dimension(:), intent(inout), target, optional :: work + logical, intent(in), optional :: reverse + end subroutine int16_radix_sort + + pure module subroutine int32_radix_sort(array, work, reverse) + integer(kind=int32), dimension(:), intent(inout) :: array + integer(kind=int32), dimension(:), intent(inout), target, optional :: work + logical, intent(in), optional :: reverse + end subroutine int32_radix_sort + + pure module subroutine int64_radix_sort(array, work, reverse) + integer(kind=int64), dimension(:), intent(inout) :: array + integer(kind=int64), dimension(:), intent(inout), target, optional :: work + logical, intent(in), optional :: reverse + end subroutine int64_radix_sort + + module subroutine sp_radix_sort(array, work, reverse) + real(kind=sp), dimension(:), intent(inout), target :: array + real(kind=sp), dimension(:), intent(inout), target, optional :: work + logical, intent(in), optional :: reverse + end subroutine sp_radix_sort + + module subroutine dp_radix_sort(array, work, reverse) + real(kind=dp), dimension(:), intent(inout), target :: array + real(kind=dp), dimension(:), intent(inout), target, optional :: work + logical, intent(in), optional :: reverse + end subroutine dp_radix_sort + end interface radix_sort + interface sort !! Version: experimental !! diff --git a/src/stdlib_sorting_radix_sort.f90 b/src/stdlib_sorting_radix_sort.f90 new file mode 100644 index 000000000..9f55626e9 --- /dev/null +++ b/src/stdlib_sorting_radix_sort.f90 @@ -0,0 +1,454 @@ +submodule(stdlib_sorting) stdlib_sorting_radix_sort + implicit none + + integer, parameter :: radix_bits = 8 + integer, parameter :: radix_mask = 255 + integer(kind=int16), parameter :: radix_bits_i16 = 8_int16 + integer(kind=int16), parameter :: radix_mask_i16 = 255_int16 + integer(kind=int32), parameter :: radix_bits_i32 = 8_int32 + integer(kind=int32), parameter :: radix_mask_i32 = 255_int32 + integer(kind=int64), parameter :: radix_bits_i64 = 8_int64 + integer(kind=int64), parameter :: radix_mask_i64 = 255_int64 + +contains + pure subroutine radix_sort_u8_helper(N, arr) + integer(kind=int64), intent(in) :: N + integer(kind=int8), dimension(N), intent(inout) :: arr + integer(kind=int64) :: i + integer :: bin_idx + integer(kind=int64), dimension(-128:127) :: counts + counts(:) = 0 + do i = 1, N + bin_idx = arr(i) + counts(bin_idx) = counts(bin_idx) + 1 + end do + i = 1 + do bin_idx = -128, 127 + do while (counts(bin_idx) > 0) + arr(i) = int(bin_idx, kind=int8) + i = i + 1 + counts(bin_idx) = counts(bin_idx) - 1 + end do + end do + end subroutine + + pure subroutine radix_sort_u16_helper(N, arr, buf) + integer(kind=int64), intent(in) :: N + integer(kind=int16), dimension(N), intent(inout) :: arr + integer(kind=int16), dimension(N), intent(inout) :: buf + integer(kind=int64) :: i + integer :: b, b0, b1 + integer(kind=int64), dimension(0:radix_mask) :: c0, c1 + c0(:) = 0 + c1(:) = 0 + do i = 1, N + b0 = iand(arr(i), radix_mask_i16) + b1 = ishft(arr(i), -radix_bits_i16) + c0(b0) = c0(b0) + 1 + c1(b1) = c1(b1) + 1 + end do + do b = 1, radix_mask + c0(b) = c0(b) + c0(b - 1) + c1(b) = c1(b) + c1(b - 1) + end do + do i = N, 1, -1 + b0 = iand(arr(i), radix_mask_i16) + buf(c0(b0)) = arr(i) + c0(b0) = c0(b0) - 1 + end do + do i = N, 1, -1 + b1 = ishft(buf(i), -radix_bits_i16) + arr(c1(b1)) = buf(i) + c1(b1) = c1(b1) - 1 + end do + end subroutine + + pure subroutine radix_sort_u32_helper(N, arr, buf) + integer(kind=int64), intent(in) :: N + integer(kind=int32), dimension(N), intent(inout) :: arr + integer(kind=int32), dimension(N), intent(inout) :: buf + integer(kind=int64) :: i + integer :: b, b0, b1, b2, b3 + integer(kind=int64), dimension(0:radix_mask) :: c0, c1, c2, c3 + c0(:) = 0 + c1(:) = 0 + c2(:) = 0 + c3(:) = 0 + do i = 1, N + b0 = iand(arr(i), radix_mask_i32) + b1 = iand(ishft(arr(i), -radix_bits_i32), radix_mask_i32) + b2 = iand(ishft(arr(i), -2*radix_bits_i32), radix_mask_i32) + b3 = ishft(arr(i), -3*radix_bits_i32) + c0(b0) = c0(b0) + 1 + c1(b1) = c1(b1) + 1 + c2(b2) = c2(b2) + 1 + c3(b3) = c3(b3) + 1 + end do + do b = 1, radix_mask + c0(b) = c0(b) + c0(b - 1) + c1(b) = c1(b) + c1(b - 1) + c2(b) = c2(b) + c2(b - 1) + c3(b) = c3(b) + c3(b - 1) + end do + do i = N, 1, -1 + b0 = iand(arr(i), radix_mask_i32) + buf(c0(b0)) = arr(i) + c0(b0) = c0(b0) - 1 + end do + do i = N, 1, -1 + b1 = iand(ishft(buf(i), -radix_bits_i32), radix_mask_i32) + arr(c1(b1)) = buf(i) + c1(b1) = c1(b1) - 1 + end do + do i = N, 1, -1 + b2 = iand(ishft(arr(i), -2*radix_bits_i32), radix_mask_i32) + buf(c2(b2)) = arr(i) + c2(b2) = c2(b2) - 1 + end do + do i = N, 1, -1 + b3 = ishft(buf(i), -3*radix_bits_i32) + arr(c3(b3)) = buf(i) + c3(b3) = c3(b3) - 1 + end do + end subroutine radix_sort_u32_helper + + pure subroutine radix_sort_u64_helper(N, arr, buffer) + integer(kind=int64), intent(in) :: N + integer(kind=int64), dimension(N), intent(inout) :: arr + integer(kind=int64), dimension(N), intent(inout) :: buffer + integer(kind=int64) :: i + integer(kind=int64) :: b, b0, b1, b2, b3, b4, b5, b6, b7 + integer(kind=int64), dimension(0:radix_mask) :: c0, c1, c2, c3, c4, c5, c6, c7 + c0(:) = 0 + c1(:) = 0 + c2(:) = 0 + c3(:) = 0 + c4(:) = 0 + c5(:) = 0 + c6(:) = 0 + c7(:) = 0 + do i = 1, N + b0 = iand(arr(i), radix_mask_i64) + b1 = iand(ishft(arr(i), -radix_bits_i64), radix_mask_i64) + b2 = iand(ishft(arr(i), -2*radix_bits_i64), radix_mask_i64) + b3 = iand(ishft(arr(i), -3*radix_bits_i64), radix_mask_i64) + b4 = iand(ishft(arr(i), -4*radix_bits_i64), radix_mask_i64) + b5 = iand(ishft(arr(i), -5*radix_bits_i64), radix_mask_i64) + b6 = iand(ishft(arr(i), -6*radix_bits_i64), radix_mask_i64) + b7 = ishft(arr(i), -7*radix_bits_i64) + c0(b0) = c0(b0) + 1 + c1(b1) = c1(b1) + 1 + c2(b2) = c2(b2) + 1 + c3(b3) = c3(b3) + 1 + c4(b4) = c4(b4) + 1 + c5(b5) = c5(b5) + 1 + c6(b6) = c6(b6) + 1 + c7(b7) = c7(b7) + 1 + end do + do b = 1, radix_mask + c0(b) = c0(b) + c0(b - 1) + c1(b) = c1(b) + c1(b - 1) + c2(b) = c2(b) + c2(b - 1) + c3(b) = c3(b) + c3(b - 1) + c4(b) = c4(b) + c4(b - 1) + c5(b) = c5(b) + c5(b - 1) + c6(b) = c6(b) + c6(b - 1) + c7(b) = c7(b) + c7(b - 1) + end do + do i = N, 1, -1 + b0 = iand(arr(i), radix_mask_i64) + buffer(c0(b0)) = arr(i) + c0(b0) = c0(b0) - 1 + end do + do i = N, 1, -1 + b1 = iand(ishft(buffer(i), -radix_bits_i64), radix_mask_i64) + arr(c1(b1)) = buffer(i) + c1(b1) = c1(b1) - 1 + end do + do i = N, 1, -1 + b2 = iand(ishft(arr(i), -2*radix_bits_i64), radix_mask_i64) + buffer(c2(b2)) = arr(i) + c2(b2) = c2(b2) - 1 + end do + do i = N, 1, -1 + b3 = iand(ishft(buffer(i), -3*radix_bits_i64), radix_mask_i64) + arr(c3(b3)) = buffer(i) + c3(b3) = c3(b3) - 1 + end do + do i = N, 1, -1 + b4 = iand(ishft(arr(i), -4*radix_bits_i64), radix_mask_i64) + buffer(c4(b4)) = arr(i) + c4(b4) = c4(b4) - 1 + end do + do i = N, 1, -1 + b5 = iand(ishft(buffer(i), -5*radix_bits_i64), radix_mask_i64) + arr(c5(b5)) = buffer(i) + c5(b5) = c5(b5) - 1 + end do + do i = N, 1, -1 + b6 = iand(ishft(arr(i), -6*radix_bits_i64), radix_mask_i64) + buffer(c6(b6)) = arr(i) + c6(b6) = c6(b6) - 1 + end do + do i = N, 1, -1 + b7 = ishft(buffer(i), -7*radix_bits_i64) + arr(c7(b7)) = buffer(i) + c7(b7) = c7(b7) - 1 + end do + end subroutine radix_sort_u64_helper + + pure module subroutine int8_radix_sort(array, reverse) + integer(kind=int8), dimension(:), intent(inout) :: array + logical, intent(in), optional :: reverse + integer(kind=int8) :: buffer_item + integer(kind=int64) :: i, N + N = size(array, kind=int64) + call radix_sort_u8_helper(N, array) + if (optval(reverse, .false.)) then + do i = 1, N/2 + buffer_item = array(i) + array(i) = array(N - i + 1) + array(N - i + 1) = buffer_item + end do + end if + end subroutine int8_radix_sort + + pure module subroutine int16_radix_sort(array, work, reverse) + integer(kind=int16), dimension(:), intent(inout) :: array + integer(kind=int16), dimension(:), intent(inout), target, optional :: work + logical, intent(in), optional :: reverse + integer(kind=int64) :: i, N, start, middle, end + integer(kind=int16), dimension(:), pointer :: buffer + integer(kind=int16) :: item + logical :: use_internal_buffer + N = size(array, kind=int64) + if (present(work)) then + if (size(work, kind=int64) < N) then + error stop "int16_radix_sort: work array is too small." + end if + use_internal_buffer = .false. + buffer => work + else + use_internal_buffer = .true. + allocate (buffer(N)) + end if + call radix_sort_u16_helper(N, array, buffer) + if (array(1) >= 0 .and. array(N) < 0) then + start = 1 + end = N + middle = (1 + N)/2 + do while (.true.) + if (array(middle) >= 0) then + start = middle + 1 + else + end = middle + end if + middle = (start + end)/2 + if (start == end) exit + end do + buffer(1:(N - middle + 1)) = array(middle:N) + buffer(N - middle + 2:N) = array(1:middle - 1) + array(:) = buffer(:) + end if + if (optval(reverse, .false.)) then + do i = 1, N/2 + item = array(i) + array(i) = array(N - i + 1) + array(N - i + 1) = item + end do + end if + if (use_internal_buffer) then + deallocate (buffer) + end if + end subroutine int16_radix_sort + + pure module subroutine int32_radix_sort(array, work, reverse) + integer(kind=int32), dimension(:), intent(inout) :: array + integer(kind=int32), dimension(:), intent(inout), target, optional :: work + logical, intent(in), optional :: reverse + integer(kind=int64) :: i, N, start, middle, end + integer(kind=int32), dimension(:), pointer :: buffer + integer(kind=int32) :: item + logical :: use_internal_buffer + N = size(array, kind=int64) + if (present(work)) then + if (size(work, kind=int64) < N) then + error stop "int32_radix_sort: work array is too small." + end if + use_internal_buffer = .false. + buffer => work + else + use_internal_buffer = .true. + allocate (buffer(N)) + end if + call radix_sort_u32_helper(N, array, buffer) + if (array(1) >= 0 .and. array(N) < 0) then + start = 1 + end = N + middle = (1 + N)/2 + do while (.true.) + if (array(middle) >= 0) then + start = middle + 1 + else + end = middle + end if + middle = (start + end)/2 + if (start == end) exit + end do + buffer(1:(N - middle + 1)) = array(middle:N) + buffer(N - middle + 2:N) = array(1:middle - 1) + array(:) = buffer(:) + end if + if (optval(reverse, .false.)) then + do i = 1, N/2 + item = array(i) + array(i) = array(N - i + 1) + array(N - i + 1) = item + end do + end if + if (use_internal_buffer) then + deallocate (buffer) + end if + end subroutine int32_radix_sort + + module subroutine sp_radix_sort(array, work, reverse) + use iso_c_binding + real(kind=sp), dimension(:), intent(inout), target :: array + real(kind=sp), dimension(:), intent(inout), target, optional :: work + logical, intent(in), optional :: reverse + integer(kind=int64) :: i, N, pos, rev_pos + integer(kind=int32), dimension(:), pointer :: arri32 + integer(kind=int32), dimension(:), pointer :: buffer + real(kind=sp) :: item + logical :: use_internal_buffer + N = size(array, kind=int64) + if (present(work)) then + if (size(work, kind=int64) < N) then + error stop "sp_radix_sort: work array is too small." + end if + use_internal_buffer = .false. + call c_f_pointer(c_loc(work), buffer, [N]) + else + use_internal_buffer = .true. + allocate (buffer(N)) + end if + call c_f_pointer(c_loc(array), arri32, [N]) + call radix_sort_u32_helper(N, arri32, buffer) + if (arri32(1) >= 0 .and. arri32(N) < 0) then + pos = 1 + rev_pos = N + do while (arri32(rev_pos) < 0) + buffer(pos) = arri32(rev_pos) + pos = pos + 1 + rev_pos = rev_pos - 1 + end do + buffer(pos:N) = arri32(1:rev_pos) + arri32(:) = buffer(:) + end if + if (optval(reverse, .false.)) then + do i = 1, N/2 + item = array(i) + array(i) = array(N - i + 1) + array(N - i + 1) = item + end do + end if + if (use_internal_buffer) then + deallocate (buffer) + end if + end subroutine sp_radix_sort + + pure module subroutine int64_radix_sort(array, work, reverse) + integer(kind=int64), dimension(:), intent(inout) :: array + integer(kind=int64), dimension(:), intent(inout), target, optional :: work + logical, intent(in), optional :: reverse + integer(kind=int64) :: i, N, start, middle, end + integer(kind=int64), dimension(:), pointer :: buffer + integer(kind=int64) :: item + logical :: use_internal_buffer + N = size(array, kind=int64) + if (present(work)) then + if (size(work, kind=int64) < N) then + error stop "int64_radix_sort: work array is too small." + end if + use_internal_buffer = .false. + buffer => work + else + use_internal_buffer = .true. + allocate (buffer(N)) + end if + call radix_sort_u64_helper(N, array, buffer) + if (array(1) >= 0 .and. array(N) < 0) then + start = 1 + end = N + middle = (1 + N)/2 + do while (.true.) + if (array(middle) >= 0) then + start = middle + 1 + else + end = middle + end if + middle = (start + end)/2 + if (start == end) exit + end do + buffer(1:(N - middle + 1)) = array(middle:N) + buffer(N - middle + 2:N) = array(1:middle - 1) + array(:) = buffer(:) + end if + if (optval(reverse, .false.)) then + do i = 1, N/2 + item = array(i) + array(i) = array(N - i + 1) + array(N - i + 1) = item + end do + end if + if (use_internal_buffer) then + deallocate (buffer) + end if + end subroutine int64_radix_sort + + module subroutine dp_radix_sort(array, work, reverse) + use iso_c_binding + real(kind=dp), dimension(:), intent(inout), target :: array + real(kind=dp), dimension(:), intent(inout), target, optional :: work + logical, intent(in), optional :: reverse + integer(kind=int64) :: i, N, pos, rev_pos + integer(kind=int64), dimension(:), pointer :: arri64 + integer(kind=int64), dimension(:), pointer :: buffer + real(kind=dp) :: item + logical :: use_internal_buffer + N = size(array, kind=int64) + if (present(work)) then + if (size(work, kind=int64) < N) then + error stop "sp_radix_sort: work array is too small." + end if + use_internal_buffer = .false. + call c_f_pointer(c_loc(work), buffer, [N]) + else + use_internal_buffer = .true. + allocate (buffer(N)) + end if + call c_f_pointer(c_loc(array), arri64, [N]) + call radix_sort_u64_helper(N, arri64, buffer) + if (arri64(1) >= 0 .and. arri64(N) < 0) then + pos = 1 + rev_pos = N + do while (arri64(rev_pos) < 0) + buffer(pos) = arri64(rev_pos) + pos = pos + 1 + rev_pos = rev_pos - 1 + end do + buffer(pos:N) = arri64(1:rev_pos) + arri64(:) = buffer(:) + end if + if (optval(reverse, .false.)) then + do i = 1, N/2 + item = array(i) + array(i) = array(N - i + 1) + array(N - i + 1) = item + end do + end if + if (use_internal_buffer) then + deallocate (buffer) + end if + end subroutine dp_radix_sort +end submodule stdlib_sorting_radix_sort diff --git a/test/sorting/test_sorting.f90 b/test/sorting/test_sorting.f90 index 23b730e3b..9538f0076 100644 --- a/test/sorting/test_sorting.f90 +++ b/test/sorting/test_sorting.f90 @@ -64,6 +64,7 @@ subroutine collect_sorting(testsuite) new_unittest('int_ord_sorts', test_int_ord_sorts), & new_unittest('char_ord_sorts', test_char_ord_sorts), & new_unittest('string_ord_sorts', test_string_ord_sorts), & + new_unittest('int_radix_sorts', test_int_radix_sorts), & new_unittest('int_sorts', test_int_sorts), & new_unittest('char_sorts', test_char_sorts), & new_unittest('string_sorts', test_string_sorts), & @@ -446,6 +447,102 @@ subroutine test_string_ord_sort( a, a_name, ltest ) end subroutine test_string_ord_sort + subroutine test_int_radix_sorts(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + integer(int64) :: i + integer, allocatable :: d1(:) + logical :: ltest + + call test_int_radix_sort( blocks, "Blocks", ltest ) + call check(error, ltest) + if (allocated(error)) return + + call test_int_radix_sort( decrease, "Decreasing", ltest ) + call check(error, ltest) + if (allocated(error)) return + + call test_int_radix_sort( identical, "Identical", ltest ) + call check(error, ltest) + if (allocated(error)) return + + call test_int_radix_sort( increase, "Increasing", ltest ) + call check(error, ltest) + if (allocated(error)) return + + call test_int_radix_sort( rand1, "Random dense", ltest ) + call check(error, ltest) + if (allocated(error)) return + + call test_int_radix_sort( rand2, "Random order", ltest ) + call check(error, ltest) + if (allocated(error)) return + + call test_int_radix_sort( rand0, "Random sparse", ltest ) + call check(error, ltest) + if (allocated(error)) return + + call test_int_radix_sort( rand3, "Random 3", ltest ) + call check(error, ltest) + if (allocated(error)) return + + call test_int_radix_sort( rand10, "Random 10", ltest ) + call check(error, ltest) + if (allocated(error)) return + + !triggered an issue in insertion + d1 = [10, 2, -3, -4, 6, -6, 7, -8, 9, 0, 1, 20] + call sort( d1 ) + call verify_sort( d1, ltest, i ) + call check(error, ltest) + + end subroutine test_int_radix_sorts + + subroutine test_int_radix_sort( a, a_name, ltest ) + integer(int32), intent(in) :: a(:) + character(*), intent(in) :: a_name + logical, intent(out) :: ltest + + integer(int64) :: t0, t1, tdiff + real(dp) :: rate + integer(int64) :: i + logical :: valid + + ltest = .true. + + tdiff = 0 + do i = 1, repeat + dummy = a + call system_clock( t0, rate ) + call radix_sort( dummy ) + call system_clock( t1, rate ) + tdiff = tdiff + t1 - t0 + end do + tdiff = tdiff/repeat + + call verify_sort( dummy, valid, i ) + ltest = (ltest .and. valid) + if ( .not. valid ) then + write( *, * ) "SORT did not sort " // a_name // "." + write(*,*) 'i = ', i + write(*,'(a12, 2i7)') 'dummy(i-1:i) = ', dummy(i-1:i) + end if + write( lun, '("| Integer |", 1x, i7, 2x, "|", 1x, a15, " |", ' // & + 'a12, " |", F10.6, " |" )' ) & + test_size, a_name, "Sort", tdiff/rate + + ! reverse + dummy = a + call radix_sort( dummy, reverse = .true.) + call verify_reverse_sort(dummy, valid, i) + ltest = (ltest .and. valid) + if ( .not. valid ) then + write( *, * ) "reverse SORT did not sort " // a_name // "." + write(*,*) 'i = ', i + write(*,'(a12, 2i7)') 'dummy(i-1:i) = ', dummy(i-1:i) + end if + + end subroutine test_int_radix_sort subroutine test_int_sorts(error) !> Error handling From 80adf2ba96207bc1501786c66257a0f8650568e0 Mon Sep 17 00:00:00 2001 From: 0382 <18322825326@163.com> Date: Sun, 7 May 2023 19:59:16 +0800 Subject: [PATCH 2/7] Apply suggestions from code review Co-authored-by: Jeremie Vandenplas --- example/sorting/example_radix_sort.f90 | 11 +---------- test/sorting/test_sorting.f90 | 2 +- 2 files changed, 2 insertions(+), 11 deletions(-) diff --git a/example/sorting/example_radix_sort.f90 b/example/sorting/example_radix_sort.f90 index ccdae2c18..abccc7b30 100644 --- a/example/sorting/example_radix_sort.f90 +++ b/example/sorting/example_radix_sort.f90 @@ -23,16 +23,7 @@ program example_radix_sort allocate (arrf64(10)) x = 0.0 ! divide zero will arise compile error - arrf64(1) = 1.0/x - arrf64(2) = 0.0 - arrf64(3) = 0.0/x - arrf64(4) = -1.0/x - arrf64(5) = -0.0 - arrf64(6) = 1.0 - arrf64(7) = -1.0 - arrf64(8) = 3.45 - arrf64(9) = -3.14 - arrf64(10) = 3.44 + arrf64 = [1.0/x, 0.0, 0.0/x, -1.0/x, -0.0, 1.0, -1.0, 3.45, -3.14, 3.44] call radix_sort(arrf64) print *, arrf64 ! In my computer, it gives diff --git a/test/sorting/test_sorting.f90 b/test/sorting/test_sorting.f90 index 9538f0076..9de0ccb91 100644 --- a/test/sorting/test_sorting.f90 +++ b/test/sorting/test_sorting.f90 @@ -529,7 +529,7 @@ subroutine test_int_radix_sort( a, a_name, ltest ) end if write( lun, '("| Integer |", 1x, i7, 2x, "|", 1x, a15, " |", ' // & 'a12, " |", F10.6, " |" )' ) & - test_size, a_name, "Sort", tdiff/rate + test_size, a_name, "Radix_sort", tdiff/rate ! reverse dummy = a From a4c8037b1340d2dc62bc5707a2f0afade7278827 Mon Sep 17 00:00:00 2001 From: "0.382" <18322825326@163.com> Date: Sun, 7 May 2023 21:52:50 +0800 Subject: [PATCH 3/7] add test_real_radix_sorts --- example/sorting/example_radix_sort.f90 | 39 ++-------- src/stdlib_sorting.fypp | 93 ++++++++++++------------ src/stdlib_sorting_radix_sort.f90 | 76 +++++++++++--------- test/sorting/test_sorting.f90 | 98 +++++++++++++++++++++++++- 4 files changed, 189 insertions(+), 117 deletions(-) diff --git a/example/sorting/example_radix_sort.f90 b/example/sorting/example_radix_sort.f90 index abccc7b30..02204a716 100644 --- a/example/sorting/example_radix_sort.f90 +++ b/example/sorting/example_radix_sort.f90 @@ -1,17 +1,10 @@ program example_radix_sort - use iso_fortran_env - use iso_c_binding - use ieee_arithmetic - use stdlib_sorting, only: radix_sort, ord_sort + use iso_fortran_env, only: int8, int16, dp => real64 + use stdlib_sorting, only: radix_sort implicit none integer(int8), allocatable :: arri8(:) integer(int16), allocatable :: arri16(:) - real(real64), allocatable, target :: arrf64(:), x - - real(real32), dimension(:, :), allocatable :: arr1, arr2 - real(real32), dimension(:), allocatable :: work - integer :: i, test_size, repeat - real :: start, stop + real(dp), allocatable, target :: arrf64(:), x arri8 = [-128, 127, 0, -1, 1] call radix_sort(arri8) @@ -22,33 +15,11 @@ program example_radix_sort print *, arri16 allocate (arrf64(10)) - x = 0.0 ! divide zero will arise compile error - arrf64 = [1.0/x, 0.0, 0.0/x, -1.0/x, -0.0, 1.0, -1.0, 3.45, -3.14, 3.44] + x = 0.0_dp ! divide zero will arise compile error + arrf64 = [1.0_dp/x, 0.0_dp, 0.0_dp/x, -1.0_dp/x, -0.0_dp, 1.0_dp, -1.0_dp, 3.45_dp, -3.14_dp, 3.44_dp] call radix_sort(arrf64) print *, arrf64 ! In my computer, it gives ! nan, -inf, -3.14, -1.0, -0.0, 0.0, 1.0, 3.44, 3.45, inf ! but the position of nan is undefined, the position of `±inf`, `±0.0` is as expected - - test_size = 100000 - repeat = 100 - allocate (arr1(test_size, repeat)) - allocate (arr2(test_size, repeat)) - call random_number(arr1) - arr1 = arr1 - 0.5 - arr2(:, :) = arr1(:, :) - allocate (work(test_size)) - call cpu_time(start) - do i = 1, repeat - call ord_sort(arr1(:, i), work) - end do - call cpu_time(stop) - print *, "ord_sort time = ", (stop - start) - call cpu_time(start) - do i = 1, repeat - call radix_sort(arr2(:, i), work) - end do - call cpu_time(stop) - print *, "radix_sort time = ", (stop - start) - print *, all(arr1 == arr2) end program example_radix_sort diff --git a/src/stdlib_sorting.fypp b/src/stdlib_sorting.fypp index 9d9316fc2..c0c034c51 100644 --- a/src/stdlib_sorting.fypp +++ b/src/stdlib_sorting.fypp @@ -236,6 +236,51 @@ module stdlib_sorting !! ! Process the sorted data !! call array_search( array, values ) !! ... +!!``` + + public radix_sort +!! Version: experimental +!! +!! The generic subroutine implementing the LSD radix sort algorithm to return +!! an input array with its elements sorted in order of (non-)decreasing +!! value. Its use has the syntax: +!! +!! call radix_sort( array[, work, reverse] ) +!! +!! with the arguments: +!! +!! * array: the rank 1 array to be sorted. It is an `intent(inout)` +!! argument of any of the types `integer(int8)`, `integer(int16)`, +!! `integer(int32)`, `integer(int64)`, `real(real32)`, `real(real64)`. +!! If both the type of `array` is real and at least one of the +!! elements is a `NaN`, then the ordering of the result is undefined. +!! Otherwise it is defined to be the original elements in +!! non-decreasing order. Especially, -0.0 is lesser than 0.0. +!! +!! * work (optional): shall be a rank 1 array of the same type as +!! `array`, and shall have at least `size(array)` elements. It is an +!! `intent(inout)` argument to be used as buffer. Its value on return is +!! undefined. If it is not present, `radix_sort` will allocate a +!! buffer for use, and deallocate it before return. If you do several +!! similar `radix_sort`s, reusing the `work` array is a good parctice. +!! This argument is not present for `int8_radix_sort` because it use +!! counting sort, so no buffer is needed. +!! +!! * `reverse` (optional): shall be a scalar of type default logical. It +!! is an `intent(in)` argument. If present with a value of `.true.` then +!! `array` will be sorted in order of non-increasing values in stable +!! order. Otherwise index will sort `array` in order of non-decreasing +!! values in stable order. +!! +!!#### Example +!! +!!```fortran +!! ... +!! ! Read random data from a file +!! call read_file( 'dummy_file', array ) +!! ! Sort the random data +!! call radix_sort( array ) +!! ... !!``` public sort_index @@ -379,54 +424,6 @@ module stdlib_sorting #:endfor end interface ord_sort - - public radix_sort -!! Version: experimental -!! -!! The generic subroutine implementing the LSD radix sort algorithm to return -!! an input array with its elements sorted in order of (non-)decreasing -!! value. Its use has the syntax: -!! -!! call radix_sort( array[, work, reverse] ) -!! -!! with the arguments: -!! -!! * array: the rank 1 array to be sorted. It is an `intent(inout)` -!! argument of any of the types `integer(int8)`, `integer(int16)`, -!! `integer(int32)`, `integer(int64)`, `real(real32)`, `real(real64)`. -!! If both the type of `array` is real and at least one of the -!! elements is a `NaN`, then the ordering of the result is undefined. -!! Otherwise it is defined to be the original elements in -!! non-decreasing order. Especially, -0.0 is lesser than 0.0. -!! -!! * work (optional): shall be a rank 1 array of the same type as -!! `array`, and shall have at least `size(array)` elements. It is an -!! `intent(inout)` argument to be used as buffer. Its value on return is -!! undefined. If it is not present, `radix_sort` will allocate a -!! buffer for use, and deallocate it bufore return. If you do several -!! similar `radix_sort`, reuse the `work` array is a good parctice. -!! This argument is not present for `int8_radix_sort` because it use -!! counting sort, so no buffer is needed. -!! -!! * `reverse` (optional): shall be a scalar of type default logical. It -!! is an `intent(in)` argument. If present with a value of `.true.` then -!! `array` will be sorted in order of non-increasing values in stable -!! order. Otherwise index will sort `array` in order of non-decreasing -!! values in stable order. -!! -!!#### Example -!! -!!```fortran -!! ... -!! ! Read random data from a file -!! call read_file( 'dummy_file', array ) -!! ! Sort the random data -!! call radix_sort( array ) -!! ! Process the sorted data -!! call array_search( array, values ) -!! ... -!!``` - interface radix_sort !! Version: experimental !! diff --git a/src/stdlib_sorting_radix_sort.f90 b/src/stdlib_sorting_radix_sort.f90 index 9f55626e9..a21e87185 100644 --- a/src/stdlib_sorting_radix_sort.f90 +++ b/src/stdlib_sorting_radix_sort.f90 @@ -11,12 +11,13 @@ integer(kind=int64), parameter :: radix_mask_i64 = 255_int64 contains +! For int8, radix sort becomes counting sort, so buffer is not needed pure subroutine radix_sort_u8_helper(N, arr) - integer(kind=int64), intent(in) :: N + integer(kind=int_size), intent(in) :: N integer(kind=int8), dimension(N), intent(inout) :: arr - integer(kind=int64) :: i + integer(kind=int_size) :: i integer :: bin_idx - integer(kind=int64), dimension(-128:127) :: counts + integer(kind=int_size), dimension(-128:127) :: counts counts(:) = 0 do i = 1, N bin_idx = arr(i) @@ -33,12 +34,12 @@ pure subroutine radix_sort_u8_helper(N, arr) end subroutine pure subroutine radix_sort_u16_helper(N, arr, buf) - integer(kind=int64), intent(in) :: N + integer(kind=int_size), intent(in) :: N integer(kind=int16), dimension(N), intent(inout) :: arr integer(kind=int16), dimension(N), intent(inout) :: buf - integer(kind=int64) :: i + integer(kind=int_size) :: i integer :: b, b0, b1 - integer(kind=int64), dimension(0:radix_mask) :: c0, c1 + integer(kind=int_size), dimension(0:radix_mask) :: c0, c1 c0(:) = 0 c1(:) = 0 do i = 1, N @@ -64,12 +65,12 @@ pure subroutine radix_sort_u16_helper(N, arr, buf) end subroutine pure subroutine radix_sort_u32_helper(N, arr, buf) - integer(kind=int64), intent(in) :: N + integer(kind=int_size), intent(in) :: N integer(kind=int32), dimension(N), intent(inout) :: arr integer(kind=int32), dimension(N), intent(inout) :: buf - integer(kind=int64) :: i + integer(kind=int_size) :: i integer :: b, b0, b1, b2, b3 - integer(kind=int64), dimension(0:radix_mask) :: c0, c1, c2, c3 + integer(kind=int_size), dimension(0:radix_mask) :: c0, c1, c2, c3 c0(:) = 0 c1(:) = 0 c2(:) = 0 @@ -113,12 +114,12 @@ pure subroutine radix_sort_u32_helper(N, arr, buf) end subroutine radix_sort_u32_helper pure subroutine radix_sort_u64_helper(N, arr, buffer) - integer(kind=int64), intent(in) :: N + integer(kind=int_size), intent(in) :: N integer(kind=int64), dimension(N), intent(inout) :: arr integer(kind=int64), dimension(N), intent(inout) :: buffer - integer(kind=int64) :: i + integer(kind=int_size) :: i integer(kind=int64) :: b, b0, b1, b2, b3, b4, b5, b6, b7 - integer(kind=int64), dimension(0:radix_mask) :: c0, c1, c2, c3, c4, c5, c6, c7 + integer(kind=int_size), dimension(0:radix_mask) :: c0, c1, c2, c3, c4, c5, c6, c7 c0(:) = 0 c1(:) = 0 c2(:) = 0 @@ -200,15 +201,15 @@ end subroutine radix_sort_u64_helper pure module subroutine int8_radix_sort(array, reverse) integer(kind=int8), dimension(:), intent(inout) :: array logical, intent(in), optional :: reverse - integer(kind=int8) :: buffer_item - integer(kind=int64) :: i, N - N = size(array, kind=int64) + integer(kind=int8) :: item + integer(kind=int_size) :: i, N + N = size(array, kind=int_size) call radix_sort_u8_helper(N, array) if (optval(reverse, .false.)) then do i = 1, N/2 - buffer_item = array(i) + item = array(i) array(i) = array(N - i + 1) - array(N - i + 1) = buffer_item + array(N - i + 1) = item end do end if end subroutine int8_radix_sort @@ -217,13 +218,13 @@ pure module subroutine int16_radix_sort(array, work, reverse) integer(kind=int16), dimension(:), intent(inout) :: array integer(kind=int16), dimension(:), intent(inout), target, optional :: work logical, intent(in), optional :: reverse - integer(kind=int64) :: i, N, start, middle, end + integer(kind=int_size) :: i, N, start, middle, end integer(kind=int16), dimension(:), pointer :: buffer integer(kind=int16) :: item logical :: use_internal_buffer - N = size(array, kind=int64) + N = size(array, kind=int_size) if (present(work)) then - if (size(work, kind=int64) < N) then + if (size(work, kind=int_size) < N) then error stop "int16_radix_sort: work array is too small." end if use_internal_buffer = .false. @@ -233,6 +234,9 @@ pure module subroutine int16_radix_sort(array, work, reverse) allocate (buffer(N)) end if call radix_sort_u16_helper(N, array, buffer) +! After calling `radix_sort_u_helper. The array is sorted as unsigned integers. +! In the view of signed arrsy, the negative numbers are sorted but in the tail of the array. +! Use binary search to find the boundary, and move then to correct position. if (array(1) >= 0 .and. array(N) < 0) then start = 1 end = N @@ -266,13 +270,13 @@ pure module subroutine int32_radix_sort(array, work, reverse) integer(kind=int32), dimension(:), intent(inout) :: array integer(kind=int32), dimension(:), intent(inout), target, optional :: work logical, intent(in), optional :: reverse - integer(kind=int64) :: i, N, start, middle, end + integer(kind=int_size) :: i, N, start, middle, end integer(kind=int32), dimension(:), pointer :: buffer integer(kind=int32) :: item logical :: use_internal_buffer - N = size(array, kind=int64) + N = size(array, kind=int_size) if (present(work)) then - if (size(work, kind=int64) < N) then + if (size(work, kind=int_size) < N) then error stop "int32_radix_sort: work array is too small." end if use_internal_buffer = .false. @@ -316,14 +320,14 @@ module subroutine sp_radix_sort(array, work, reverse) real(kind=sp), dimension(:), intent(inout), target :: array real(kind=sp), dimension(:), intent(inout), target, optional :: work logical, intent(in), optional :: reverse - integer(kind=int64) :: i, N, pos, rev_pos + integer(kind=int_size) :: i, N, pos, rev_pos integer(kind=int32), dimension(:), pointer :: arri32 integer(kind=int32), dimension(:), pointer :: buffer real(kind=sp) :: item logical :: use_internal_buffer - N = size(array, kind=int64) + N = size(array, kind=int_size) if (present(work)) then - if (size(work, kind=int64) < N) then + if (size(work, kind=int_size) < N) then error stop "sp_radix_sort: work array is too small." end if use_internal_buffer = .false. @@ -334,6 +338,14 @@ module subroutine sp_radix_sort(array, work, reverse) end if call c_f_pointer(c_loc(array), arri32, [N]) call radix_sort_u32_helper(N, arri32, buffer) +! After calling `radix_sort_u_helper. The array is sorted as unsigned integers. +! The positive real number is sorted, guaranteed by IEEE-754 standard. +! But the negative real number is sorted in a reversed order, and also in the tail of array. +! Remark that -0.0 is the minimum nagative integer, so using radix sort, -0.0 is naturally lesser than 0.0. +! In IEEE-754 standard, the bit representation of `Inf` is greater than all positive real numbers, +! and the `-Inf` is lesser than all negative real numbers. So the order of `Inf, -Inf` is ok. +! The bit representation of `NaN` may be positive or negative integers in different machine, +! thus if the array contains `NaN`, the result is undefined. if (arri32(1) >= 0 .and. arri32(N) < 0) then pos = 1 rev_pos = N @@ -361,13 +373,13 @@ pure module subroutine int64_radix_sort(array, work, reverse) integer(kind=int64), dimension(:), intent(inout) :: array integer(kind=int64), dimension(:), intent(inout), target, optional :: work logical, intent(in), optional :: reverse - integer(kind=int64) :: i, N, start, middle, end + integer(kind=int_size) :: i, N, start, middle, end integer(kind=int64), dimension(:), pointer :: buffer integer(kind=int64) :: item logical :: use_internal_buffer - N = size(array, kind=int64) + N = size(array, kind=int_size) if (present(work)) then - if (size(work, kind=int64) < N) then + if (size(work, kind=int_size) < N) then error stop "int64_radix_sort: work array is too small." end if use_internal_buffer = .false. @@ -411,14 +423,14 @@ module subroutine dp_radix_sort(array, work, reverse) real(kind=dp), dimension(:), intent(inout), target :: array real(kind=dp), dimension(:), intent(inout), target, optional :: work logical, intent(in), optional :: reverse - integer(kind=int64) :: i, N, pos, rev_pos + integer(kind=int_size) :: i, N, pos, rev_pos integer(kind=int64), dimension(:), pointer :: arri64 integer(kind=int64), dimension(:), pointer :: buffer real(kind=dp) :: item logical :: use_internal_buffer - N = size(array, kind=int64) + N = size(array, kind=int_size) if (present(work)) then - if (size(work, kind=int64) < N) then + if (size(work, kind=int_size) < N) then error stop "sp_radix_sort: work array is too small." end if use_internal_buffer = .false. diff --git a/test/sorting/test_sorting.f90 b/test/sorting/test_sorting.f90 index 9de0ccb91..03b4ac5d1 100644 --- a/test/sorting/test_sorting.f90 +++ b/test/sorting/test_sorting.f90 @@ -27,6 +27,7 @@ module test_sorting rand2(0:test_size-1), & rand3(0:test_size-1), & rand10(0:test_size-1) + real(sp) :: rand_r32(0:test_size-1) character(len=4) :: & char_decrease(0:char_size-1), & char_increase(0:char_size-1), & @@ -37,6 +38,7 @@ module test_sorting string_rand(0:string_size-1) integer(int32) :: dummy(0:test_size-1) + real(sp) :: real_dummy(0:test_size-1) character(len=4) :: char_dummy(0:char_size-1) type(string_type) :: string_dummy(0:string_size-1) integer(int_size) :: index(0:max(test_size, char_size, string_size)-1) @@ -65,6 +67,7 @@ subroutine collect_sorting(testsuite) new_unittest('char_ord_sorts', test_char_ord_sorts), & new_unittest('string_ord_sorts', test_string_ord_sorts), & new_unittest('int_radix_sorts', test_int_radix_sorts), & + new_unittest('real_radix_sorts', test_real_radix_sorts), & new_unittest('int_sorts', test_int_sorts), & new_unittest('char_sorts', test_char_sorts), & new_unittest('string_sorts', test_string_sorts), & @@ -117,6 +120,9 @@ subroutine initialize_tests() rand10(i) = int( floor( arand * test_size ), kind=int32 ) end do + call random_number(rand_r32) + rand_r32 = rand_r32 - 0.5 ! to test both positive and negative numbers + count = 0 do i=0, char_set_size-1 do j=0, char_set_size-1 @@ -523,13 +529,13 @@ subroutine test_int_radix_sort( a, a_name, ltest ) call verify_sort( dummy, valid, i ) ltest = (ltest .and. valid) if ( .not. valid ) then - write( *, * ) "SORT did not sort " // a_name // "." + write( *, * ) "RADIX_SORT did not sort " // a_name // "." write(*,*) 'i = ', i write(*,'(a12, 2i7)') 'dummy(i-1:i) = ', dummy(i-1:i) end if write( lun, '("| Integer |", 1x, i7, 2x, "|", 1x, a15, " |", ' // & 'a12, " |", F10.6, " |" )' ) & - test_size, a_name, "Radix_sort", tdiff/rate + test_size, a_name, "Radix_Sort", tdiff/rate ! reverse dummy = a @@ -537,13 +543,69 @@ subroutine test_int_radix_sort( a, a_name, ltest ) call verify_reverse_sort(dummy, valid, i) ltest = (ltest .and. valid) if ( .not. valid ) then - write( *, * ) "reverse SORT did not sort " // a_name // "." + write( *, * ) "reverse RADIX_SORT did not sort " // a_name // "." write(*,*) 'i = ', i write(*,'(a12, 2i7)') 'dummy(i-1:i) = ', dummy(i-1:i) end if end subroutine test_int_radix_sort + subroutine test_real_radix_sort( a, a_name, ltest ) + real(sp), intent(in) :: a(:) + character(*), intent(in) :: a_name + logical, intent(out) :: ltest + + integer(int64) :: t0, t1, tdiff + real(dp) :: rate + integer(int64) :: i + logical :: valid + + ltest = .true. + + tdiff = 0 + do i = 1, repeat + real_dummy = a + call system_clock( t0, rate ) + call radix_sort( real_dummy ) + call system_clock( t1, rate ) + tdiff = tdiff + t1 - t0 + end do + tdiff = tdiff/repeat + + call verify_real_sort( real_dummy, valid, i ) + ltest = (ltest .and. valid) + if ( .not. valid ) then + write( *, * ) "RADIX_SORT did not sort " // a_name // "." + write(*,*) 'i = ', i + write(*,'(a12, 2f12.5)') 'real_dummy(i-1:i) = ', real_dummy(i-1:i) + end if + write( lun, '("| Integer |", 1x, i7, 2x, "|", 1x, a15, " |", ' // & + 'a12, " |", F10.6, " |" )' ) & + test_size, a_name, "Radix_Sort", tdiff/rate + + ! reverse + real_dummy = a + call radix_sort( real_dummy, reverse = .true.) + call verify_real_reverse_sort(real_dummy, valid, i) + ltest = (ltest .and. valid) + if ( .not. valid ) then + write( *, * ) "reverse RADIX_SORT did not sort " // a_name // "." + write(*,*) 'i = ', i + write(*,'(a12, 2f12.5)') 'real_dummy(i-1:i) = ', real_dummy(i-1:i) + end if + + end subroutine test_real_radix_sort + + subroutine test_real_radix_sorts(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + logical :: ltest + + call test_real_radix_sort( rand_r32, "rand-real32", ltest ) + call check(error, ltest) + if (allocated(error)) return + end subroutine test_real_radix_sorts + subroutine test_int_sorts(error) !> Error handling type(error_type), allocatable, intent(out) :: error @@ -1000,7 +1062,21 @@ subroutine verify_sort( a, valid, i ) end subroutine verify_sort + subroutine verify_real_sort( a, valid, i ) + real(sp), intent(in) :: a(0:) + logical, intent(out) :: valid + integer(int64), intent(out) :: i + integer(int64) :: n + + n = size( a, kind=int64 ) + valid = .false. + do i=1, n-1 + if ( a(i-1) > a(i) ) return + end do + valid = .true. + + end subroutine verify_real_sort subroutine verify_string_sort( a, valid, i ) type(string_type), intent(in) :: a(0:) @@ -1066,6 +1142,22 @@ subroutine verify_reverse_sort( a, valid, i ) end subroutine verify_reverse_sort + subroutine verify_real_reverse_sort( a, valid, i ) + real(sp), intent(in) :: a(0:) + logical, intent(out) :: valid + integer(int64), intent(out) :: i + + integer(int64) :: n + + n = size( a, kind=int64 ) + valid = .false. + do i=1, n-1 + if ( a(i-1) < a(i) ) return + end do + valid = .true. + + end subroutine verify_real_reverse_sort + subroutine verify_string_reverse_sort( a, valid, i ) type(string_type), intent(in) :: a(0:) logical, intent(out) :: valid From d49178a29995a68e68a7e6450a3f40c08dcb9dd3 Mon Sep 17 00:00:00 2001 From: 0382 <18322825326@163.com> Date: Wed, 17 May 2023 20:18:27 +0800 Subject: [PATCH 4/7] Apply suggestions from code review Co-authored-by: Jeremie Vandenplas --- example/sorting/example_radix_sort.f90 | 7 ++++--- src/stdlib_sorting_radix_sort.f90 | 4 ++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/example/sorting/example_radix_sort.f90 b/example/sorting/example_radix_sort.f90 index 02204a716..21b355c8f 100644 --- a/example/sorting/example_radix_sort.f90 +++ b/example/sorting/example_radix_sort.f90 @@ -4,7 +4,8 @@ program example_radix_sort implicit none integer(int8), allocatable :: arri8(:) integer(int16), allocatable :: arri16(:) - real(dp), allocatable, target :: arrf64(:), x + real(dp) :: x + real(dp), allocatable :: arrf64(:) arri8 = [-128, 127, 0, -1, 1] call radix_sort(arri8) @@ -19,7 +20,7 @@ program example_radix_sort arrf64 = [1.0_dp/x, 0.0_dp, 0.0_dp/x, -1.0_dp/x, -0.0_dp, 1.0_dp, -1.0_dp, 3.45_dp, -3.14_dp, 3.44_dp] call radix_sort(arrf64) print *, arrf64 - ! In my computer, it gives + ! Expect output: ! nan, -inf, -3.14, -1.0, -0.0, 0.0, 1.0, 3.44, 3.45, inf - ! but the position of nan is undefined, the position of `±inf`, `±0.0` is as expected + ! Note: the position of nan is undefined end program example_radix_sort diff --git a/src/stdlib_sorting_radix_sort.f90 b/src/stdlib_sorting_radix_sort.f90 index a21e87185..9b23562f7 100644 --- a/src/stdlib_sorting_radix_sort.f90 +++ b/src/stdlib_sorting_radix_sort.f90 @@ -235,8 +235,8 @@ pure module subroutine int16_radix_sort(array, work, reverse) end if call radix_sort_u16_helper(N, array, buffer) ! After calling `radix_sort_u_helper. The array is sorted as unsigned integers. -! In the view of signed arrsy, the negative numbers are sorted but in the tail of the array. -! Use binary search to find the boundary, and move then to correct position. +! In the view of signed array, the negative numbers are sorted but in the tail of the array. +! Use binary search to find the boundary, and move them to correct position. if (array(1) >= 0 .and. array(N) < 0) then start = 1 end = N From 4931ecb78677c2ffcb7c90f2f905bc86833ef941 Mon Sep 17 00:00:00 2001 From: "0.382" <18322825326@163.com> Date: Thu, 18 May 2023 16:48:05 +0800 Subject: [PATCH 5/7] add specs file for radix_sort --- doc/specs/stdlib_sorting.md | 288 ++++++++++++++++++++++------------ src/stdlib_sorting.fypp | 1 + test/sorting/test_sorting.f90 | 2 +- 3 files changed, 187 insertions(+), 104 deletions(-) diff --git a/doc/specs/stdlib_sorting.md b/doc/specs/stdlib_sorting.md index 59c484f1d..9b9af9d7e 100644 --- a/doc/specs/stdlib_sorting.md +++ b/doc/specs/stdlib_sorting.md @@ -25,9 +25,9 @@ module's `string_type` type. ## Overview of the module The module `stdlib_sorting` defines several public entities, one -default integer parameter, `int_size`, and three overloaded -subroutines: `ORD_SORT`, `SORT`, and `SORT_INDEX`. The -overloaded subroutines also each have seven specific names for +default integer parameter, `int_size`, and four overloaded +subroutines: `ORD_SORT`, `SORT`, `RADIX_SORT` and `SORT_INDEX`. The +overloaded subroutines also each have several specific names for versions corresponding to different types of array arguments. ### The `int_size` parameter @@ -47,9 +47,11 @@ data: * `SORT_INDEX` is based on `ORD_SORT`, but in addition to sorting the input array, it returns indices that map the original array to its sorted version. This enables related arrays to be re-ordered in the - same way; and + same way; * `SORT` is intended to sort simple arrays of intrinsic data - that are effectively unordered before the sort. + that are effectively unordered before the sort; +* `RADIX_SORT` is intended to sort fixed width intrinsic data + types (integers and reals). #### Licensing @@ -196,6 +198,18 @@ magnitude slower than `ORD_SORT`. Its memory requirements are also low, being of order O(Ln(N)), while the memory requirements of `ORD_SORT` and `SORT_INDEX` are of order O(N). +#### The `RADIX_SORT` subroutine + +`RADIX_SORT` is a implementation of LSD [radix sort](https://en.wikipedia.org/wiki/Radix_sort), +using `256` (one byte) as the radix. It only works for fixed width data, +thus integers and reals. `RADIX_SORT` is always of O(N) runtime performance +for any input data. For large and random data, it is about five (or more) +times faster than other sort subroutines. + +The `RADIX_SORT` needs a buffer that have same size of the input data. +Your can provide it using `work` arguement, if not the subroutine will +allocate the buffer and deallocate before return. + ### Specifications of the `stdlib_sorting` procedures #### `ord_sort` - sorts an input array @@ -318,6 +332,55 @@ element of `array` is a `NaN`. Sorting of `CHARACTER(*)` and {!example/sorting/example_sort.f90!} ``` +#### `radix_sort` - sorts an input array + +##### Status + +Experimental + +##### Description + +Returns an input array with the elements sorted in order of increasing, or +decreasing, value. + +##### Syntax + +`call [[stdlib_sorting(module):radix_sort(interface)]]( array[, work, reverse] )` + +##### Class + +Generic subroutine. + +##### Arguments + +`array` : shall be a rank one array of any of the types: +`integer(int8)`, `integer(int16)`, `integer(int32)`, `integer(int64)`, +`real(sp)`, `real(dp)`. It is an `intent(inout)` argument. On return its +input elements will be sorted in order of non-decreasing value. + +`work` (optional): shall be a rank one array of the same type as +array, and shall have at least `size(array)` elements. It is an +`intent(inout)` argument, and its contents on return are undefined. + +`reverse` (optional): shall be a scalar of type default logical. It +is an `intent(in)` argument. If present with a value of `.true.` then +`array` will be sorted in order of non-increasing values in unstable +order. Otherwise index will sort `array` in order of non-decreasing +values in unstable order. + +##### Notes + +`SORT` implements a LSD radix sort algorithm with a `256` radix. For any +input data it provides `O(N)` run time performance. If `array` is of +any type `REAL` the order of its elements on return undefined if any +element of `array` is a `NaN`. + +##### Example + +```fortran +{!example/sorting/example_radix_sort.f90!} +``` + #### `sort_index` - creates an array of sorting indices for an input array, while also sorting the array. ##### Status @@ -496,107 +559,126 @@ of size `16**3`, with characters drawn from the set "a"-"p": random order. These benchmarks have been performed on two different compilers, both -on a MacBook Pro, featuring a 2.3 GHz Quad-Core Intel Core i5, with 8 -GB 2133 MHz LPDDR3 memory. The first compiler was GNU Fortran -(GCC) 10.2.0, with the following results: +on WSL with Ubuntu-20.04, Intel(R) Core(TM) i7-10700 CPU @ 2.9GHz, with +32 GB DDR4 memory. The first compiler is GNU Fortran (GCC) 9.4.0, with +the following results. | Type | Elements | Array Name | Method | Time (s) | |-------------|----------|-----------------|-------------|-----------| -| Integer | 65536 | Blocks | Ord_Sort | 0.000579 | -| Integer | 65536 | Decreasing | Ord_Sort | 0.000212 | -| Integer | 65536 | Identical | Ord_Sort | 0.000165 | -| Integer | 65536 | Increasing | Ord_Sort | 0.000131 | -| Integer | 65536 | Random dense | Ord_Sort | 0.009991 | -| Integer | 65536 | Random order | Ord_Sort | 0.008574 | -| Integer | 65536 | Random sparse | Ord_Sort | 0.010504 | -| Integer | 65536 | Random 3 | Ord_Sort | 0.000532 | -| Integer | 65536 | Random 10 | Ord_Sort | 0.000315 | -| Character | 65536 | Char. Decrease | Ord_Sort | 0.001041 | -| Character | 65536 | Char. Increase | Ord_Sort | 0.000584 | -| Character | 65536 | Char. Random | Ord_Sort | 0.026273 | -| String_type | 4096 | String Decrease | Ord_Sort | 0.001202 | -| String_type | 4096 | String Increase | Ord_Sort | 0.000758 | -| String_type | 4096 | String Random | Ord_Sort | 0.018180 | -| Integer | 65536 | Blocks | Sort | 0.005073 | -| Integer | 65536 | Decreasing | Sort | 0.005830 | -| Integer | 65536 | Identical | Sort | 0.007372 | -| Integer | 65536 | Increasing | Sort | 0.002421 | -| Integer | 65536 | Random dense | Sort | 0.007006 | -| Integer | 65536 | Random order | Sort | 0.007211 | -| Integer | 65536 | Random sparse | Sort | 0.007109 | -| Integer | 65536 | Random 3 | Sort | 0.012232 | -| Integer | 65536 | Random 10 | Sort | 0.017345 | -| Character | 65536 | Char. Decrease | Sort | 0.031350 | -| Character | 65536 | Char. Increase | Sort | 0.011606 | -| Character | 65536 | Char. Random | Sort | 0.022440 | -| String_type | 4096 | String Decrease | Sort | 0.026539 | -| String_type | 4096 | String Increase | Sort | 0.009755 | -| String_type | 4096 | String Random | Sort | 0.016218 | -| Integer | 65536 | Blocks | Sort_Index | 0.000953 | -| Integer | 65536 | Decreasing | Sort_Index | 0.000418 | -| Integer | 65536 | Identical | Sort_Index | 0.000264 | -| Integer | 65536 | Increasing | Sort_Index | 0.000262 | -| Integer | 65536 | Random dense | Sort_Index | 0.009563 | -| Integer | 65536 | Random order | Sort_Index | 0.009592 | -| Integer | 65536 | Random sparse | Sort_Index | 0.009691 | -| Integer | 65536 | Random 3 | Sort_Index | 0.000781 | -| Integer | 65536 | Random 10 | Sort_Index | 0.000455 | -| Character | 65536 | Char. Decrease | Sort_Index | 0.001189 | -| Character | 65536 | Char. Increase | Sort_Index | 0.000752 | -| Character | 65536 | Char. Random | Sort_Index | 0.025767 | -| String_type | 4096 | String Decrease | Sort_Index | 0.001411 | -| String_type | 4096 | String Increase | Sort_Index | 0.000761 | -| String_type | 4096 | String Random | Sort_Index | 0.018202 | - -The second compiler was Intel(R) Fortran Intel(R) 64 Compiler Classic -for applications running on Intel(R) 64, Version 2021.2.0 Build -20210228_000000, with the following results: +| Integer | 65536 | Blocks | Ord_Sort | 0.001048 | +| Integer | 65536 | Decreasing | Ord_Sort | 0.000204 | +| Integer | 65536 | Identical | Ord_Sort | 0.000097 | +| Integer | 65536 | Increasing | Ord_Sort | 0.000096 | +| Integer | 65536 | Random dense | Ord_Sort | 0.006580 | +| Integer | 65536 | Random order | Ord_Sort | 0.006886 | +| Integer | 65536 | Random sparse | Ord_Sort | 0.006821 | +| Integer | 65536 | Random 3 | Ord_Sort | 0.000461 | +| Integer | 65536 | Random 10 | Ord_Sort | 0.000226 | +| Character | 65536 | Char. Decrease | Ord_Sort | 0.000824 | +| Character | 65536 | Char. Increase | Ord_Sort | 0.000370 | +| Character | 65536 | Char. Random | Ord_Sort | 0.016020 | +| String_type | 4096 | String Decrease | Ord_Sort | 0.000465 | +| String_type | 4096 | String Increase | Ord_Sort | 0.000169 | +| String_type | 4096 | String Random | Ord_Sort | 0.004194 | +| Integer | 65536 | Blocks | Radix_Sort | 0.001610 | +| Integer | 65536 | Decreasing | Radix_Sort | 0.001076 | +| Integer | 65536 | Identical | Radix_Sort | 0.001074 | +| Integer | 65536 | Increasing | Radix_Sort | 0.001060 | +| Integer | 65536 | Random dense | Radix_Sort | 0.001161 | +| Integer | 65536 | Random order | Radix_Sort | 0.001069 | +| Integer | 65536 | Random sparse | Radix_Sort | 0.001005 | +| Integer | 65536 | Random 3 | Radix_Sort | 0.001057 | +| Integer | 65536 | Random 10 | Radix_Sort | 0.001046 | +| Integer | 65536 | rand-real32 | Radix_Sort | 0.001429 | +| Integer | 65536 | Blocks | Sort | 0.004269 | +| Integer | 65536 | Decreasing | Sort | 0.005108 | +| Integer | 65536 | Identical | Sort | 0.006257 | +| Integer | 65536 | Increasing | Sort | 0.002093 | +| Integer | 65536 | Random dense | Sort | 0.006032 | +| Integer | 65536 | Random order | Sort | 0.006030 | +| Integer | 65536 | Random sparse | Sort | 0.006126 | +| Integer | 65536 | Random 3 | Sort | 0.007930 | +| Integer | 65536 | Random 10 | Sort | 0.014729 | +| Character | 65536 | Char. Decrease | Sort | 0.020623 | +| Character | 65536 | Char. Increase | Sort | 0.008028 | +| Character | 65536 | Char. Random | Sort | 0.014258 | +| String_type | 4096 | String Decrease | Sort | 0.005542 | +| String_type | 4096 | String Increase | Sort | 0.001987 | +| String_type | 4096 | String Random | Sort | 0.003267 | +| Integer | 65536 | Blocks | Sort_Index | 0.000686 | +| Integer | 65536 | Decreasing | Sort_Index | 0.000529 | +| Integer | 65536 | Identical | Sort_Index | 0.000218 | +| Integer | 65536 | Increasing | Sort_Index | 0.000214 | +| Integer | 65536 | Random dense | Sort_Index | 0.008044 | +| Integer | 65536 | Random order | Sort_Index | 0.008042 | +| Integer | 65536 | Random sparse | Sort_Index | 0.008148 | +| Integer | 65536 | Random 3 | Sort_Index | 0.000677 | +| Integer | 65536 | Random 10 | Sort_Index | 0.000387 | +| Character | 65536 | Char. Decrease | Sort_Index | 0.000932 | +| Character | 65536 | Char. Increase | Sort_Index | 0.000487 | +| Character | 65536 | Char. Random | Sort_Index | 0.017231 | +| String_type | 4096 | String Decrease | Sort_Index | 0.000489 | +| String_type | 4096 | String Increase | Sort_Index | 0.000183 | +| String_type | 4096 | String Random | Sort_Index | 0.004102 | + +The second compiler is Intel(R) Fortran Intel(R) 64 Compiler Classic +for applications running on Intel(R) 64, Version 2021.7.0 Build +20220726_000000, with the following results: | Type | Elements | Array Name | Method | Time (s) | |-------------|----------|-----------------|-------------|-----------| -| Integer | 65536 | Blocks | Ord_Sort | 0.000267 | -| Integer | 65536 | Decreasing | Ord_Sort | 0.000068 | -| Integer | 65536 | Identical | Ord_Sort | 0.000056 | -| Integer | 65536 | Increasing | Ord_Sort | 0.000056 | -| Integer | 65536 | Random dense | Ord_Sort | 0.004615 | -| Integer | 65536 | Random order | Ord_Sort | 0.006325 | -| Integer | 65536 | Random sparse | Ord_Sort | 0.004601 | -| Integer | 65536 | Random 3 | Ord_Sort | 0.000193 | -| Integer | 65536 | Random 10 | Ord_Sort | 0.000101 | -| Character | 65536 | Char. Decrease | Ord_Sort | 0.001009 | -| Character | 65536 | Char. Increase | Ord_Sort | 0.000529 | -| Character | 65536 | Char. Random | Ord_Sort | 0.024547 | -| String_type | 4096 | String Decrease | Ord_Sort | 0.003381 | -| String_type | 4096 | String Increase | Ord_Sort | 0.000133 | -| String_type | 4096 | String Random | Ord_Sort | 0.051985 | -| Integer | 65536 | Blocks | Sort | 0.001614 | -| Integer | 65536 | Decreasing | Sort | 0.001783 | -| Integer | 65536 | Identical | Sort | 0.002111 | -| Integer | 65536 | Increasing | Sort | 0.000674 | -| Integer | 65536 | Random dense | Sort | 0.003574 | -| Integer | 65536 | Random order | Sort | 0.003296 | -| Integer | 65536 | Random sparse | Sort | 0.003380 | -| Integer | 65536 | Random 3 | Sort | 0.003623 | -| Integer | 65536 | Random 10 | Sort | 0.006839 | -| Character | 65536 | Char. Decrease | Sort | 0.032564 | -| Character | 65536 | Char. Increase | Sort | 0.012346 | -| Character | 65536 | Char. Random | Sort | 0.022932 | -| String_type | 4096 | String Decrease | Sort | 0.082140 | -| String_type | 4096 | String Increase | Sort | 0.029591 | -| String_type | 4096 | String Random | Sort | 0.043078 | -| Integer | 65536 | Blocks | Sort_Index | 0.000848 | -| Integer | 65536 | Decreasing | Sort_Index | 0.000103 | -| Integer | 65536 | Identical | Sort_Index | 0.000102 | -| Integer | 65536 | Increasing | Sort_Index | 0.000066 | -| Integer | 65536 | Random dense | Sort_Index | 0.006434 | -| Integer | 65536 | Random order | Sort_Index | 0.005941 | -| Integer | 65536 | Random sparse | Sort_Index | 0.005957 | -| Integer | 65536 | Random 3 | Sort_Index | 0.000326 | -| Integer | 65536 | Random 10 | Sort_Index | 0.000175 | -| Character | 65536 | Char. Decrease | Sort_Index | 0.001082 | -| Character | 65536 | Char. Increase | Sort_Index | 0.000468 | -| Character | 65536 | Char. Random | Sort_Index | 0.023100 | -| String_type | 4096 | String Decrease | Sort_Index | 0.003292 | -| String_type | 4096 | String Increase | Sort_Index | 0.000122 | -| String_type | 4096 | String Random | Sort_Index | 0.049155 | - +| Integer | 65536 | Blocks | Ord_Sort | 0.000135 | +| Integer | 65536 | Decreasing | Ord_Sort | 0.000053 | +| Integer | 65536 | Identical | Ord_Sort | 0.000033 | +| Integer | 65536 | Increasing | Ord_Sort | 0.000034 | +| Integer | 65536 | Random dense | Ord_Sort | 0.003291 | +| Integer | 65536 | Random order | Ord_Sort | 0.003546 | +| Integer | 65536 | Random sparse | Ord_Sort | 0.003313 | +| Integer | 65536 | Random 3 | Ord_Sort | 0.000145 | +| Integer | 65536 | Random 10 | Ord_Sort | 0.000070 | +| Character | 65536 | Char. Decrease | Ord_Sort | 0.000696 | +| Character | 65536 | Char. Increase | Ord_Sort | 0.000338 | +| Character | 65536 | Char. Random | Ord_Sort | 0.015255 | +| String_type | 4096 | String Decrease | Ord_Sort | 0.001276 | +| String_type | 4096 | String Increase | Ord_Sort | 0.000153 | +| String_type | 4096 | String Random | Ord_Sort | 0.024705 | +| Integer | 65536 | Blocks | Radix_Sort | 0.001038 | +| Integer | 65536 | Decreasing | Radix_Sort | 0.000910 | +| Integer | 65536 | Identical | Radix_Sort | 0.000441 | +| Integer | 65536 | Increasing | Radix_Sort | 0.000803 | +| Integer | 65536 | Random dense | Radix_Sort | 0.000363 | +| Integer | 65536 | Random order | Radix_Sort | 0.000741 | +| Integer | 65536 | Random sparse | Radix_Sort | 0.000384 | +| Integer | 65536 | Random 3 | Radix_Sort | 0.000877 | +| Integer | 65536 | Random 10 | Radix_Sort | 0.000801 | +| Integer | 65536 | rand-real32 | Radix_Sort | 0.000604 | +| Integer | 65536 | Blocks | Sort | 0.001342 | +| Integer | 65536 | Decreasing | Sort | 0.001391 | +| Integer | 65536 | Identical | Sort | 0.001485 | +| Integer | 65536 | Increasing | Sort | 0.000447 | +| Integer | 65536 | Random dense | Sort | 0.002778 | +| Integer | 65536 | Random order | Sort | 0.002896 | +| Integer | 65536 | Random sparse | Sort | 0.003136 | +| Integer | 65536 | Random 3 | Sort | 0.002996 | +| Integer | 65536 | Random 10 | Sort | 0.005752 | +| Character | 65536 | Char. Decrease | Sort | 0.021973 | +| Character | 65536 | Char. Increase | Sort | 0.008391 | +| Character | 65536 | Char. Random | Sort | 0.015155 | +| String_type | 4096 | String Decrease | Sort | 0.034014 | +| String_type | 4096 | String Increase | Sort | 0.010464 | +| String_type | 4096 | String Random | Sort | 0.015748 | +| Integer | 65536 | Blocks | Sort_Index | 0.000381 | +| Integer | 65536 | Decreasing | Sort_Index | 0.000085 | +| Integer | 65536 | Identical | Sort_Index | 0.000046 | +| Integer | 65536 | Increasing | Sort_Index | 0.000046 | +| Integer | 65536 | Random dense | Sort_Index | 0.004020 | +| Integer | 65536 | Random order | Sort_Index | 0.004059 | +| Integer | 65536 | Random sparse | Sort_Index | 0.004073 | +| Integer | 65536 | Random 3 | Sort_Index | 0.000215 | +| Integer | 65536 | Random 10 | Sort_Index | 0.000101 | +| Character | 65536 | Char. Decrease | Sort_Index | 0.000680 | +| Character | 65536 | Char. Increase | Sort_Index | 0.000356 | +| Character | 65536 | Char. Random | Sort_Index | 0.016231 | +| String_type | 4096 | String Decrease | Sort_Index | 0.001219 | +| String_type | 4096 | String Increase | Sort_Index | 0.000125 | +| String_type | 4096 | String Random | Sort_Index | 0.018631 | diff --git a/src/stdlib_sorting.fypp b/src/stdlib_sorting.fypp index c0c034c51..7420816e2 100644 --- a/src/stdlib_sorting.fypp +++ b/src/stdlib_sorting.fypp @@ -430,6 +430,7 @@ module stdlib_sorting !! The generic subroutine interface implementing the LSD radix sort algorithm, !! see https://en.wikipedia.org/wiki/Radix_sort for more details. !! It is always O(N) in sorting random data, but need a O(N) buffer. +!! ([Specification](../page/specs/stdlib_sorting.html#radix_sort-sorts-an-input-array)) !! pure module subroutine int8_radix_sort(array, reverse) diff --git a/test/sorting/test_sorting.f90 b/test/sorting/test_sorting.f90 index 03b4ac5d1..c3de24c2f 100644 --- a/test/sorting/test_sorting.f90 +++ b/test/sorting/test_sorting.f90 @@ -38,7 +38,7 @@ module test_sorting string_rand(0:string_size-1) integer(int32) :: dummy(0:test_size-1) - real(sp) :: real_dummy(0:test_size-1) + real(sp) :: real_dummy(0:test_size-1) character(len=4) :: char_dummy(0:char_size-1) type(string_type) :: string_dummy(0:string_size-1) integer(int_size) :: index(0:max(test_size, char_size, string_size)-1) From 2ca2e4c6588905ffde554fe5f9abe65af1a30b1c Mon Sep 17 00:00:00 2001 From: 0382 <18322825326@163.com> Date: Sun, 21 May 2023 10:07:11 +0800 Subject: [PATCH 6/7] Apply suggestions from code review Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_sorting.md | 10 +++++----- example/sorting/example_radix_sort.f90 | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/specs/stdlib_sorting.md b/doc/specs/stdlib_sorting.md index 9b9af9d7e..47ea9c665 100644 --- a/doc/specs/stdlib_sorting.md +++ b/doc/specs/stdlib_sorting.md @@ -201,13 +201,13 @@ low, being of order O(Ln(N)), while the memory requirements of #### The `RADIX_SORT` subroutine `RADIX_SORT` is a implementation of LSD [radix sort](https://en.wikipedia.org/wiki/Radix_sort), -using `256` (one byte) as the radix. It only works for fixed width data, +using `256` as the radix. It only works for fixed width data, thus integers and reals. `RADIX_SORT` is always of O(N) runtime performance for any input data. For large and random data, it is about five (or more) times faster than other sort subroutines. The `RADIX_SORT` needs a buffer that have same size of the input data. -Your can provide it using `work` arguement, if not the subroutine will +Your can provide it using `work` argument, if not the subroutine will allocate the buffer and deallocate before return. ### Specifications of the `stdlib_sorting` procedures @@ -362,7 +362,7 @@ input elements will be sorted in order of non-decreasing value. array, and shall have at least `size(array)` elements. It is an `intent(inout)` argument, and its contents on return are undefined. -`reverse` (optional): shall be a scalar of type default logical. It +`reverse` (optional): shall be a scalar of type default `logical`. It is an `intent(in)` argument. If present with a value of `.true.` then `array` will be sorted in order of non-increasing values in unstable order. Otherwise index will sort `array` in order of non-decreasing @@ -370,9 +370,9 @@ values in unstable order. ##### Notes -`SORT` implements a LSD radix sort algorithm with a `256` radix. For any +`radix_sort` implements a LSD radix sort algorithm with a `256` radix. For any input data it provides `O(N)` run time performance. If `array` is of -any type `REAL` the order of its elements on return undefined if any +any type `real` the order of its elements on return undefined if any element of `array` is a `NaN`. ##### Example diff --git a/example/sorting/example_radix_sort.f90 b/example/sorting/example_radix_sort.f90 index 21b355c8f..10b8dcd13 100644 --- a/example/sorting/example_radix_sort.f90 +++ b/example/sorting/example_radix_sort.f90 @@ -20,7 +20,7 @@ program example_radix_sort arrf64 = [1.0_dp/x, 0.0_dp, 0.0_dp/x, -1.0_dp/x, -0.0_dp, 1.0_dp, -1.0_dp, 3.45_dp, -3.14_dp, 3.44_dp] call radix_sort(arrf64) print *, arrf64 - ! Expect output: + ! Expected output: ! nan, -inf, -3.14, -1.0, -0.0, 0.0, 1.0, 3.44, 3.45, inf ! Note: the position of nan is undefined end program example_radix_sort From de88ad21b515bc3df81277bfc37b6990fe4c5281 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Thu, 15 Jun 2023 14:02:21 -0400 Subject: [PATCH 7/7] Update doc/specs/stdlib_sorting.md Co-authored-by: 0382 <18322825326@163.com> --- doc/specs/stdlib_sorting.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_sorting.md b/doc/specs/stdlib_sorting.md index 47ea9c665..f37cefaec 100644 --- a/doc/specs/stdlib_sorting.md +++ b/doc/specs/stdlib_sorting.md @@ -200,7 +200,7 @@ low, being of order O(Ln(N)), while the memory requirements of #### The `RADIX_SORT` subroutine -`RADIX_SORT` is a implementation of LSD [radix sort](https://en.wikipedia.org/wiki/Radix_sort), +`RADIX_SORT` is a implementation of LSD [radix sort](https://www.growingwiththeweb.com/sorting/radix-sort-lsd/), using `256` as the radix. It only works for fixed width data, thus integers and reals. `RADIX_SORT` is always of O(N) runtime performance for any input data. For large and random data, it is about five (or more)