From 2c3a25e20a6e924de0ac63ba6a353c1db01455b9 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 May 2021 10:45:37 +0200 Subject: [PATCH] Cleaned sorting --- src/utils/intel.f90 | 125 ++++++++++++++++++++++++++--- src/utils/sort.irp.f | 187 ++++++++++++++----------------------------- 2 files changed, 173 insertions(+), 139 deletions(-) diff --git a/src/utils/intel.f90 b/src/utils/intel.f90 index 681adde9..4c18af8d 100644 --- a/src/utils/intel.f90 +++ b/src/utils/intel.f90 @@ -1,13 +1,5 @@ module intel use, intrinsic :: iso_c_binding - interface - subroutine ippsSortRadixIndexGetBufferSize(len, dataType, pBufSize) bind(C, name='ippsSortRadixIndexGetBufferSize') - use iso_c_binding - integer, intent(in), value :: len - integer, intent(in), value :: dataType - integer, intent(out) :: pBufSize - end - end interface interface subroutine ippsSortAscend_32s_I(pSrc, len) bind(C, name='ippsSortAscend_32s_I') use iso_c_binding @@ -15,12 +7,86 @@ module intel integer, intent(inout) :: pSrc(len) end end interface + interface + subroutine ippsSortAscend_32f_I(pSrc, len) bind(C, name='ippsSortAscend_32f_I') + use iso_c_binding + integer, intent(in), value :: len + real, intent(inout) :: pSrc(len) + end + end interface + interface + subroutine ippsSortAscend_64s_I(pSrc, len) bind(C, name='ippsSortAscend_64s_I') + use iso_c_binding + integer*8, intent(in), value :: len + integer, intent(inout) :: pSrc(len) + end + end interface + interface + subroutine ippsSortAscend_64f_I(pSrc, len) bind(C, name='ippsSortAscend_64f_I') + use iso_c_binding + double precision, intent(in), value :: len + real, intent(inout) :: pSrc(len) + end + end interface + + interface + subroutine ippsSortRadixIndexGetBufferSize(len, dataType, pBufSize) bind(C, name='ippsSortRadixIndexGetBufferSize') + use iso_c_binding + integer, intent(in), value :: len + integer, intent(in), value :: dataType + integer, intent(out) :: pBufSize + end + end interface + + interface + subroutine ippsSortRadixAscend_16s_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_16s_I') + use iso_c_binding + integer, intent(in), value :: len + integer*2, intent(inout) :: pSrc(len) + character, intent(inout) :: pTmp(len) + end + end interface interface subroutine ippsSortRadixAscend_32s_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_32s_I') use iso_c_binding integer, intent(in), value :: len integer, intent(inout) :: pSrc(len) - integer, intent(inout) :: pTmp(len) + character, intent(inout) :: pTmp(len) + end + end interface + interface + subroutine ippsSortRadixAscend_32f_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_32f_I') + use iso_c_binding + integer, intent(in), value :: len + real, intent(inout) :: pSrc(len) + character, intent(inout) :: pTmp(len) + end + end interface + interface + subroutine ippsSortRadixAscend_64s_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_64s_I') + use iso_c_binding + integer, intent(in), value :: len + integer*8, intent(inout) :: pSrc(len) + character, intent(inout) :: pTmp(len) + end + end interface + interface + subroutine ippsSortRadixAscend_64f_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_64f_I') + use iso_c_binding + integer, intent(in), value :: len + double precision, intent(inout) :: pSrc(len) + character, intent(inout) :: pTmp(len) + end + end interface + + interface + subroutine ippsSortRadixIndexAscend_16s(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C, name='ippsSortRadixIndexAscend_16s') + use iso_c_binding + integer, intent(in), value :: len + integer*2, intent(inout) :: pSrc(len) + integer, intent(in), value :: srcStrideBytes + integer, intent(inout) :: pDstIndx(len) + character, intent(inout) :: pTmpIndx(len) end end interface interface @@ -30,7 +96,7 @@ module intel integer, intent(inout) :: pSrc(len) integer, intent(in), value :: srcStrideBytes integer, intent(inout) :: pDstIndx(len) - integer, intent(inout) :: pTmpIndx(len) + character, intent(inout) :: pTmpIndx(len) end end interface interface @@ -40,9 +106,30 @@ module intel real , intent(inout) :: pSrc(len) integer, intent(in), value :: srcStrideBytes integer, intent(inout) :: pDstIndx(len) - integer, intent(inout) :: pTmpIndx(len) + character, intent(inout) :: pTmpIndx(len) end end interface + interface + subroutine ippsSortRadixIndexAscend_64s(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C, name='ippsSortRadixIndexAscend_64s') + use iso_c_binding + integer, intent(in), value :: len + integer*8, intent(inout) :: pSrc(len) + integer, intent(in), value :: srcStrideBytes + integer, intent(inout) :: pDstIndx(len) + character, intent(inout) :: pTmpIndx(len) + end + end interface + interface + subroutine ippsSortRadixIndexAscend_64f(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C,name='ippsSortRadixIndexAscend_64f') + use iso_c_binding + integer, intent(in), value :: len + real*8 , intent(inout) :: pSrc(len) + integer, intent(in), value :: srcStrideBytes + integer, intent(inout) :: pDstIndx(len) + character, intent(inout) :: pTmpIndx(len) + end + end interface + interface subroutine ippsSortIndexAscend_32f_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_32f_I') use iso_c_binding @@ -67,4 +154,20 @@ module intel integer(4), intent(in), value :: len end end interface + interface + subroutine ippsSortIndexAscend_64s_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_64s_I') + use iso_c_binding + integer(8), intent(in) :: pSrcDst(*) + integer(4), intent(inout) :: pDstIndx(*) + integer(4), intent(in), value :: len + end + end interface + interface + subroutine ippsSortIndexAscend_16s_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_16s_I') + use iso_c_binding + integer(2), intent(in) :: pSrcDst(*) + integer(4), intent(inout) :: pDstIndx(*) + integer(4), intent(in), value :: len + end + end interface end module diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index b4e30db4..d40b7d1d 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -262,83 +262,13 @@ SUBST [ X, type ] i2 ; integer*2 ;; END_TEMPLATE + +!---------------------- INTEL IRP_IF INTEL - subroutine sort(x,iorder,isize) - use intel - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - real ,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n - call ippsSortIndexAscend_32f_I(x, iorder, isize) - iorder(:) = iorder(:)+1 - end subroutine sort - - subroutine dsort(x,iorder,isize) - use intel - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - real(8) ,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n - call ippsSortIndexAscend_64f_I(x, iorder, isize) - iorder(:) = iorder(:)+1 - end subroutine dsort - - subroutine isort(x,iorder,isize) - use intel - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - integer ,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n - integer, allocatable :: iorder1(:) - allocate(iorder1(isize*2)) - n=4 - call ippsSortRadixIndexAscend_32s(x, n, iorder, isize, iorder1) - iorder(1:isize) = iorder(1:isize)+1 - deallocate(iorder1) - call iset_order(x,iorder,isize) - end subroutine isort - - subroutine isort_noidx(x,isize) - use intel - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - integer ,intent(inout) :: x(isize) - integer, allocatable :: iorder1(:) - integer :: n - call ippsSortRadixIndexGetBufferSize(isize, 11, n) - n = n/4 - allocate(iorder1(n)) - call ippsSortRadixAscend_32s_I(x, isize, iorder1) - deallocate(iorder1) - end subroutine isort_noidx - - BEGIN_TEMPLATE subroutine $Xsort(x,iorder,isize) + use intel implicit none BEGIN_DOC ! Sort array x(isize). @@ -349,17 +279,46 @@ BEGIN_TEMPLATE $type,intent(inout) :: x(isize) integer,intent(inout) :: iorder(isize) integer :: n -! call $Xradix_sort(x,iorder,isize,-1) - call quick_$Xsort(x,iorder,isize) - end subroutine $Xsort + character, allocatable :: tmp(:) + if (isize < 2) return + call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n) + allocate(tmp(n)) + call ippsSortRadixIndexAscend_$ityp(x, $n, iorder, isize, tmp) + deallocate(tmp) + iorder(1:isize) = iorder(1:isize)+1 + call $Xset_order(x,iorder,isize) + end -SUBST [ X, type ] - i8 ; integer*8 ;; - i2 ; integer*2 ;; + subroutine $Xsort_noidx(x,isize) + use intel + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer :: n + character, allocatable :: tmp(:) + if (isize < 2) return + call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n) + allocate(tmp(n)) + call ippsSortRadixAscend_$ityp_I(x, isize, tmp) + deallocate(tmp) + end + +SUBST [ X, type, ityp, n, ippsz ] + ; real ; 32f ; 4 ; 13 ;; + d ; double precision ; 64f ; 8 ; 19;; + i ; integer ; 32s ; 4 ; 11 ;; + i8 ; integer*8 ; 64s ; 8 ; 17;; + i2 ; integer*2 ; 16s ; 2 ; 7 ;; END_TEMPLATE +!---------------------- END INTEL IRP_ELSE - +!---------------------- NON-INTEL BEGIN_TEMPLATE subroutine $Xsort(x,iorder,isize) implicit none @@ -387,50 +346,34 @@ BEGIN_TEMPLATE endif end subroutine $Xsort -SUBST [ X, type ] - ; real ;; - d ; double precision ;; -END_TEMPLATE - -BEGIN_TEMPLATE - subroutine $Xsort(x,iorder,isize) + subroutine $Xsort_noidx(x,isize) implicit none BEGIN_DOC ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. END_DOC integer,intent(in) :: isize $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n -! call $Xradix_sort(x,iorder,isize,-1) - call quick_$Xsort(x,iorder,isize) - end subroutine $Xsort + integer, allocatable :: iorder(:) + integer :: i + allocate(iorder(isize)) + do i=1,isize + iorder(i)=i + enddo + call $Xsort(x,iorder,isize) + deallocate(iorder) + end subroutine $Xsort_noidx SUBST [ X, type ] + ; real ;; + d ; double precision ;; i ; integer ;; i8 ; integer*8 ;; i2 ; integer*2 ;; END_TEMPLATE - subroutine isort_noidx(x,isize) - implicit none - BEGIN_DOC - ! Sort array x(isize). - END_DOC - integer,intent(in) :: isize - integer,intent(inout) :: x(isize) - integer, allocatable :: iorder(:) - integer :: i - allocate(iorder(isize)) - do i=1,isize - iorder(i)=i - enddo - call iradix_sort(x,iorder,isize,-1) - deallocate(iorder) - end subroutine isort_noidx IRP_ENDIF +!---------------------- END NON-INTEL + BEGIN_TEMPLATE @@ -534,9 +477,6 @@ END_TEMPLATE BEGIN_TEMPLATE recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) -IRP_IF INTEL - use intel -IRP_ENDIF implicit none BEGIN_DOC @@ -570,15 +510,6 @@ IRP_ENDIF stop endif - IRP_IF INTEL - if ( ($type == 4).and.($integer_size == 32).and.($is_big == .False.) ) then - $intel - iorder(:) = iorder(:)+1 - return - endif - IRP_ENDIF - - i1=1_$int_type i2=1_$int_type do i=1_$int_type,isize @@ -768,12 +699,12 @@ IRP_ENDIF end -SUBST [ X, type, integer_size, is_big, big, int_type, intel ] - i ; 4 ; 32 ; .False. ; ; 4 ; call ippsSortRadixIndexAscend_32s(x, 4, iorder, isize, iorder1) ;; - i8 ; 8 ; 64 ; .False. ; ; 4 ; ;; - i2 ; 2 ; 16 ; .False. ; ; 4 ; ;; - i ; 4 ; 32 ; .True. ; _big ; 8 ; ;; - i8 ; 8 ; 64 ; .True. ; _big ; 8 ; ;; +SUBST [ X, type, integer_size, is_big, big, int_type ] + i ; 4 ; 32 ; .False. ; ; 4 ;; + i8 ; 8 ; 64 ; .False. ; ; 4 ;; + i2 ; 2 ; 16 ; .False. ; ; 4 ;; + i ; 4 ; 32 ; .True. ; _big ; 8 ;; + i8 ; 8 ; 64 ; .True. ; _big ; 8 ;; END_TEMPLATE