Fortran中的排序功能

Question

是否有一个Fortran库，它具有一个排序函数的实现，即函数ordering(list)（如Mathematica中的Ordering[]），它给出list中的位置，在该位置出现排序列表的每个连续元素？

我可以实现它，但我不想重新发明轮子（我的轮子可能远非完美......）。由于它是如此基本，我正在搜索包含此类列表操作但未能找到的列表。

你有什么建议吗？

Answer 1

由于我已经很久以前已经实现了这一点（它依赖于Bill Press等人的数字食谱书，并且大量借用），这里是Fortran中一个独立的实现：

module index_mod
    use, intrinsic :: iso_fortran_env, only: IK=>int32, RK=>real64
    implicit none
contains
    subroutine indexArrayReal(n,Array,Index)
        implicit none
        integer(IK), intent(in)  :: n
        real(RK)   , intent(in)  :: Array(n)
        integer(IK), intent(out) :: Index(n)
        integer(IK), parameter   :: nn=15, nstack=50
        integer(IK)              :: k,i,j,indext,jstack,l,r
        integer(IK)              :: istack(nstack)
        real(RK)                 :: a
        do j = 1,n
            Index(j) = j
        end do
        jstack=0
        l=1
        r=n
        do
            if (r-l < nn) then
                do j=l+1,r
                    indext=Index(j)
                    a=Array(indext)
                    do i=j-1,l,-1
                        if (Array(Index(i)) <= a) exit
                        Index(i+1)=Index(i)
                    end do
                    Index(i+1)=indext
                end do
                if (jstack == 0) return
                r=istack(jstack)
                l=istack(jstack-1)
                jstack=jstack-2
            else
                k=(l+r)/2
                call swap(Index(k),Index(l+1))
                call exchangeIndex(Index(l),Index(r))
                call exchangeIndex(Index(l+1),Index(r))
                call exchangeIndex(Index(l),Index(l+1))
                i=l+1
                j=r
                indext=Index(l+1)
                a=Array(indext)
                do
                    do
                        i=i+1
                        if (Array(Index(i)) >= a) exit
                    end do
                    do
                        j=j-1
                        if (Array(Index(j)) <= a) exit
                    end do
                    if (j < i) exit
                    call swap(Index(i),Index(j))
                end do
                Index(l+1)=Index(j)
                Index(j)=indext
                jstack=jstack+2
                if (jstack > nstack) then
                    write(*,*) 'NSTACK too small in indexArrayReal()'   ! xxx
                    error stop
                end if
                if (r-i+1 >= j-l) then
                    istack(jstack)=r
                    istack(jstack-1)=i
                    r=j-1
                else
                    istack(jstack)=j-1
                    istack(jstack-1)=l
                    l=i
                end if
            end if
        end do
    contains
        subroutine exchangeIndex(i,j)
            integer(IK), intent(inout) :: i,j
            integer(IK)                :: swp
            if (Array(j) < Array(i)) then
                swp=i
                i=j
                j=swp
            end if
        end subroutine exchangeIndex
        pure elemental subroutine swap(a,b)
            implicit none
            integer(IK), intent(inout) :: a,b
            integer(IK) :: dum
            dum=a
            a=b
            b=dum
        end subroutine swap
    end subroutine indexArrayReal
end module Index_mod

program Index_prog
    use Index_mod, only: IK, RK, indexArrayReal
    implicit none
    integer(IK), parameter  :: n = 5
    integer(IK)             :: Index(n)
    real(RK)                :: Array(n) = [ 1.,3.,4.,2.,-1. ]
    call indexArrayReal(n,Array,Index)
    write(*,*) "Index: ", Index
    write(*,*) "Array(Index): ", Array(Index)
end program Index_prog

使用GFortran 2008编译，这是输出：

$gfortran -std=f2008 *.f95 -o main
$main
 Index:            5           1           4           2           3
 Array(Index):   -1.0000000000000000        1.0000000000000000            2.0000000000000000        3.0000000000000000        4.0000000000000000     

上面的例程用于排序实值数组。要对整数数组进行排序，只需将real(RK) :: Array(n)接口中的subroutine indexArrayReal()更改为integer(IK) :: Array(n)即可。