diff --git a/doc/specs/stdlib_sorting.md b/doc/specs/stdlib_sorting.md
index 3a44d84f8..dbf395a10 100644
--- a/doc/specs/stdlib_sorting.md
+++ b/doc/specs/stdlib_sorting.md
@@ -46,6 +46,9 @@ data:
* `ORD_SORT` is intended to sort simple arrays of intrinsic data
that have significant sections that were partially ordered before
the sort;
+* `SORT_ADJOINT` is based on `ORD_SORT`, but in addition to sorting the
+ input array, it re-orders a second array of the same size
+ according to the same permutations;
* `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
@@ -60,10 +63,10 @@ data:
The Fortran Standard Library is distributed under the MIT
License. However components of the library may be based on code with
additional licensing restrictions. In particular `ORD_SORT`,
-`SORT_INDEX`, and `SORT` are translations of codes with their
+`SORT_ADJOINT`, `SORT_INDEX`, and `SORT` are translations of codes with their
own distribution restrictions.
-The `ORD_SORT` and `SORT_INDEX` subroutines are essentially
+The `ORD_SORT`, `SORT_ADJOINT` and `SORT_INDEX` subroutines are essentially
translations to Fortran 2008 of the `"Rust" sort` of the Rust Language
distributed as part of
[`slice.rs`](https://github.com/rust-lang/rust/blob/90eb44a5897c39e3dff9c7e48e3973671dcd9496/src/liballoc/slice.rs).
@@ -140,6 +143,23 @@ argument or allocated internally on the stack.
Arrays can be also sorted in a decreasing order by providing the argument `reverse
= .true.`.
+#### The `SORT_ADJOINT` subroutine
+
+The `SORT` and `ORD_SORT` subroutines can sort rank 1 isolated
+arrays of intrinsic types, but do nothing for the coordinated sorting
+of related data, e.g., a related rank 1 array. Therefore the module
+provides a subroutine, `SORT_ADJOINT`, that re-order such a rank 1 array
+according to the same permutations as for the input array based on the `ORD_SORT` algorithm,
+in addition to sorting the input array.
+
+The logic of `SORT_ADJOINT` parallels that of `ORD_SORT`, with
+additional housekeeping to keep the associated array consistent with
+the sorted positions of the input array. Because of this additional
+housekeeping it has slower runtime performance than `ORD_SORT`.
+`SORT_ADJOINT` requires the use of two "scratch" arrays, that may be
+provided as optional `work` and `iwork` arguments or allocated
+internally on the stack.
+
#### The `SORT_INDEX` subroutine
The `SORT` and `ORD_SORT` subroutines can sort rank 1 isolated
@@ -198,7 +218,7 @@ factor of six. Still, even when it shows enhanced performance, its
performance on partially sorted data is typically an order of
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).
+`ORD_SORT`, `SORT_ADJOINT` and `SORT_INDEX` are of order O(N).
#### The `RADIX_SORT` subroutine
@@ -385,6 +405,85 @@ element of `array` is a `NaN`.
{!example/sorting/example_radix_sort.f90!}
```
+#### `sort_adjoint` - sorts an associated array
+according to the same permutations as for the input array.
+
+##### Status
+
+Experimental
+
+##### Description
+
+Returns the input `array` sorted in the direction requested while
+retaining order stability, and an associated array whose elements are
+sorted according to the same permutations as for the input `array`.
+
+##### Syntax
+
+`call ` [[stdlib_sorting(module):sort_adjoint(interface)]] `( array, adjoint_array[, work, iwork, 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)`, `real(qp)`, `character(*)`, `type(string_type)`,
+`type(bitset_64)`, or `type(bitset_large)`.
+It is an `intent(inout)` argument. On input it
+will be an array whose sorting indices are to be determined. On return
+it will be the sorted array.
+
+`adjoint_array`: shall be a rank one `integer` or `real` array of
+the size of `array`. It is an `intent(inout)` argument. On return it
+shall have values that are the indices needed to sort the original
+array in the desired direction.
+
+`work` (optional): shall be a rank one array of any of the same type as
+`array`, and shall have at least `size(array)/2` elements. It is an
+`intent(out)` argument. It is intended to be used as "scratch"
+memory for internal record keeping. If associated with an array in
+static storage, its use can significantly reduce the stack memory
+requirements for the code. Its contents on return are undefined.
+
+`iwork` (optional): shall be a rank one integer array of the same kind
+of the array `adjoint_array`, and shall have at least `size(array)/2` elements. It
+is an `intent(out)` argument. It is intended to be used as "scratch"
+memory for internal record keeping. If associated with an array in
+static storage, its use can significantly reduce the stack memory
+requirements for the code. 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 stable
+order. Otherwise `array` will be sorted in order of non-decreasing
+values in stable order.
+
+##### Notes
+
+`SORT_ADJOINT` implements the hybrid sorting algorithm of `ORD_SORT`,
+keeping the values of `adjoint_array` consistent with the elements of `array`
+as it is sorted. As a `merge sort` based algorithm, it is a stable
+sorting comparison algorithm. The optional `work` and `iwork` arrays
+replace "scratch" memory that would otherwise be allocated on the
+stack. If `array` is of any kind of `REAL` the order of the elements in
+`adjoint_array` and `array` on return are undefined if any element of `array`
+is a `NaN`. Sorting of `CHARACTER(*)` and `STRING_TYPE` arrays are
+based on the operator `>`, and not on the function `LGT`.
+
+It should be emphasized that the order of `array` will typically be
+different on return
+
+##### Examples
+
+Sorting a rank one array with `sort_adjoint`:
+
+```Fortran
+{!example/sorting/example_sort_adjoint.f90!}
+```
+
#### `sort_index` - creates an array of sorting indices for an input array, while also sorting the array.
##### Status
diff --git a/example/sorting/CMakeLists.txt b/example/sorting/CMakeLists.txt
index 4628ce20c..6d64ea2f1 100644
--- a/example/sorting/CMakeLists.txt
+++ b/example/sorting/CMakeLists.txt
@@ -1,5 +1,6 @@
ADD_EXAMPLE(ord_sort)
ADD_EXAMPLE(sort)
+ADD_EXAMPLE(sort_adjoint)
ADD_EXAMPLE(sort_index)
ADD_EXAMPLE(radix_sort)
ADD_EXAMPLE(sort_bitset)
diff --git a/example/sorting/example_sort_adjoint.f90 b/example/sorting/example_sort_adjoint.f90
new file mode 100644
index 000000000..d311fea4e
--- /dev/null
+++ b/example/sorting/example_sort_adjoint.f90
@@ -0,0 +1,15 @@
+program example_sort_adjoint
+ use stdlib_sorting, only: sort_adjoint
+ implicit none
+ integer, allocatable :: array(:)
+ real, allocatable :: adjoint(:)
+
+ array = [5, 4, 3, 1, 10, 4, 9]
+ allocate(adjoint, source=real(array))
+
+ call sort_adjoint(array, adjoint)
+
+ print *, array !print [1, 3, 4, 4, 5, 9, 10]
+ print *, adjoint !print [1., 3., 4., 4., 5., 9., 10.]
+
+end program example_sort_adjoint
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 579b70c72..5b92b7eff 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -38,7 +38,7 @@ set(fppFiles
stdlib_sorting.fypp
stdlib_sorting_ord_sort.fypp
stdlib_sorting_sort.fypp
- stdlib_sorting_sort_index.fypp
+ stdlib_sorting_sort_adjoint.fypp
stdlib_specialfunctions_gamma.fypp
stdlib_stats.fypp
stdlib_stats_corr.fypp
diff --git a/src/stdlib_sorting.fypp b/src/stdlib_sorting.fypp
index e0bb93827..c675e5f3f 100644
--- a/src/stdlib_sorting.fypp
+++ b/src/stdlib_sorting.fypp
@@ -13,6 +13,7 @@
#! This approach allows us to have the same code for all input types.
#:set IRSCB_TYPES_ALT_NAME = INT_TYPES_ALT_NAME + REAL_TYPES_ALT_NAME + STRING_TYPES_ALT_NAME + CHAR_TYPES_ALT_NAME &
& + BITSET_TYPES_ALT_NAME
+#:set IR_INDEX_TYPES_ALT_NAME = INT_TYPES_ALT_NAME + REAL_TYPES_ALT_NAME
!! Licensing:
!!
@@ -292,6 +293,76 @@ module stdlib_sorting
!! ! Sort the random data
!! call radix_sort( array )
!! ...
+!!```
+
+ public sort_adjoint
+!! Version: experimental
+!!
+!! The generic subroutine implementing the `SORT_ADJ` algorithm to
+!! return an adjoint array whose elements are sorted in the same order
+!! as the input array in the
+!! desired direction. It is primarily intended to be used to sort a
+!! rank 1 `integer` or `real` array based on the values of a component of the array.
+!! Its use has the syntax:
+!!
+!! call sort_adjoint( array, adjoint_array[, work, iwork, 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)`,
+!! `real(real128)`, `character(*)`, `type(string_type)`,
+!! `type(bitset_64)`, `type(bitset_large)`. If both the
+!! type of `array` is real and at least one of the elements is a `NaN`,
+!! then the ordering of the `array` and `adjoint_array` results is undefined.
+!! Otherwise it is defined to be as specified by reverse.
+!!
+!! * adjoint_array: a rank 1 `integer` or `real` array. It is an `intent(inout)`
+!! argument. Its size shall be the
+!! same as `array`. On return, its elements are sorted in the same order
+!! as the input `array` in the direction specified by `reverse`.
+!!
+!! * work (optional): shall be a rank 1 array of the same type as
+!! `array`, and shall have at least `size(array)/2` elements. It is an
+!! `intent(out)` argument to be used as "scratch" memory
+!! for internal record keeping. If associated with an array in static
+!! storage, its use can significantly reduce the stack memory requirements
+!! for the code. Its value on return is undefined.
+!!
+!! * iwork (optional): shall be a rank 1 integer array of the same type as `adjoint_array`,
+!! and shall have at least `size(array)/2` elements. It is an
+!! `intent(out)` argument to be used as "scratch" memory
+!! for internal record keeping. If associated with an array in static
+!! storage, its use can significantly reduce the stack memory requirements
+!! for the code. Its value on return is 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 stable
+!! order. Otherwise `array` will be sorted in order of non-decreasing
+!! values in stable order.
+!!
+!!#### Examples
+!!
+!! Sorting a related rank one array:
+!!
+!!```Fortran
+!!program example_sort_adjoint
+!! use stdlib_sorting, only: sort_adjoint
+!! implicit none
+!! integer, allocatable :: array(:)
+!! real, allocatable :: adj(:)
+!!
+!! array = [5, 4, 3, 1, 10, 4, 9]
+!! allocate(adj, source=real(array))
+!!
+!! call sort_adjoint(array, adj)
+!!
+!! print *, array !print [1, 3, 4, 4, 5, 9, 10]
+!! print *, adj !print [1., 3., 4., 4., 5., 9., 10.]
+!!
+!!end program example_sort_adjoint
!!```
public sort_index
@@ -505,6 +576,43 @@ module stdlib_sorting
end interface sort
+ interface sort_adjoint
+!! Version: experimental
+!!
+!! The generic subroutine interface implementing the `SORT_ADJ` algorithm,
+!! based on the `"Rust" sort` algorithm found in `slice.rs`
+!! https://github.com/rust-lang/rust/blob/90eb44a5897c39e3dff9c7e48e3973671dcd9496/src/liballoc/slice.rs#L2159
+!! but modified to return an array of indices that would provide a stable
+!! sort of the rank one `ARRAY` input.
+!! ([Specification](../page/specs/stdlib_sorting.html#sort_adjoint-creates-an-array-of-sorting-indices-for-an-input-array-while-also-sorting-the-array))
+!!
+!! The indices by default correspond to a
+!! non-decreasing sort, but if the optional argument `REVERSE` is present
+!! with a value of `.TRUE.` the indices correspond to a non-increasing sort.
+
+#:for ti, tii, namei in IR_INDEX_TYPES_ALT_NAME
+ #:for t1, t2, name1 in IRSCB_TYPES_ALT_NAME
+ module subroutine ${name1}$_${namei}$_sort_adjoint( array, adjoint_array, work, iwork, &
+ reverse )
+!! Version: experimental
+!!
+!! `${name1}$_${namei}$_sort_adjoint( array, adjoint_array[, work, iwork, reverse] )` sorts
+!! an input `ARRAY` of type `${t1}$`
+!! using a hybrid sort based on the `"Rust" sort` algorithm found in `slice.rs`
+!! and returns the sorted `ARRAY` and an array `INDEX` of indices in the
+!! order that would sort the input `ARRAY` in the desired direction.
+ ${t1}$, intent(inout) :: array(0:)
+ ${ti}$, intent(inout) :: adjoint_array(0:)
+ ${t2}$, intent(out), optional :: work(0:)
+ ${ti}$, intent(out), optional :: iwork(0:)
+ logical, intent(in), optional :: reverse
+ end subroutine ${name1}$_${namei}$_sort_adjoint
+
+ #:endfor
+#:endfor
+
+ end interface sort_adjoint
+
interface sort_index
!! Version: experimental
!!
@@ -521,7 +629,24 @@ module stdlib_sorting
#:for ki, ti, namei in INT_INDEX_TYPES_ALT_NAME
#:for t1, t2, name1 in IRSCB_TYPES_ALT_NAME
- module subroutine ${name1}$_sort_index_${namei}$( array, index, work, iwork, &
+!> Version: experimental
+!>
+!> `${name1}$_sort_index_${namei}$( array, index[, work, iwork, reverse] )` sorts
+!> an input `ARRAY` of type `${t1}$`
+!> using a hybrid sort based on the `"Rust" sort` algorithm found in `slice.rs`
+!> and returns the sorted `ARRAY` and an array `INDEX` of indices in the
+!> order that would sort the input `ARRAY` in the desired direction.
+ module procedure ${name1}$_sort_index_${namei}$
+ #:endfor
+#:endfor
+
+ end interface sort_index
+
+contains
+
+#:for ki, ti, namei in INT_INDEX_TYPES_ALT_NAME
+ #:for t1, t2, name1 in IRSCB_TYPES_ALT_NAME
+ subroutine ${name1}$_sort_index_${namei}$( array, index, work, iwork, &
reverse )
!! Version: experimental
!!
@@ -535,12 +660,29 @@ module stdlib_sorting
${t2}$, intent(out), optional :: work(0:)
${ti}$, intent(out), optional :: iwork(0:)
logical, intent(in), optional :: reverse
+
+ integer(int_index) :: array_size, i
+
+ array_size = size(array, kind=int_index)
+
+ if ( array_size > huge(index)) then
+ error stop "Too many entries for the kind of index."
+ end if
+
+ if ( array_size > size(index, kind=int_index) ) then
+ error stop "Too many entries for the size of index."
+ end if
+
+ do i = 0, array_size-1
+ index(i) = int(i+1, kind=${ki}$)
+ end do
+
+ call sort_adjoint(array, index, work, iwork, reverse)
+
end subroutine ${name1}$_sort_index_${namei}$
#:endfor
#:endfor
- end interface sort_index
-
end module stdlib_sorting
diff --git a/src/stdlib_sorting_ord_sort.fypp b/src/stdlib_sorting_ord_sort.fypp
index b96ea295a..c77e1c797 100644
--- a/src/stdlib_sorting_ord_sort.fypp
+++ b/src/stdlib_sorting_ord_sort.fypp
@@ -117,6 +117,8 @@ contains
integer :: stat
array_size = size( array, kind=int_index )
+
+! If necessary allocate buffers to serve as scratch memory.
if ( present(work) ) then
if ( size(work, kind=int_index) < array_size/2 ) then
error stop "${name1}$_${sname}$_ord_sort: work array is too small."
diff --git a/src/stdlib_sorting_sort_index.fypp b/src/stdlib_sorting_sort_adjoint.fypp
similarity index 84%
rename from src/stdlib_sorting_sort_index.fypp
rename to src/stdlib_sorting_sort_adjoint.fypp
index cc2afe9cf..a8b99b034 100644
--- a/src/stdlib_sorting_sort_index.fypp
+++ b/src/stdlib_sorting_sort_adjoint.fypp
@@ -5,13 +5,12 @@
#:set CHAR_TYPES_ALT_NAME = list(zip(["character(len=*)"], ["character(len=:)"], ["character(len=len(array))"], ["char"]))
#:set BITSET_TYPES_ALT_NAME = list(zip(BITSET_TYPES, BITSET_TYPES, BITSET_TYPES, BITSET_KINDS))
-#:set INT_INDEX_TYPES_ALT_NAME = list(zip(["int_index", "int_index_low"], ["integer(int_index)", "integer(int_index_low)"], ["default", "low"]))
-
#! For better code reuse in fypp, make lists that contain the input types,
#! with each having output types and a separate name prefix for subroutines
#! This approach allows us to have the same code for all input types.
#:set IRSCB_TYPES_ALT_NAME = INT_TYPES_ALT_NAME + REAL_TYPES_ALT_NAME + STRING_TYPES_ALT_NAME + CHAR_TYPES_ALT_NAME &
& + BITSET_TYPES_ALT_NAME
+#:set IR_INDEX_TYPES_ALT_NAME = INT_TYPES_ALT_NAME + REAL_TYPES_ALT_NAME
!! Licensing:
!!
@@ -43,7 +42,7 @@
!! TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
!! SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!!
-!! The generic subroutine, `SORT_INDEX`, is substantially a translation to
+!! The generic subroutine, `SORT_ADJ`, is substantially a translation to
!! Fortran 2008 of the `"Rust" sort` sorting routines in
!! [`slice.rs`](https://github.com/rust-lang/rust/blob/90eb44a5897c39e3dff9c7e48e3973671dcd9496/src/liballoc/slice.rs)
!! The `rust sort` implementation is distributed with the header:
@@ -62,16 +61,16 @@
!! of modified versions of the code in the Fortran Standard Library under
!! the MIT license.
-submodule(stdlib_sorting) stdlib_sorting_sort_index
+submodule(stdlib_sorting) stdlib_sorting_sort_adjoint
implicit none
contains
-#:for ki, ti, namei in INT_INDEX_TYPES_ALT_NAME
+#:for ki, ti, tii, namei in IR_INDEX_TYPES_ALT_NAME
#:for t1, t2, t3, name1 in IRSCB_TYPES_ALT_NAME
- module subroutine ${name1}$_sort_index_${namei}$( array, index, work, iwork, reverse )
+ module subroutine ${name1}$_${namei}$_sort_adjoint( array, adjoint_array, work, iwork, reverse )
! A modification of `${name1}$_ord_sort` to return an array of indices that
! would perform a stable sort of the `ARRAY` as input, and also sort `ARRAY`
! as desired. The indices by default
@@ -96,33 +95,21 @@ contains
! estimation of the optimal `run size` as suggested in Tim Peters'
! original `listsort.txt`, and the optional `work` and `iwork` arrays to be
! used as scratch memory.
-
${t1}$, intent(inout) :: array(0:)
- ${ti}$, intent(out) :: index(0:)
+ ${ti}$, intent(inout) :: adjoint_array(0:)
${t3}$, intent(out), optional :: work(0:)
${ti}$, intent(out), optional :: iwork(0:)
logical, intent(in), optional :: reverse
${t2}$, allocatable :: buf(:)
${ti}$, allocatable :: ibuf(:)
- integer(int_index) :: array_size, i, stat
+ integer(int_index) :: array_size
+ integer(int_index) :: stat
array_size = size(array, kind=int_index)
- if ( array_size > huge(index)) then
- error stop "Too many entries for the kind of index."
- end if
-
- if ( array_size > size(index, kind=int_index) ) then
- error stop "Too many entries for the size of index."
- end if
-
- do i = 0, array_size-1
- index(i) = int(i+1, kind=${ki}$)
- end do
-
if ( optval(reverse, .false.) ) then
- call reverse_segment( array, index )
+ call reverse_segment( array, adjoint_array )
end if
! If necessary allocate buffers to serve as scratch memory.
@@ -134,13 +121,14 @@ contains
if ( size(iwork, kind=int_index) < array_size/2 ) then
error stop "iwork array is too small."
endif
- call merge_sort( array, index, work, iwork )
+ call merge_sort( array, adjoint_array, work, iwork )
else
allocate( ibuf(0:array_size/2-1), stat=stat )
- if ( stat /= 0 ) error stop "Allocation of index buffer failed."
- call merge_sort( array, index, work, ibuf )
+ if ( stat /= 0 ) error stop "Allocation of adjoint_array buffer failed."
+ call merge_sort( array, adjoint_array, work, ibuf )
end if
else
+! Allocate a buffer to use as scratch memory.
#:if t1[0:4] == "char"
allocate( ${t3}$ :: buf(0:array_size/2-1), &
stat=stat )
@@ -152,16 +140,16 @@ contains
if ( size(iwork, kind=int_index) < array_size/2 ) then
error stop "iwork array is too small."
endif
- call merge_sort( array, index, buf, iwork )
+ call merge_sort( array, adjoint_array, buf, iwork )
else
allocate( ibuf(0:array_size/2-1), stat=stat )
- if ( stat /= 0 ) error stop "Allocation of index buffer failed."
- call merge_sort( array, index, buf, ibuf )
+ if ( stat /= 0 ) error stop "Allocation of adjoint_array buffer failed."
+ call merge_sort( array, adjoint_array, buf, ibuf )
end if
end if
if ( optval(reverse, .false.) ) then
- call reverse_segment( array, index )
+ call reverse_segment( array, adjoint_array )
end if
contains
@@ -187,28 +175,28 @@ contains
end function calc_min_run
- pure subroutine insertion_sort( array, index )
+ pure subroutine insertion_sort( array, adjoint_array )
! Sorts `ARRAY` using an insertion sort, while maintaining consistency in
! location of the indices in `INDEX` to the elements of `ARRAY`.
${t1}$, intent(inout) :: array(0:)
- ${ti}$, intent(inout) :: index(0:)
+ ${ti}$, intent(inout) :: adjoint_array(0:)
integer(int_index) :: i, j
- ${ti}$ :: key_index
+ ${ti}$ :: key_adjoint_array
${t3}$ :: key
do j=1, size(array, kind=int_index)-1
key = array(j)
- key_index = index(j)
+ key_adjoint_array = adjoint_array(j)
i = j - 1
do while( i >= 0 )
if ( array(i) <= key ) exit
array(i+1) = array(i)
- index(i+1) = index(i)
+ adjoint_array(i+1) = adjoint_array(i)
i = i - 1
end do
array(i+1) = key
- index(i+1) = key_index
+ adjoint_array(i+1) = key_adjoint_array
end do
end subroutine insertion_sort
@@ -265,36 +253,36 @@ contains
end function collapse
- pure subroutine insert_head( array, index )
+ pure subroutine insert_head( array, adjoint_array )
! Inserts `array(0)` into the pre-sorted sequence `array(1:)` so that the
! whole `array(0:)` becomes sorted, copying the first element into
! a temporary variable, iterating until the right place for it is found.
! copying every traversed element into the slot preceding it, and finally,
! copying data from the temporary variable into the resulting hole.
-! Consistency of the indices in `index` with the elements of `array`
+! Consistency of the indices in `adjoint_array` with the elements of `array`
! are maintained.
${t1}$, intent(inout) :: array(0:)
- ${ti}$, intent(inout) :: index(0:)
+ ${ti}$, intent(inout) :: adjoint_array(0:)
${t3}$ :: tmp
integer(int_index) :: i
- ${ti}$ :: tmp_index
+ ${ti}$ :: tmp_adjoint_array
tmp = array(0)
- tmp_index = index(0)
+ tmp_adjoint_array = adjoint_array(0)
find_hole: do i=1, size(array, kind=int_index)-1
if ( array(i) >= tmp ) exit find_hole
array(i-1) = array(i)
- index(i-1) = index(i)
+ adjoint_array(i-1) = adjoint_array(i)
end do find_hole
array(i-1) = tmp
- index(i-1) = tmp_index
+ adjoint_array(i-1) = tmp_adjoint_array
end subroutine insert_head
- subroutine merge_sort( array, index, buf, ibuf )
+ subroutine merge_sort( array, adjoint_array, buf, ibuf )
! The Rust merge sort borrows some (but not all) of the ideas from TimSort,
! which is described in detail at
! (http://svn.python.org/projects/python/trunk/Objects/listsort.txt).
@@ -312,11 +300,11 @@ contains
! runs(i - 1)%len + runs(i)%len`
!
! The invariants ensure that the total running time is `O(n log n)`
-! worst-case. Consistency of the indices in `index` with the elements of
+! worst-case. Consistency of the indices in `adjoint_array` with the elements of
! `array` are maintained.
${t1}$, intent(inout) :: array(0:)
- ${ti}$, intent(inout) :: index(0:)
+ ${ti}$, intent(inout) :: adjoint_array(0:)
${t3}$, intent(inout) :: buf(0:)
${ti}$, intent(inout) :: ibuf(0:)
@@ -331,7 +319,7 @@ contains
min_run = calc_min_run( array_size )
if ( array_size <= min_run ) then
- if ( array_size >= 2 ) call insertion_sort( array, index )
+ if ( array_size >= 2 ) call insertion_sort( array, adjoint_array )
return
end if
@@ -354,7 +342,7 @@ contains
start = start - 1
end do Descending
call reverse_segment( array(start:finish), &
- index(start:finish) )
+ adjoint_array(start:finish) )
else
Ascending: do while( start > 0 )
if ( array(start) < array(start-1) ) exit Ascending
@@ -367,7 +355,7 @@ contains
Insert: do while ( start > 0 )
if ( finish - start >= min_run - 1 ) exit Insert
start = start - 1
- call insert_head( array(start:finish), index(start:finish) )
+ call insert_head( array(start:finish), adjoint_array(start:finish) )
end do Insert
if ( start == 0 .and. finish == array_size - 1 ) return
@@ -386,7 +374,7 @@ contains
call merge( array( left % base: &
right % base + right % len - 1 ), &
left % len, buf, &
- index( left % base: &
+ adjoint_array( left % base: &
right % base + right % len - 1 ), ibuf )
runs(r) = run_type( base = left % base, &
@@ -402,7 +390,7 @@ contains
end subroutine merge_sort
- pure subroutine merge( array, mid, buf, index, ibuf )
+ pure subroutine merge( array, mid, buf, adjoint_array, ibuf )
! Merges the two non-decreasing runs `ARRAY(0:MID-1)` and `ARRAY(MID:)`
! using `BUF` as temporary storage, and stores the merged runs into
! `ARRAY(0:)`. `MID` must be > 0, and < `SIZE(ARRAY)-1`. Buffer `BUF`
@@ -410,7 +398,7 @@ contains
${t1}$, intent(inout) :: array(0:)
integer(int_index), intent(in) :: mid
${t3}$, intent(inout) :: buf(0:)
- ${ti}$, intent(inout) :: index(0:)
+ ${ti}$, intent(inout) :: adjoint_array(0:)
${ti}$, intent(inout) :: ibuf(0:)
integer(int_index) :: array_len, i, j, k
@@ -424,44 +412,44 @@ contains
if ( mid <= array_len - mid ) then ! The left run is shorter.
buf(0:mid-1) = array(0:mid-1)
- ibuf(0:mid-1) = index(0:mid-1)
+ ibuf(0:mid-1) = adjoint_array(0:mid-1)
i = 0
j = mid
merge_lower: do k = 0, array_len-1
if ( buf(i) <= array(j) ) then
array(k) = buf(i)
- index(k) = ibuf(i)
+ adjoint_array(k) = ibuf(i)
i = i + 1
if ( i >= mid ) exit merge_lower
else
array(k) = array(j)
- index(k) = index(j)
+ adjoint_array(k) = adjoint_array(j)
j = j + 1
if ( j >= array_len ) then
array(k+1:) = buf(i:mid-1)
- index(k+1:) = ibuf(i:mid-1)
+ adjoint_array(k+1:) = ibuf(i:mid-1)
exit merge_lower
end if
end if
end do merge_lower
else ! The right run is shorter
buf(0:array_len-mid-1) = array(mid:array_len-1)
- ibuf(0:array_len-mid-1) = index(mid:array_len-1)
+ ibuf(0:array_len-mid-1) = adjoint_array(mid:array_len-1)
i = mid - 1
j = array_len - mid -1
merge_upper: do k = array_len-1, 0, -1
if ( buf(j) >= array(i) ) then
array(k) = buf(j)
- index(k) = ibuf(j)
+ adjoint_array(k) = ibuf(j)
j = j - 1
if ( j < 0 ) exit merge_upper
else
array(k) = array(i)
- index(k) = index(i)
+ adjoint_array(k) = adjoint_array(i)
i = i - 1
if ( i < 0 ) then
array(0:k-1) = buf(0:j)
- index(0:k-1) = ibuf(0:j)
+ adjoint_array(0:k-1) = ibuf(0:j)
exit merge_upper
end if
end if
@@ -470,10 +458,10 @@ contains
end subroutine merge
- pure subroutine reverse_segment( array, index )
+ pure subroutine reverse_segment( array, adjoint_array )
! Reverse a segment of an array in place
${t1}$, intent(inout) :: array(0:)
- ${ti}$, intent(inout) :: index(0:)
+ ${ti}$, intent(inout) :: adjoint_array(0:)
${ti}$ :: itemp
integer(int_index) :: lo, hi
@@ -485,18 +473,19 @@ contains
temp = array(lo)
array(lo) = array(hi)
array(hi) = temp
- itemp = index(lo)
- index(lo) = index(hi)
- index(hi) = itemp
+ itemp = adjoint_array(lo)
+ adjoint_array(lo) = adjoint_array(hi)
+ adjoint_array(hi) = itemp
lo = lo + 1
hi = hi - 1
end do
end subroutine reverse_segment
- end subroutine ${name1}$_sort_index_${namei}$
+ end subroutine ${name1}$_${namei}$_sort_adjoint
#:endfor
#:endfor
-end submodule stdlib_sorting_sort_index
+end submodule stdlib_sorting_sort_adjoint
+
diff --git a/test/sorting/test_sorting.fypp b/test/sorting/test_sorting.fypp
index 1418a032f..7e74927ea 100644
--- a/test/sorting/test_sorting.fypp
+++ b/test/sorting/test_sorting.fypp
@@ -1,11 +1,14 @@
#:include "common.fypp"
+#:set INT_TYPES_ALT_NAME = list(zip(INT_TYPES, INT_TYPES, INT_KINDS))
+#:set REAL_TYPES_ALT_NAME = list(zip(REAL_TYPES, REAL_TYPES, REAL_KINDS))
+
#:set INT_INDEX_TYPES_ALT_NAME = list(zip(["int_index", "int_index_low"], ["integer(int_index)", "integer(int_index_low)"], ["default", "low"]))
module test_sorting
use, intrinsic :: iso_fortran_env, only: compiler_version, error_unit
- use stdlib_kinds, only: int32, int64, dp, sp
- use stdlib_sorting
+ use stdlib_kinds, only: int8, int16, int32, int64, dp, sp
+ use stdlib_sorting, only: sort, sort_index, sort_adjoint, ord_sort, radix_sort, int_index, int_index_low
use stdlib_string_type, only: string_type, assignment(=), operator(>), &
operator(<), write(formatted)
use stdlib_bitsets, only: bitset_64, bitset_large, &
@@ -105,6 +108,13 @@ contains
new_unittest('string_sort_indexes_${namei}$', test_string_sort_indexes_${namei}$), &
new_unittest('bitset_large_sort_indexes_${namei}$', test_bitsetl_sort_indexes_${namei}$), &
new_unittest('bitset_64_sort_indexes_${namei}$', test_bitset64_sort_indexes_${namei}$), &
+#:endfor
+#:for ki, ti, namei in INT_TYPES_ALT_NAME
+ new_unittest('int_sort_adjointes_${namei}$', test_int_sort_adjointes_${namei}$), &
+ new_unittest('char_sort_adjointes_${namei}$', test_char_sort_adjointes_${namei}$), &
+ new_unittest('string_sort_adjointes_${namei}$', test_string_sort_adjointes_${namei}$), &
+ new_unittest('bitset_large_sort_adjointes_${namei}$', test_bitsetl_sort_adjointes_${namei}$), &
+ new_unittest('bitset_64_sort_adjointes_${namei}$', test_bitset64_sort_adjointes_${namei}$), &
#:endfor
new_unittest('int_ord_sorts', test_int_ord_sorts) &
]
@@ -1542,6 +1552,350 @@ contains
end subroutine test_bitset64_sort_index_${namei}$
#:endfor
+#:for ki, ti, namei in INT_TYPES_ALT_NAME
+ subroutine test_int_sort_adjointes_${namei}$(error)
+ !> Error handling
+ type(error_type), allocatable, intent(out) :: error
+ integer(int64) :: i
+ integer(int32), allocatable :: d1(:)
+ ${ti}$, allocatable :: adjoint(:)
+ logical :: ltest
+
+ call test_int_sort_adjoint_${namei}$( blocks, "Blocks", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_int_sort_adjoint_${namei}$( decrease, "Decreasing", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_int_sort_adjoint_${namei}$( identical, "Identical", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_int_sort_adjoint_${namei}$( increase, "Increasing", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_int_sort_adjoint_${namei}$( rand1, "Random dense", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_int_sort_adjoint_${namei}$( rand2, "Random order", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_int_sort_adjoint_${namei}$( rand0, "Random sparse", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_int_sort_adjoint_${namei}$( rand3, "Random 3", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_int_sort_adjoint_${namei}$( rand10, "Random 10", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ d1 = [10, 2, -3, -4, 6, -6, 7, -8, 9, 0, 1, 20]
+ allocate( adjoint(size(d1)))
+ adjoint = int([(i, i=1, size(d1))], kind=${namei}$)
+ call sort_adjoint( d1, adjoint )
+ call verify_sort( d1, ltest, i )
+ call check(error, ltest)
+
+ end subroutine test_int_sort_adjointes_${namei}$
+
+ subroutine test_int_sort_adjoint_${namei}$( a, a_name, ltest )
+ integer(int32), intent(inout) :: a(:)
+ character(*), intent(in) :: a_name
+ logical, intent(out) :: ltest
+
+ integer(int64) :: t0, t1, tdiff
+ real(dp) :: rate
+ ${ti}$ :: adjoint(size(a))
+ ${ti}$ :: iwork(size(a))
+ integer(int64) :: i, j
+ logical :: valid
+
+ ltest = .true.
+
+ tdiff = 0
+ do i = 1, repeat
+ dummy = a
+ adjoint = int([(j, j=1_int64, size(a, kind=1_int64))], kind=${namei}$)
+ call system_clock( t0, rate )
+ call sort_adjoint( dummy, adjoint, work, iwork )
+ 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_ADJOINT did not sort " // a_name // "."
+ write(*,*) 'i = ', i
+ write(*,'(a18, 2i7)') 'a(i-1:i) = ', a(i-1:i)
+ end if
+ write( lun, '("| Integer |", 1x, i7, 2x, "|", 1x, a15, " |", ' // &
+ 'a12, " |", F10.6, " |" )' ) &
+ test_size, a_name, "Sort_adjoint", tdiff/rate
+
+ !reverse
+ dummy = a
+ adjoint = int([(j, j=1_int64, size(a, kind=1_int64))], kind=${namei}$)
+ call sort_adjoint( dummy, adjoint, work, iwork, reverse=.true. )
+
+ call verify_reverse_sort( dummy, valid, i )
+ ltest = (ltest .and. valid)
+ if ( .not. valid ) then
+ write( *, * ) "SORT_ADJOINT did not reverse sort " // &
+ a_name // "."
+ write(*,*) 'i = ', i
+ write(*,'(a18, 2i7)') 'a(i-1:i) = ', a(i-1:i)
+ end if
+
+ end subroutine test_int_sort_adjoint_${namei}$
+
+ subroutine test_char_sort_adjointes_${namei}$(error)
+ !> Error handling
+ type(error_type), allocatable, intent(out) :: error
+ logical :: ltest
+
+ call test_char_sort_adjoint_${namei}$( char_decrease, "Char. Decrease", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_char_sort_adjoint_${namei}$( char_increase, "Char. Increase", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_char_sort_adjoint_${namei}$( char_rand, "Char. Random", ltest )
+ call check(error, ltest)
+
+ end subroutine test_char_sort_adjointes_${namei}$
+
+ subroutine test_char_sort_adjoint_${namei}$( a, a_name, ltest )
+ character(len=4), intent(in) :: a(0:)
+ character(*), intent(in) :: a_name
+ logical, intent(out) :: ltest
+
+ ${ti}$ :: adjoint(size(a))
+ ${ti}$ :: iwork(size(a))
+ integer(int64) :: t0, t1, tdiff
+ real(dp) :: rate
+ integer(int64) :: i, j
+ logical :: valid
+
+ ltest = .true.
+
+ tdiff = 0
+ do i = 1, repeat
+ char_dummy = a
+ adjoint = int([(j, j=1_int64, size(a, kind=1_int64))], kind=${namei}$)
+ call system_clock( t0, rate )
+
+ call sort_adjoint( char_dummy, adjoint, char_work, iwork )
+
+ call system_clock( t1, rate )
+
+ tdiff = tdiff + t1 - t0
+ end do
+ tdiff = tdiff/repeat
+
+ call verify_char_sort( char_dummy, valid, i )
+
+ ltest = (ltest .and. valid)
+ if ( .not. valid ) then
+ write( *, * ) "SORT_ADJ did not sort " // a_name // "."
+ write(*,*) 'i = ', i
+ write(*,'(a17, 2(1x,a4))') 'char_dummy(i-1:i) = ', char_dummy(i-1:i)
+ end if
+ write( lun, '("| Character |", 1x, i7, 2x, "|", 1x, a15, " |", ' // &
+ 'a12, " |", F10.6, " |" )' ) &
+ char_size, a_name, "Sort_adjoint", tdiff/rate
+
+ end subroutine test_char_sort_adjoint_${namei}$
+
+ subroutine test_string_sort_adjointes_${namei}$(error)
+ !> Error handling
+ type(error_type), allocatable, intent(out) :: error
+ logical :: ltest
+
+ call test_string_sort_adjoint_${namei}$( string_decrease, "String Decrease", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_string_sort_adjoint_${namei}$( string_increase, "String Increase", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_string_sort_adjoint_${namei}$( string_rand, "String Random", ltest )
+ call check(error, ltest)
+
+ end subroutine test_string_sort_adjointes_${namei}$
+
+ subroutine test_string_sort_adjoint_${namei}$( a, a_name, ltest )
+ type(string_type), intent(in) :: a(0:)
+ character(*), intent(in) :: a_name
+ logical, intent(out) :: ltest
+
+ ${ti}$ :: adjoint(size(a))
+ ${ti}$ :: iwork(size(a))
+ integer(int64) :: t0, t1, tdiff
+ real(dp) :: rate
+ integer(int64) :: i, j
+ logical :: valid
+
+ ltest = .true.
+
+ tdiff = 0
+ do i = 1, repeat
+ string_dummy = a
+ adjoint = int([(j, j=1_int64, size(a, kind=1_int64))], kind=${namei}$)
+ call system_clock( t0, rate )
+ call sort_adjoint( string_dummy, adjoint, string_work, iwork )
+ call system_clock( t1, rate )
+ tdiff = tdiff + t1 - t0
+ end do
+ tdiff = tdiff/repeat
+
+ call verify_string_sort( string_dummy, valid, i )
+ ltest = (ltest .and. valid)
+ if ( .not. valid ) then
+ write( *, * ) "SORT_ADJOINT did not sort " // a_name // "."
+ write(*,*) 'i = ', i
+ write(*,'(a17, 2(1x,a4))') 'string_dummy(i-1:i) = ', &
+ string_dummy(i-1:i)
+ end if
+ write( lun, '("| String_type |", 1x, i7, 2x, "|", 1x, a15, " |", ' // &
+ 'a12, " |", F10.6, " |" )' ) &
+ string_size, a_name, "Sort_adjoint", tdiff/rate
+
+ end subroutine test_string_sort_adjoint_${namei}$
+
+ subroutine test_bitsetl_sort_adjointes_${namei}$(error)
+ !> Error handling
+ type(error_type), allocatable, intent(out) :: error
+ logical :: ltest
+
+ call test_bitsetl_sort_adjoint_${namei}$( bitsetl_decrease, "Bitset Decrease", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_bitsetl_sort_adjoint_${namei}$( bitsetl_increase, "Bitset Increase", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_bitsetl_sort_adjoint_${namei}$( bitsetl_rand, "Bitset Random", ltest )
+ call check(error, ltest)
+
+ end subroutine test_bitsetl_sort_adjointes_${namei}$
+
+ subroutine test_bitsetl_sort_adjoint_${namei}$( a, a_name, ltest )
+ type(bitset_large), intent(in) :: a(0:)
+ character(*), intent(in) :: a_name
+ logical, intent(out) :: ltest
+
+ ${ti}$ :: adjoint(size(a))
+ ${ti}$ :: iwork(size(a))
+ integer(int64) :: t0, t1, tdiff
+ real(dp) :: rate
+ integer(int64) :: i, j
+ logical :: valid
+ character(:), allocatable :: bin_im1, bin_i
+
+ ltest = .true.
+
+ tdiff = 0
+ do i = 1, repeat
+ bitsetl_dummy = a
+ adjoint = int([(j, j=1_int64, size(a, kind=1_int64))], kind=${namei}$)
+ call system_clock( t0, rate )
+ call sort_adjoint( bitsetl_dummy, adjoint, bitsetl_work, iwork )
+ call system_clock( t1, rate )
+ tdiff = tdiff + t1 - t0
+ end do
+ tdiff = tdiff/repeat
+
+ call verify_bitsetl_sort( bitsetl_dummy, valid, i )
+ ltest = (ltest .and. valid)
+ if ( .not. valid ) then
+ write( *, * ) "SORT_ADJOINT did not sort " // a_name // "."
+ write(*,*) 'i = ', i
+ call bitsetl_dummy(i-1)%to_string(bin_im1)
+ call bitsetl_dummy(i)%to_string(bin_i)
+ write(*,'(a, 2(a:,1x))') 'bitsetl_dummy(i-1:i) = ', &
+ bin_im1, bin_i
+ end if
+ write( lun, '("| Bitset_large |", 1x, i7, 2x, "|", 1x, a15, " |", ' // &
+ 'a12, " |", F10.6, " |" )' ) &
+ bitset_size, a_name, "Sort_adjoint", tdiff/rate
+
+ end subroutine test_bitsetl_sort_adjoint_${namei}$
+
+ subroutine test_bitset64_sort_adjointes_${namei}$(error)
+ !> Error handling
+ type(error_type), allocatable, intent(out) :: error
+ logical :: ltest
+
+ call test_bitset64_sort_adjoint_${namei}$( bitset64_decrease, "Bitset Decrease", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_bitset64_sort_adjoint_${namei}$( bitset64_increase, "Bitset Increase", ltest )
+ call check(error, ltest)
+ if (allocated(error)) return
+
+ call test_bitset64_sort_adjoint_${namei}$( bitset64_rand, "Bitset Random", ltest )
+ call check(error, ltest)
+
+ end subroutine test_bitset64_sort_adjointes_${namei}$
+
+ subroutine test_bitset64_sort_adjoint_${namei}$( a, a_name, ltest )
+ type(bitset_64), intent(in) :: a(0:)
+ character(*), intent(in) :: a_name
+ logical, intent(out) :: ltest
+
+ ${ti}$ :: adjoint(size(a))
+ ${ti}$ :: iwork(size(a))
+ integer(int64) :: t0, t1, tdiff
+ real(dp) :: rate
+ integer(int64) :: i, j
+ logical :: valid
+ character(:), allocatable :: bin_im1, bin_i
+
+ ltest = .true.
+
+ tdiff = 0
+ do i = 1, repeat
+ bitset64_dummy = a
+ adjoint = int([(j, j=1_int64, size(a, kind=1_int64))], kind=${namei}$)
+ call system_clock( t0, rate )
+ call sort_adjoint( bitset64_dummy, adjoint, bitset64_work, iwork )
+ call system_clock( t1, rate )
+ tdiff = tdiff + t1 - t0
+ end do
+ tdiff = tdiff/repeat
+
+ call verify_bitset64_sort( bitset64_dummy, valid, i )
+ ltest = (ltest .and. valid)
+ if ( .not. valid ) then
+ write( *, * ) "SORT_ADJOINT did not sort " // a_name // "."
+ write(*,*) 'i = ', i
+ call bitset64_dummy(i-1)%to_string(bin_im1)
+ call bitset64_dummy(i)%to_string(bin_i)
+ write(*,'(a, 2(a:,1x))') 'bitset64_dummy(i-1:i) = ', &
+ bin_im1, bin_i
+ end if
+ write( lun, '("| Bitset_64 |", 1x, i7, 2x, "|", 1x, a15, " |", ' // &
+ 'a12, " |", F10.6, " |" )' ) &
+ bitset_size, a_name, "Sort_adjoint", tdiff/rate
+
+ end subroutine test_bitset64_sort_adjoint_${namei}$
+#:endfor
+
subroutine verify_sort( a, valid, i )
integer(int32), intent(in) :: a(0:)
logical, intent(out) :: valid