add test_real_radix_sorts

Author: 0382

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

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
 !!

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<width>_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<width>_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.