[Fortran] 纯文本查看 复制代码
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=[3,2,5,1,4]
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