Add radix_sort

example/sorting/CMakeLists.txt

@@ -1,2 +1,3 @@
 ADD_EXAMPLE(ord_sort)
 ADD_EXAMPLE(sort)
+ADD_EXAMPLE(radix_sort)

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

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

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