diff --git a/doc/specs/stdlib_sorting.md b/doc/specs/stdlib_sorting.md index 59c484f1d..f37cefaec 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://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) +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` argument, 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 + +`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 +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/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..10b8dcd13 --- /dev/null +++ b/example/sorting/example_radix_sort.f90 @@ -0,0 +1,26 @@ +program example_radix_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(dp) :: x + real(dp), allocatable :: arrf64(:) + + 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_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 + ! 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 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..7420816e2 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,6 +424,50 @@ module stdlib_sorting #:endfor end interface ord_sort + 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. +!! ([Specification](../page/specs/stdlib_sorting.html#radix_sort-sorts-an-input-array)) +!! + + 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..9b23562f7 --- /dev/null +++ b/src/stdlib_sorting_radix_sort.f90 @@ -0,0 +1,466 @@ +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 +! For int8, radix sort becomes counting sort, so buffer is not needed + pure subroutine radix_sort_u8_helper(N, arr) + integer(kind=int_size), intent(in) :: N + integer(kind=int8), dimension(N), intent(inout) :: arr + integer(kind=int_size) :: i + integer :: bin_idx + integer(kind=int_size), 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=int_size), intent(in) :: N + integer(kind=int16), dimension(N), intent(inout) :: arr + integer(kind=int16), dimension(N), intent(inout) :: buf + integer(kind=int_size) :: i + integer :: b, b0, b1 + integer(kind=int_size), 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=int_size), intent(in) :: N + integer(kind=int32), dimension(N), intent(inout) :: arr + integer(kind=int32), dimension(N), intent(inout) :: buf + integer(kind=int_size) :: i + integer :: b, b0, b1, b2, b3 + integer(kind=int_size), 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=int_size), intent(in) :: N + integer(kind=int64), dimension(N), intent(inout) :: arr + integer(kind=int64), dimension(N), intent(inout) :: buffer + integer(kind=int_size) :: i + integer(kind=int64) :: b, b0, b1, b2, b3, b4, b5, b6, b7 + integer(kind=int_size), 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) :: 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 + item = array(i) + array(i) = array(N - i + 1) + array(N - i + 1) = 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=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=int_size) + if (present(work)) 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. + buffer => work + else + use_internal_buffer = .true. + 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 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 + 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=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=int_size) + if (present(work)) 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. + 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=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=int_size) + if (present(work)) 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. + 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) +! 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 + 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=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=int_size) + if (present(work)) 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. + 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=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=int_size) + if (present(work)) 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. + 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..c3de24c2f 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) @@ -64,6 +66,8 @@ 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('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), & @@ -116,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 @@ -446,6 +453,158 @@ 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( *, * ) "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 + + ! 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 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 @@ -903,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:) @@ -969,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