rargsort.f90 Source File

module mod_rargsort

    implicit none
    private

    public :: rargsort

contains

! https://github.com/certik/fortran-utils/blob/master/src/sorting.f90

    subroutine rargsort(a,b,n_a)
        ! Returns the indices that would sort an array.
        !
        ! Arguments
        ! ---------
        !
        integer n_a
        real, intent(in):: a(n_a)   ! array of numbers
        integer, intent(out) :: b(n_a)         ! indices into the array 'a' that sort it
        !integer :: rargsort(size(a))
        !
        ! Example
        ! -------
        !
        ! rargsort([4.1_dp, 2.1_dp, 2.05_dp, -1.5_dp, 4.2_dp]) ! Returns [4, 3, 2, 1, 5]

        integer :: N                           ! number of numbers/vectors
        integer :: i,imin                      ! indices: i, i of smallest
        integer :: temp1                       ! temporary
        real :: temp2
        real :: a2(n_a)

        !write(*,*) a
        a2 = a
        N=size(a)
        do i = 1, N
            b(i) = i
        end do
        do i = 1, N-1
            ! find ith smallest in 'a'
            imin = minloc(a2(i:),1) + i - 1
            ! swap to position i in 'a' and 'b', if not already there
            if (imin /= i) then
                temp2 = a2(i); a2(i) = a2(imin); a2(imin) = temp2
                temp1 = b(i); b(i) = b(imin); b(imin) = temp1
            end if
        end do
        !rargsort=b
    end subroutine

    function rargsort_original(a) result(b)
        ! Returns the indices that would sort an array.
        !
        ! Arguments
        ! ---------
        !
        real, intent(in):: a(:)   ! array of numbers
        integer :: b(size(a))         ! indices into the array 'a' that sort it
        !
        ! Example
        ! -------
        !
        ! rargsort([4.1_dp, 2.1_dp, 2.05_dp, -1.5_dp, 4.2_dp]) ! Returns [4, 3, 2, 1, 5]

        integer :: N                           ! number of numbers/vectors
        integer :: i,imin                      ! indices: i, i of smallest
        integer :: temp1                       ! temporary
        real :: temp2
        real :: a2(size(a))
        a2 = a
        N=size(a)
        do i = 1, N
            b(i) = i
        end do
        do i = 1, N-1
            ! find ith smallest in 'a'
            imin = minloc(a2(i:),1) + i - 1
            ! swap to position i in 'a' and 'b', if not already there
            if (imin /= i) then
                temp2 = a2(i); a2(i) = a2(imin); a2(imin) = temp2
                temp1 = b(i); b(i) = b(imin); b(imin) = temp1
            end if
        end do
    end function

end module mod_rargsort