论坛 › Fortran算法研究 › Fortran调用C语言排序函数qsort

weixing1531 发表于 2024-5-26 00:44:36

Fortran调用C语言排序函数qsort

本帖最后由 weixing1531 于 2024-5-26 01:03 编辑

众所周知， Fortran标准目前没有排序子程序，但可通过4种途径进行排序：
1.祖传代码，如NR中sort.f90；
2.编译器扩展功能；
3.第三方库，如Dislin提供的sortr1子程序；
4.混合编程。

C语言标准有排序函数qsort，头文件为stdlib.h（绝对不公平）。
Fortran通过iso_c_binding模块调用C语言qsort函数 （傍大款）。

module sort
    implicit none

    interface swap !数据交换子程序通用接口 公有
      procedure :: swap_int !交换整数
      procedure :: swap_real !交换实数
    end interface
contains
!NR快速排序法 默认升序
subroutine sort_quick_real( arr, id ) !模块方法
    real, dimension(:), intent(inout) :: arr
    integer, dimension(:), intent(out), allocatable, optional :: id !排序前对应的索引值
    !最小值对应排序前下标为id(1) 最大值对应排序前下标为id(size(id))
    integer, parameter :: NN=15, NSTACK=50 !第3版修改nn=7, nstack=64 第2版为nn=15, nstack=50
    real :: a !临时变量
    integer :: i,j,k,jstack,l,r
    integer, dimension(NSTACK) :: istack
    integer, pointer :: b => null() !临时变量

    if ( present(id) ) then
      allocate(b,source=0) !分配内存
      id = [(i, i=1, size(arr))] !初值即排序前下标
    end if

    jstack=0
    l=1
    r=size(arr)

    do
      if (r-l < NN) then
      do j=l+1,r
          a=arr(j)
          if ( present(id) ) b=id(j) !新增
          do i=j-1,l,-1
            if (arr(i) < a) exit
            arr(i+1)=arr(i)
            if ( present(id) ) id(i+1)=id(i) !新增
          end do
          arr(i+1)=a
          if ( present(id) ) id(i+1)=b !新增
      end do
      if (jstack == 0) return
      r=istack(jstack)
      l=istack(jstack-1)
      jstack=jstack-2
      else
      k=(l+r)/2
      call swap(arr(k),arr(l+1))
      if ( present(id) ) call swap(id(k),id(l+1)) !交换元素索引值

      if(arr(r)<arr(l))then !小于号已重载
          call swap(arr(l),arr(r))
          if ( present(id) ) call swap(id(l),id(r)) !交换元素索引值
      end if

      if(arr(r)<arr(l+1))then !小于号已重载
          call swap(arr(l+1),arr(r))
          if ( present(id) ) call swap(id(l+1),id(r)) !交换元素索引值
      end if

      if(arr(l+1)<arr(l))then !小于号已重载
          call swap(arr(l),arr(l+1))
          if ( present(id) ) call swap(id(l),id(l+1)) !交换元素索引值
      end if
      i=l+1
      j=r
      a=arr(l+1)
      if ( present(id) ) b=id(l+1) !新增
      do
          do
            i=i+1
            if (a < arr(i)) exit !小于号已重载
          end do
          do
            j=j-1
            if (arr(j) < a) exit !小于号已重载
          end do
          if (j < i) exit
          call swap(arr(i),arr(j))
          if ( present(id) ) call swap(id(i),id(j)) !交换元素索引值
      end do
      arr(l+1)=arr(j)
      if ( present(id) ) id(l+1)=id(j) !新增
      arr(j)=a
      if ( present(id) ) id(j)=b !新增
      jstack=jstack+2
      if (jstack > NSTACK) error stop 'sort: NSTACK too small'

      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

    if ( present(id) ) deallocate(b) !释放内存
end subroutine sort_quick_real

subroutine swap_int(a,b) !交换整数
    integer,intent(inout)::a,b
    integer::temp

    temp=b
    b=a
    a=temp
end subroutine swap_int

subroutine swap_real(a,b) !交换浮点数
    real,intent(inout)::a,b
    real::temp

    temp=b
    b=a
    a=temp
end subroutine swap_real
end module sort

program test_sort
    use sort
    use iso_c_binding
    implicit none

    interface
      !void* qsort(void* arr,c_size_t n,c_size_t size,int (*comp)(const void* ,const void*));
      subroutine qsort(arr,n,size,comp) bind(c) !C语言qsort函数
            use iso_c_binding
            type(c_ptr),value::arr !void* arr
            integer(c_size_t),value::n,size !元素个数 元素所占大小
            type(c_funptr),intent(in),value::comp !C语言函数指针 用于比较
      end subroutine

      function comp_real(a,b) bind(c) !浮点数比较
            use iso_c_binding
            real(c_float),intent(in)::a,b
            integer(c_int)::comp_real
      end function

      function comp_int(a,b) bind(c) !整数比较
            use iso_c_binding
            integer(c_int),intent(in)::a,b
            integer(c_int)::comp_int
      end function
    end interface

    integer,parameter::N=5
    real::d(N) !原始数组
    real,target::s(N) !排序后数组 c_loc()必需
    integer,target::iarry(N) !排序后数组 c_loc()必需
    real::t1,t2
    type(c_funptr)::f_ptr !c语言函数指针

    call init_seed(999) !固定种子 每次计算时随机数组相同
    call RANDOM_NUMBER(d) !生成随机数组
    !以下实数类型
    s=d
    call cpu_time(t1)
    call sort_quick_real(s) !NR快速排序法
    call cpu_time(t2)
    write(*,*)"NR quick_sort:",t2-t1
    write(*,*)s

    s=d
    f_ptr=c_funloc(comp_real) !获取Fortran函数的C地址
    call cpu_time(t1)
    call qsort(c_loc(s),int(N,kind=c_size_t),c_sizeof(s(1)),f_ptr) !C语言标准库qsort 混合编程
    call cpu_time(t2)
    write(*,*)"qsort_float:",t2-t1
    write(*,*)s

    iarry=
    f_ptr=c_funloc(comp_int) !获取Fortran函数的C地址
    call cpu_time(t1)
    call qsort(c_loc(iarry),int(N,kind=c_size_t),c_sizeof(iarry(1)),f_ptr) !C语言标准库qsort 混合编程
    call cpu_time(t2)
    write(*,*)"qsort_int:",t2-t1
    write(*,*)iarry

    read(*,*)
contains
subroutine init_seed(i) !固定种子
    integer,intent(in)::i
    integer,allocatable::seed(:)
    integer::nn

    call random_seed(size=nn)
    allocate(seed(nn))
    seed = i
    call random_seed(put=seed)
    deallocate(seed)
end subroutine
end program test_sort

function comp_real(a,b) bind(c) !浮点数比较
    use iso_c_binding
    implicit none
    real(c_float),intent(in)::a,b
    integer(c_int)::comp_real

    comp_real=0
    if(a>b)then
      comp_real=1 !升序
    else if(a<b)then
      comp_real=-1
    end if
end function

function comp_int(a,b) bind(c) !整数比较
    use iso_c_binding
    implicit none
    integer(c_int),intent(in)::a,b
    integer(c_int)::comp_int

    comp_int=0
    if(a>b)then
      comp_int=1 !升序
    else if(a<b)then
      comp_int=-1
    end if
end function

li913 发表于 2024-5-27 09:13:15

intel fortran编译器提供 qsort

weixing1531 发表于 2024-5-27 13:27:22

li913 发表于 2024-5-27 09:13
intel fortran编译器提供 qsort

编译器扩展功能
GFortran没有

查看完整版本: Fortran调用C语言排序函数qsort