Addition of sort_adj to sort an rank 1 array in the same order as an input array (#849)

Author: perazz

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

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)

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

src/CMakeLists.txt

@@ -47,7 +47,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_sparse_constants.fypp
     stdlib_sparse_conversion.fypp
     stdlib_sparse_kinds.fypp

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

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."