LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-22 20:35:19 +01:00
Code Issues Releases Wiki Activity

Radix sort for negative numbers

Browse Source
This commit is contained in:
Anthony Scemama 2017-04-17 22:55:59 +02:00
parent dc2481c966
commit fd882fc0c9
2 changed files with 112 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/Determinants/spindeterminants.irp.f
View File

@ -64,9 +64,9 @@ BEGIN_TEMPLATE
integer :: i,j,k integer :: i,j,k
integer, allocatable :: iorder(:) integer, allocatable :: iorder(:)
integer(8), allocatable :: bit_tmp(:) integer*8, allocatable :: bit_tmp(:)
integer(8) :: last_key integer*8 :: last_key
integer(8), external :: spin_det_search_key integer*8, external :: spin_det_search_key
logical,allocatable :: duplicate(:) logical,allocatable :: duplicate(:)
allocate ( iorder(N_det), bit_tmp(N_det), duplicate(N_det) ) allocate ( iorder(N_det), bit_tmp(N_det), duplicate(N_det) )
@ -514,7 +514,7 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_
enddo enddo
!$OMP ENDDO !$OMP ENDDO
!$OMP END PARALLEL !$OMP END PARALLEL
call i8sort(to_sort, psi_bilinear_matrix_transp_order, N_det) call i8radix_sort(to_sort, psi_bilinear_matrix_transp_order, N_det,-1)
call iset_order(psi_bilinear_matrix_transp_rows,psi_bilinear_matrix_transp_order,N_det) call iset_order(psi_bilinear_matrix_transp_rows,psi_bilinear_matrix_transp_order,N_det)
call iset_order(psi_bilinear_matrix_transp_columns,psi_bilinear_matrix_transp_order,N_det) call iset_order(psi_bilinear_matrix_transp_columns,psi_bilinear_matrix_transp_order,N_det)
do l=1,N_states do l=1,N_states

144
src/Utils/sort.irp.f
View File

@ -186,6 +186,15 @@ BEGIN_TEMPLATE
end end
SUBST [ X, type ]
; real ;;
d ; double precision ;;
i ; integer ;;
i8 ; integer*8 ;;
i2 ; integer*2 ;;
END_TEMPLATE
BEGIN_TEMPLATE
subroutine $Xsort(x,iorder,isize) subroutine $Xsort(x,iorder,isize)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -208,6 +217,24 @@ BEGIN_TEMPLATE
SUBST [ X, type ] SUBST [ X, type ]
; real ;; ; real ;;
d ; double precision ;; d ; double precision ;;
END_TEMPLATE
BEGIN_TEMPLATE
subroutine $Xsort(x,iorder,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)
end subroutine $Xsort
SUBST [ X, type ]
i ; integer ;; i ; integer ;;
i8 ; integer*8 ;; i8 ; integer*8 ;;
i2 ; integer*2 ;; i2 ; integer*2 ;;
@ -328,33 +355,78 @@ BEGIN_TEMPLATE
integer, intent(in) :: iradix integer, intent(in) :: iradix
integer :: iradix_new integer :: iradix_new
integer*$type, allocatable :: x2(:), x1(:) integer*$type, allocatable :: x2(:), x1(:)
integer*$type :: i4 integer*$type :: i4 ! data type
integer*$int_type, allocatable :: iorder1(:),iorder2(:) integer*$int_type, allocatable :: iorder1(:),iorder2(:)
integer*$int_type :: i0, i1, i2, i3, i integer*$int_type :: i0, i1, i2, i3, i ! index type
integer, parameter :: integer_size=$octets
integer*$type :: mask integer*$type :: mask
integer :: nthreads, omp_get_num_threads integer :: err
!DIR$ ATTRIBUTES ALIGN : 128 :: iorder1,iorder2, x2, x1 !DIR$ ATTRIBUTES ALIGN : 128 :: iorder1,iorder2, x2, x1
if (iradix == -1) then if (iradix == -1) then ! Sort Positive and negative
allocate(x1(isize),iorder1(isize), x2(isize),iorder2(isize),stat=err)
if (err /= 0) then
print *, irp_here, ': Unable to allocate arrays'
stop
endif
i1=1_$int_type
i2=1_$int_type
do i=1_$int_type,isize
if (x(i) < 0_$type) then
iorder1(i1) = iorder(i)
x1(i1) = -x(i)
i1 = i1+1_$int_type
else
iorder2(i2) = iorder(i)
x2(i2) = x(i)
i2 = i2+1_$int_type
endif
enddo
i1=i1-1_$int_type
i2=i2-1_$int_type
do i=1_$int_type,i2
iorder(i1+i) = iorder2(i)
x(i1+i) = x2(i)
enddo
deallocate(x2,iorder2,stat=err)
if (err /= 0) then
print *, irp_here, ': Unable to deallocate arrays x2, iorder2'
stop
endif
if (i1 > 1_$int_type) then
call $Xradix_sort$big(x1,iorder1,i1,-2)
do i=1_$int_type,i1
x(i) = -x1(1_$int_type+i1-i)
iorder(i) = iorder1(1_$int_type+i1-i)
enddo
endif
deallocate(x1,iorder1,stat=err)
if (err /= 0) then
print *, irp_here, ': Unable to deallocate arrays x1, iorder1'
stop
endif
if (i2>1_$int_type) then
call $Xradix_sort$big(x(i1+1_$int_type),iorder(i1+1_$int_type),i2,-2)
endif
return
else if (iradix == -2) then ! Positive
! Find most significant bit ! Find most significant bit
i0 = 0_$int_type i0 = 0_$int_type
i4 = -1_$type i4 = maxval(x)
do i=1,isize iradix_new = $integer_size-1-leadz(i4)
i4 = max(i4,x(i))
enddo
i3 = int(i4,$int_type)
iradix_new = integer_size-1-leadz(i3)
mask = ibset(0_$type,iradix_new) mask = ibset(0_$type,iradix_new)
nthreads = 1
! nthreads = 1+ishft(omp_get_num_threads(),-1)
integer :: err allocate(x1(isize),iorder1(isize), x2(isize),iorder2(isize),stat=err)
allocate(x1(isize/nthreads+1),iorder1(isize/nthreads+1),x2(isize/nthreads+1),iorder2(isize/nthreads+1),stat=err)
if (err /= 0) then if (err /= 0) then
print *, irp_here, ': Unable to allocate arrays' print *, irp_here, ': Unable to allocate arrays'
stop stop
@ -363,7 +435,7 @@ BEGIN_TEMPLATE
i1=1_$int_type i1=1_$int_type
i2=1_$int_type i2=1_$int_type
do i=1,isize do i=1_$int_type,isize
if (iand(mask,x(i)) == 0_$type) then if (iand(mask,x(i)) == 0_$type) then
iorder1(i1) = iorder(i) iorder1(i1) = iorder(i)
x1(i1) = x(i) x1(i1) = x(i)
@ -377,7 +449,7 @@ BEGIN_TEMPLATE
i1=i1-1_$int_type i1=i1-1_$int_type
i2=i2-1_$int_type i2=i2-1_$int_type
do i=1,i1 do i=1_$int_type,i1
iorder(i0+i) = iorder1(i) iorder(i0+i) = iorder1(i)
x(i0+i) = x1(i) x(i0+i) = x1(i)
enddo enddo
@ -390,7 +462,7 @@ BEGIN_TEMPLATE
endif endif
do i=1,i2 do i=1_$int_type,i2
iorder(i0+i) = iorder2(i) iorder(i0+i) = iorder2(i)
x(i0+i) = x2(i) x(i0+i) = x2(i)
enddo enddo
@ -402,12 +474,12 @@ BEGIN_TEMPLATE
endif endif
if (i3>1) then if (i3>1_$int_type) then
call $Xradix_sort$big(x,iorder,i3,iradix_new-1) call $Xradix_sort$big(x,iorder,i3,iradix_new-1)
endif endif
if (isize-i3>1) then if (isize-i3>1_$int_type) then
call $Xradix_sort$big(x(i3+1),iorder(i3+1),isize-i3,iradix_new-1) call $Xradix_sort$big(x(i3+1_$int_type),iorder(i3+1_$int_type),isize-i3,iradix_new-1)
endif endif
return return
@ -429,24 +501,24 @@ BEGIN_TEMPLATE
mask = ibset(0_$type,iradix) mask = ibset(0_$type,iradix)
i0=1 i0=1_$int_type
i1=1 i1=1_$int_type
do i=1,isize do i=1_$int_type,isize
if (iand(mask,x(i)) == 0_$type) then if (iand(mask,x(i)) == 0_$type) then
iorder(i0) = iorder(i) iorder(i0) = iorder(i)
x(i0) = x(i) x(i0) = x(i)
i0 = i0+1 i0 = i0+1_$int_type
else else
iorder2(i1) = iorder(i) iorder2(i1) = iorder(i)
x2(i1) = x(i) x2(i1) = x(i)
i1 = i1+1 i1 = i1+1_$int_type
endif endif
enddo enddo
i0=i0-1 i0=i0-1_$int_type
i1=i1-1 i1=i1-1_$int_type
do i=1,i1 do i=1_$int_type,i1
iorder(i0+i) = iorder2(i) iorder(i0+i) = iorder2(i)
x(i0+i) = x2(i) x(i0+i) = x2(i)
enddo enddo
@ -463,8 +535,8 @@ BEGIN_TEMPLATE
endif endif
if (i1>1) then if (i1>1_$int_type) then
call $Xradix_sort$big(x(i0+1),iorder(i0+1),i1,iradix-1) call $Xradix_sort$big(x(i0+1_$int_type),iorder(i0+1_$int_type),i1,iradix-1)
endif endif
if (i0>1) then if (i0>1) then
call $Xradix_sort$big(x,iorder,i0,iradix-1) call $Xradix_sort$big(x,iorder,i0,iradix-1)
@ -472,11 +544,11 @@ BEGIN_TEMPLATE
end end
SUBST [ X, type, octets, is_big, big, int_type ] SUBST [ X, type, integer_size, is_big, big, int_type ]
i ; 4 ; 32 ; .False. ; ; 4 ;; i ; 4 ; 32 ; .False. ; ; 4 ;;
i8 ; 8 ; 32 ; .False. ; ; 4 ;; i8 ; 8 ; 64 ; .False. ; ; 4 ;;
i2 ; 2 ; 32 ; .False. ; ; 4 ;; i2 ; 2 ; 16 ; .False. ; ; 4 ;;
i ; 4 ; 64 ; .True. ; _big ; 8 ;; i ; 4 ; 32 ; .True. ; _big ; 8 ;;
i8 ; 8 ; 64 ; .True. ; _big ; 8 ;; i8 ; 8 ; 64 ; .True. ; _big ; 8 ;;
END_TEMPLATE END_TEMPLATE