MODULE Sorting_Ranking USE Precision USE Error_Handling USE System_Monitors USE Error_Function,ONLY:Erf,InvErf USE Random_Numbers IMPLICIT NONE PUBLIC::TestPermutation,TestRanking,DisorderPermutation,ShellInsertionRank,& RadixRank,HashRank,QuickRank,MergeRank PRIVATE LOGICAL,SAVE,PUBLIC::debug=.FALSE. INTEGER(KIND=i_dp),SAVE,PUBLIC::n_swaps=0 INTEGER,PARAMETER,PUBLIC::r_hv=r_sp INTERFACE TestRanking MODULE PROCEDURE TestRanking_1_i_sp MODULE PROCEDURE TestRanking_1_i_dp MODULE PROCEDURE TestRanking_1_r_sp MODULE PROCEDURE TestRanking_1_r_dp MODULE PROCEDURE TestRanking_2_i_sp MODULE PROCEDURE TestRanking_2_i_dp END INTERFACE INTERFACE ShellInsertionRank MODULE PROCEDURE ShellRankSingle_i_byte MODULE PROCEDURE ShellRankSingle_i_sp MODULE PROCEDURE ShellRankSingle_i_dp MODULE PROCEDURE ShellRankSingle_r_sp MODULE PROCEDURE ShellRankSingle_r_dp MODULE PROCEDURE ShellRankMultiple_i_sp MODULE PROCEDURE ShellRankMultiple_i_dp END INTERFACE INTERFACE RadixRank MODULE PROCEDURE RadixRankSingle_i_32 MODULE PROCEDURE RadixRankSingle_i_64 MODULE PROCEDURE RadixRankSingle_r_32 MODULE PROCEDURE RadixRankSingle_r_64 END INTERFACE INTERFACE HashRank MODULE PROCEDURE HashRankSingle_i_sp MODULE PROCEDURE HashRankSingle_i_dp MODULE PROCEDURE HashRankSingle_r_sp MODULE PROCEDURE HashRankSingle_r_dp END INTERFACE INTERFACE QuickRank MODULE PROCEDURE QuickRankSingle_i_byte MODULE PROCEDURE QuickRankSingle_i_sp MODULE PROCEDURE QuickRankSingle_i_dp MODULE PROCEDURE QuickRankMultiple_i_sp MODULE PROCEDURE QuickRankSingle_r_sp MODULE PROCEDURE QuickRankSingle_r_dp END INTERFACE CONTAINS FUNCTION TestPermutation(permutation)RESULT(permutation_test) IMPLICIT NONE INTEGER(KIND=i_wp),DIMENSION(:),INTENT(IN)::permutation LOGICAL::permutation_test LOGICAL(KIND=l_short),DIMENSION(SIZE(permutation))::count_array INTEGER(KIND=i_wp)::index,element count_array=.FALSE. permutation_test=.TRUE. TestDuplicates:DO index=1,SIZE(permutation) element=permutation(index) IF((1<=element).AND.(element<=SIZE(permutation)))THEN IF(.NOT.count_array(element))THEN count_array(element)=.TRUE. ELSE IF(debug)WRITE(message_print_unit,*)"Duplicate found in permutation" permutation_test=.FALSE. EXIT TestDuplicates END IF ELSE IF(debug)WRITE(message_print_unit,*)"Out of range index found in permutation" permutation_test=.FALSE. EXIT TestDuplicates END IF END DO TestDuplicates RETURN END FUNCTION TestPermutation FUNCTION TestRanking_1_i_sp(array,permutation)RESULT(ranking_test) INTEGER(KIND=i_sp),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(IN)::permutation LOGICAL::ranking_test LOGICAL::comparison_test INTEGER(KIND=i_wp)::index,key INTEGER(KIND=i_sp)::temp_compare ranking_test=.TRUE. TestMisplaced:DO index=1,SIZE(permutation)-1 comparison_test=array(permutation(index))<=array(permutation(index+1)) IF(.NOT.comparison_test)THEN IF(debug)WRITE(message_print_unit,*)"Increase in key found in ranking" ranking_test=.FALSE. EXIT TestMisplaced END IF END DO TestMisplaced RETURN END FUNCTION FUNCTION TestRanking_1_i_dp(array,permutation)RESULT(ranking_test) INTEGER(KIND=i_dp),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(IN)::permutation LOGICAL::ranking_test LOGICAL::comparison_test INTEGER(KIND=i_wp)::index,key INTEGER(KIND=i_dp)::temp_compare ranking_test=.TRUE. TestMisplaced:DO index=1,SIZE(permutation)-1 comparison_test=array(permutation(index))<=array(permutation(index+1)) IF(.NOT.comparison_test)THEN IF(debug)WRITE(message_print_unit,*)"Increase in key found in ranking" ranking_test=.FALSE. EXIT TestMisplaced END IF END DO TestMisplaced RETURN END FUNCTION FUNCTION TestRanking_1_r_sp(array,permutation)RESULT(ranking_test) REAL(KIND=r_sp),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(IN)::permutation LOGICAL::ranking_test LOGICAL::comparison_test INTEGER(KIND=i_wp)::index,key REAL(KIND=r_sp)::temp_compare ranking_test=.TRUE. TestMisplaced:DO index=1,SIZE(permutation)-1 comparison_test=array(permutation(index))<=array(permutation(index+1)) IF(.NOT.comparison_test)THEN IF(debug)WRITE(message_print_unit,*)"Increase in key found in ranking" ranking_test=.FALSE. EXIT TestMisplaced END IF END DO TestMisplaced RETURN END FUNCTION FUNCTION TestRanking_1_r_dp(array,permutation)RESULT(ranking_test) REAL(KIND=r_dp),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(IN)::permutation LOGICAL::ranking_test LOGICAL::comparison_test INTEGER(KIND=i_wp)::index,key REAL(KIND=r_dp)::temp_compare ranking_test=.TRUE. TestMisplaced:DO index=1,SIZE(permutation)-1 comparison_test=array(permutation(index))<=array(permutation(index+1)) IF(.NOT.comparison_test)THEN IF(debug)WRITE(message_print_unit,*)"Increase in key found in ranking" ranking_test=.FALSE. EXIT TestMisplaced END IF END DO TestMisplaced RETURN END FUNCTION FUNCTION TestRanking_2_i_sp(array,permutation)RESULT(ranking_test) INTEGER(KIND=i_sp),DIMENSION(:,:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(IN)::permutation LOGICAL::ranking_test LOGICAL::comparison_test INTEGER(KIND=i_wp)::index,key INTEGER(KIND=i_sp)::temp_compare ranking_test=.TRUE. TestMisplaced:DO index=1,SIZE(permutation)-1 SELECT CASE(SIZE(array,1)) CASE(1) comparison_test=(array(1,permutation(index))<=array(1,permutation(index+1))) CASE(2) temp_compare=(array(1,permutation(index))-array(1,permutation(index+1))) IF(temp_compare==0)THEN comparison_test=(array(2,permutation(index))<=array(2,permutation(index+1))) ELSE comparison_test=(temp_compare<0) END IF CASE DEFAULT comparison_test=.TRUE. CompareKey_TR:DO key=1,SIZE(array,1) temp_compare=(array(key,permutation(index))-array(key,permutation(index+1))) IF(temp_compare==0)THEN CYCLE CompareKey_TR ELSE comparison_test=(temp_compare<0) EXIT CompareKey_TR END IF END DO CompareKey_TR ENDSELECT IF(.NOT.comparison_test)THEN IF(debug)WRITE(message_print_unit,*)"Increase in key found in ranking" ranking_test=.FALSE. EXIT TestMisplaced END IF END DO TestMisplaced RETURN END FUNCTION FUNCTION TestRanking_2_i_dp(array,permutation)RESULT(ranking_test) INTEGER(KIND=i_dp),DIMENSION(:,:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(IN)::permutation LOGICAL::ranking_test LOGICAL::comparison_test INTEGER(KIND=i_wp)::index,key INTEGER(KIND=i_dp)::temp_compare ranking_test=.TRUE. TestMisplaced:DO index=1,SIZE(permutation)-1 SELECT CASE(SIZE(array,1)) CASE(1) comparison_test=(array(1,permutation(index))<=array(1,permutation(index+1))) CASE(2) temp_compare=(array(1,permutation(index))-array(1,permutation(index+1))) IF(temp_compare==0)THEN comparison_test=(array(2,permutation(index))<=array(2,permutation(index+1))) ELSE comparison_test=(temp_compare<0) END IF CASE DEFAULT comparison_test=.TRUE. CompareKey_TR:DO key=1,SIZE(array,1) temp_compare=(array(key,permutation(index))-array(key,permutation(index+1))) IF(temp_compare==0)THEN CYCLE CompareKey_TR ELSE comparison_test=(temp_compare<0) EXIT CompareKey_TR END IF END DO CompareKey_TR ENDSELECT IF(.NOT.comparison_test)THEN IF(debug)WRITE(message_print_unit,*)"Increase in key found in ranking" ranking_test=.FALSE. EXIT TestMisplaced END IF END DO TestMisplaced RETURN END FUNCTION SUBROUTINE DisorderPermutation(disorder,permutation,disorder_distribution) IMPLICIT NONE INTEGER(KIND=i_wp),DIMENSION(:),INTENT(OUT)::permutation REAL,INTENT(IN)::disorder CHARACTER(LEN=*),INTENT(IN),OPTIONAL::disorder_distribution REAL,DIMENSION(SIZE(permutation))::temp REAL::p INTEGER(KIND=i_wp)::size_permutation,index CHARACTER(LEN=1)::distribution distribution='U' IF(PRESENT(disorder_distribution))distribution=disorder_distribution(1:1) size_permutation=SIZE(permutation) p=MIN(1.0,MAX(disorder,0.0)) SELECT CASE(distribution) CASE('N','n','G','g') CALL RandomNormal(temp,mean_std=(/0.0,0.25*p*REAL(size_permutation)/)) CASE DEFAULT CALL RandomUniform(temp,range=(/0.0,p*REAL(size_permutation)/)) ENDSELECT DO index=1,size_permutation temp(index)=temp(index)+(1.0-p)*REAL(index) END DO DO index=1,size_permutation permutation(index)=index END DO CALL QuickRank(array=temp,permutation=permutation,partially_ranked=.TRUE.) END SUBROUTINE DisorderPermutation SUBROUTINE ShellRankSingle_i_byte(array,permutation,partially_ranked,method,last_gap,disorder) IMPLICIT NONE INTEGER(KIND=i_byte),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked CHARACTER(LEN=*),INTENT(IN),OPTIONAL::method INTEGER(KIND=i_wp),INTENT(IN),OPTIONAL::last_gap REAL,INTENT(IN),OPTIONAL::disorder INTEGER(KIND=i_wp)::gap,min_gap,max_gap,test_index LOGICAL::shell_sort LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status INTEGER(KIND=i_byte)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF shell_sort=.TRUE. IF(PRESENT(method))THEN IF(method(1:1)=='I'.OR.method(1:1)=='i')THEN shell_sort=.FALSE. END IF END IF IF(shell_sort)THEN IF(PRESENT(last_gap))THEN min_gap=MAX(1_i_wp,last_gap) ELSE min_gap=1_i_wp END IF IF(PRESENT(disorder))THEN max_gap=MAX(1_i_wp,MIN(INT(SIZE(array),KIND=i_wp),INT(disorder*REAL(INT(SIZE(array),KIND=i_wp)),i_wp))) ELSE max_gap=INT(SIZE(array),KIND=i_wp) END IF gap=INT(0.9102392227*LOG(2*REAL(max_gap)+1),i_wp) gap=3**gap/2 IF(debug)WRITE(message_print_unit,*)"Starting with gap:",gap SortShell_S:DO IF(gap32)n_bits=32-bit_offset offsets=0 DO element=1,INT(SIZE(array),KIND=i_wp) radix=IBITS(array(element),bit_offset,n_bits) offsets(radix)=offsets(radix)+1 END DO CheckEqual:DO radix=0,n_radices IF((offsets(radix)0))THEN EXIT CheckEqual ELSE IF(offsets(radix)==INT(SIZE(array),KIND=i_wp))THEN CYCLE BucketSort ELSE CYCLE CheckEqual END IF END DO CheckEqual IF(debug)WRITE(message_print_unit,*)"Now sorting on bits:",bit_offset,bit_offset+n_bits DO radix=1,n_radices offsets(radix)=offsets(radix)+offsets(radix-1) END DO DO index=INT(SIZE(array),KIND=i_wp),1,-1 element=permutation(index) radix=IBITS(array(element),bit_offset,n_bits) bucket_array(offsets(radix))=element offsets(radix)=offsets(radix)-1 END DO permutation=bucket_array END DO BucketSort DEALLOCATE(offsets) CALL RecordAllocation(n_elements=-n_radices-1,mold=1_i_wp) DEALLOCATE(bucket_array) CALL RecordAllocation(n_elements=-INT(SIZE(array),KIND=i_wp),mold=1_i_wp) END SUBROUTINE SUBROUTINE RadixRankSingle_i_64(array,permutation,partially_ranked,n_radix_bits,n_ignored_bits) IMPLICIT NONE INTEGER(KIND=i_64),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked INTEGER,INTENT(IN),OPTIONAL::n_radix_bits,n_ignored_bits INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::bucket_array INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::offsets INTEGER::n_bits,n_radices,n_ignored INTEGER::radix,bit_offset,bit LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status INTEGER(KIND=i_64)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF IF(PRESENT(n_radix_bits))THEN n_bits=MAX(1,MIN(BIT_SIZE(1)-1,n_radix_bits)) ELSE n_bits=8 END IF n_radices=2**n_bits-1 IF(PRESENT(n_ignored_bits))THEN n_ignored=MAX(0,n_ignored_bits) ELSE n_ignored=0 END IF ALLOCATE(offsets(0:n_radices),STAT=alloc_status) CALL RecordAllocation(n_elements=n_radices+1,mold=1_i_wp,& caller="RadixRank",alloc_status=alloc_status) ALLOCATE(bucket_array(INT(SIZE(array),KIND=i_wp)),STAT=alloc_status) CALL RecordAllocation(n_elements=INT(SIZE(array),KIND=i_wp),mold=1_i_wp,& caller="RadixRank",alloc_status=alloc_status) BucketSort:DO bit_offset=n_ignored,64-1,n_bits IF((bit_offset+n_bits)>64)n_bits=64-bit_offset offsets=0 DO element=1,INT(SIZE(array),KIND=i_wp) radix=IBITS(array(element),bit_offset,n_bits) offsets(radix)=offsets(radix)+1 END DO CheckEqual:DO radix=0,n_radices IF((offsets(radix)0))THEN EXIT CheckEqual ELSE IF(offsets(radix)==INT(SIZE(array),KIND=i_wp))THEN CYCLE BucketSort ELSE CYCLE CheckEqual END IF END DO CheckEqual IF(debug)WRITE(message_print_unit,*)"Now sorting on bits:",bit_offset,bit_offset+n_bits DO radix=1,n_radices offsets(radix)=offsets(radix)+offsets(radix-1) END DO DO index=INT(SIZE(array),KIND=i_wp),1,-1 element=permutation(index) radix=IBITS(array(element),bit_offset,n_bits) bucket_array(offsets(radix))=element offsets(radix)=offsets(radix)-1 END DO permutation=bucket_array END DO BucketSort DEALLOCATE(offsets) CALL RecordAllocation(n_elements=-n_radices-1,mold=1_i_wp) DEALLOCATE(bucket_array) CALL RecordAllocation(n_elements=-INT(SIZE(array),KIND=i_wp),mold=1_i_wp) END SUBROUTINE SUBROUTINE RadixRankSingle_r_32(array,permutation,partially_ranked,n_radix_bits,n_ignored_bits) IMPLICIT NONE REAL(KIND=r_32),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked INTEGER,INTENT(IN),OPTIONAL::n_radix_bits,n_ignored_bits INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::bucket_array INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::offsets INTEGER::n_bits,n_radices,n_ignored INTEGER::radix,bit_offset,bit INTEGER(KIND=i_32)::temp_integer_32 REAL(KIND=r_32)::temp_real_32 EQUIVALENCE(temp_integer_32,temp_real_32) LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status REAL(KIND=r_32)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF IF(PRESENT(n_radix_bits))THEN n_bits=MAX(1,MIN(BIT_SIZE(1)-1,n_radix_bits)) ELSE n_bits=8 END IF n_radices=2**n_bits-1 IF(PRESENT(n_ignored_bits))THEN n_ignored=MAX(0,n_ignored_bits) ELSE n_ignored=0 END IF ALLOCATE(offsets(0:n_radices),STAT=alloc_status) CALL RecordAllocation(n_elements=n_radices+1,mold=1_i_wp,& caller="RadixRank",alloc_status=alloc_status) ALLOCATE(bucket_array(INT(SIZE(array),KIND=i_wp)),STAT=alloc_status) CALL RecordAllocation(n_elements=INT(SIZE(array),KIND=i_wp),mold=1_i_wp,& caller="RadixRank",alloc_status=alloc_status) BucketSort:DO bit_offset=n_ignored,32-1,n_bits IF((bit_offset+n_bits)>32)n_bits=32-bit_offset offsets=0 DO element=1,INT(SIZE(array),KIND=i_wp) temp_real_32=array(element) radix=IBITS(temp_integer_32,bit_offset,n_bits) offsets(radix)=offsets(radix)+1 END DO CheckEqual:DO radix=0,n_radices IF((offsets(radix)0))THEN EXIT CheckEqual ELSE IF(offsets(radix)==INT(SIZE(array),KIND=i_wp))THEN CYCLE BucketSort ELSE CYCLE CheckEqual END IF END DO CheckEqual IF(debug)WRITE(message_print_unit,*)"Now sorting on bits:",bit_offset,bit_offset+n_bits DO radix=1,n_radices offsets(radix)=offsets(radix)+offsets(radix-1) END DO DO index=INT(SIZE(array),KIND=i_wp),1,-1 element=permutation(index) temp_real_32=array(element) radix=IBITS(temp_integer_32,bit_offset,n_bits) bucket_array(offsets(radix))=element offsets(radix)=offsets(radix)-1 END DO permutation=bucket_array END DO BucketSort DEALLOCATE(offsets) CALL RecordAllocation(n_elements=-n_radices-1,mold=1_i_wp) DEALLOCATE(bucket_array) CALL RecordAllocation(n_elements=-INT(SIZE(array),KIND=i_wp),mold=1_i_wp) END SUBROUTINE SUBROUTINE RadixRankSingle_r_64(array,permutation,partially_ranked,n_radix_bits,n_ignored_bits) IMPLICIT NONE REAL(KIND=r_64),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked INTEGER,INTENT(IN),OPTIONAL::n_radix_bits,n_ignored_bits INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::bucket_array INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::offsets INTEGER::n_bits,n_radices,n_ignored INTEGER::radix,bit_offset,bit INTEGER(KIND=i_64)::temp_integer_64 REAL(KIND=r_64)::temp_real_64 EQUIVALENCE(temp_integer_64,temp_real_64) LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status REAL(KIND=r_64)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF IF(PRESENT(n_radix_bits))THEN n_bits=MAX(1,MIN(BIT_SIZE(1)-1,n_radix_bits)) ELSE n_bits=8 END IF n_radices=2**n_bits-1 IF(PRESENT(n_ignored_bits))THEN n_ignored=MAX(0,n_ignored_bits) ELSE n_ignored=0 END IF ALLOCATE(offsets(0:n_radices),STAT=alloc_status) CALL RecordAllocation(n_elements=n_radices+1,mold=1_i_wp,& caller="RadixRank",alloc_status=alloc_status) ALLOCATE(bucket_array(INT(SIZE(array),KIND=i_wp)),STAT=alloc_status) CALL RecordAllocation(n_elements=INT(SIZE(array),KIND=i_wp),mold=1_i_wp,& caller="RadixRank",alloc_status=alloc_status) BucketSort:DO bit_offset=n_ignored,64-1,n_bits IF((bit_offset+n_bits)>64)n_bits=64-bit_offset offsets=0 DO element=1,INT(SIZE(array),KIND=i_wp) temp_real_64=array(element) radix=IBITS(temp_integer_64,bit_offset,n_bits) offsets(radix)=offsets(radix)+1 END DO CheckEqual:DO radix=0,n_radices IF((offsets(radix)0))THEN EXIT CheckEqual ELSE IF(offsets(radix)==INT(SIZE(array),KIND=i_wp))THEN CYCLE BucketSort ELSE CYCLE CheckEqual END IF END DO CheckEqual IF(debug)WRITE(message_print_unit,*)"Now sorting on bits:",bit_offset,bit_offset+n_bits DO radix=1,n_radices offsets(radix)=offsets(radix)+offsets(radix-1) END DO DO index=INT(SIZE(array),KIND=i_wp),1,-1 element=permutation(index) temp_real_64=array(element) radix=IBITS(temp_integer_64,bit_offset,n_bits) bucket_array(offsets(radix))=element offsets(radix)=offsets(radix)-1 END DO permutation=bucket_array END DO BucketSort DEALLOCATE(offsets) CALL RecordAllocation(n_elements=-n_radices-1,mold=1_i_wp) DEALLOCATE(bucket_array) CALL RecordAllocation(n_elements=-INT(SIZE(array),KIND=i_wp),mold=1_i_wp) END SUBROUTINE SUBROUTINE HashRankSingle_i_sp(array,permutation,partially_ranked,mean_value,standard_deviation,number_passes,memory_factor,& &uniformity,hash_function,distribution,store_hash_values,use_hash_pointers,full_ranking) IMPLICIT NONE INTEGER(KIND=i_sp),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked INTEGER(KIND=i_sp),INTENT(IN),OPTIONAL::mean_value,standard_deviation INTEGER,INTENT(IN),OPTIONAL::number_passes REAL,INTENT(IN),OPTIONAL::memory_factor REAL,INTENT(IN),OPTIONAL::uniformity OPTIONAL::hash_function INTERFACE FUNCTION hash_function(value)RESULT(hash_factor) USE Precision INTEGER(KIND=i_sp),INTENT(IN)::value REAL(KIND=r_sp)::hash_factor END FUNCTION hash_function END INTERFACE CHARACTER(LEN=*),INTENT(IN),OPTIONAL::distribution LOGICAL,INTENT(IN),OPTIONAL::store_hash_values,use_hash_pointers,full_ranking INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::hash_array,hash_pointers,hash_values INTEGER(KIND=i_wp)::size_array,size_hash INTEGER(KIND=i_wp)::hash_value,hash_guess,hash_range(2) INTEGER(KIND=i_wp)::hash_index,n_placed,n_hashed,n_hashed_max,n_front,n_back INTEGER::n_passes,pass REAL(KIND=r_hv)::mem_factor,extra_space,hash_factor REAL(KIND=r_hv)::mean,std,inv_width,total LOGICAL::uniform_distribution,store_hash,one_pass,use_pointers,finish_ranking LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status INTEGER(KIND=i_sp)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF IF(PRESENT(number_passes).AND.(number_passes>0))THEN n_passes=number_passes ELSE n_passes=1 END IF IF(n_passes==1)one_pass=.TRUE. IF(PRESENT(memory_factor))THEN mem_factor=REAL(MAX(1.0,memory_factor),KIND=r_hv) ELSE mem_factor=2.0_r_hv END IF IF(PRESENT(uniformity))THEN extra_space=REAL(MAX(EPSILON(0.0),MIN(1.0,uniformity)),KIND=r_hv) ELSE extra_space=0.9_r_hv END IF uniform_distribution=.FALSE. IF(PRESENT(distribution))THEN IF(distribution(1:1)=='U'.OR.distribution(1:1)=='u')& uniform_distribution=.TRUE. END IF store_hash=.FALSE. IF(PRESENT(store_hash_values))THEN IF(store_hash_values.AND.(n_passes>1))store_hash=.TRUE. END IF use_pointers=.FALSE. IF(PRESENT(use_hash_pointers))use_pointers=use_hash_pointers finish_ranking=.TRUE. IF(PRESENT(full_ranking))finish_ranking=full_ranking size_array=INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) size_hash=INT(mem_factor*REAL(size_array/n_passes,KIND=r_hv)/extra_space,i_wp)+1 n_hashed_max=INT(REAL(size_array/n_passes,KIND=r_hv)/extra_space,i_wp) ALLOCATE(hash_array(size_hash),STAT=alloc_status) CALL RecordAllocation(n_elements=size_hash,mold=1_i_wp,& caller="HashRank",alloc_status=alloc_status) IF(use_pointers.OR.(.NOT.one_pass))THEN ALLOCATE(hash_pointers(size_hash),STAT=alloc_status) CALL RecordAllocation(n_elements=size_hash,mold=1_i_wp,& caller="HashRank",alloc_status=alloc_status) END IF IF(.NOT.PRESENT(hash_function))THEN IF(PRESENT(mean_value))THEN mean=REAL(mean_value,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,size_array total=total+REAL(array(element),KIND=r_hv) END DO mean=total/REAL(size_array,KIND=r_hv) END IF IF(PRESENT(standard_deviation))THEN std=REAL(standard_deviation,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,size_array total=total+(REAL(array(element),KIND=r_hv)-mean)**2 END DO std=SQRT(total/REAL(size_array,KIND=r_hv)) END IF IF(uniform_distribution)THEN inv_width=1.0_r_hv/SQRT(3.0_r_hv)/std ELSE inv_width=1.0_r_hv/SQRT(2.0_r_hv)/std END IF END IF IF(store_hash)THEN ALLOCATE(hash_values(size_array),STAT=alloc_status) CALL RecordAllocation(n_elements=size_hash,mold=1_i_wp,& caller="HashRank",alloc_status=alloc_status) DO element=1,size_array IF(PRESENT(hash_function))THEN hash_factor=hash_function(array(element)) ELSE IF(uniform_distribution)THEN hash_factor=0.5_r_hv*((REAL(array(element),KIND=r_hv)-mean)*inv_width+1.0_r_hv) ELSE hash_factor=0.5_r_hv*(Erf((REAL(array(element),KIND=r_hv)-mean)*inv_width)+1.0_r_hv) END IF END IF hash_value=INT(extra_space*hash_factor*REAL(n_passes*size_hash,KIND=r_hv),i_wp)+1 hash_values(element)=hash_value END DO END IF n_placed=0 MakePass:DO pass=1,n_passes IF(debug)THEN WRITE(message_print_unit,"(A,I2)")" Now making pass #",pass END IF hash_array=0 IF(use_pointers)hash_pointers=0 hash_range=size_hash*(/(pass-1),pass/) hash_range=INT(extra_space*REAL(hash_range,KIND=r_hv),i_wp) n_hashed=0 TraverseArray:DO index=n_placed+1,size_array element=permutation(index) IF(store_hash)THEN hash_value=hash_values(element) ELSE IF(PRESENT(hash_function))THEN hash_factor=hash_function(array(element)) ELSE IF(uniform_distribution)THEN hash_factor=0.5_r_hv*((REAL(array(element),KIND=r_hv)-mean)*inv_width+1.0_r_hv) ELSE hash_factor=0.5_r_hv*(Erf((REAL(array(element),KIND=r_hv)-mean)*inv_width)+1.0_r_hv) END IF END IF hash_value=INT(extra_space*hash_factor*REAL(n_passes*size_hash,KIND=r_hv),i_wp)+1 END IF HashEntry:IF(hash_value=n_hashed_max)EXIT TraverseArray n_hashed=n_hashed+1 hash_guess=hash_value UsePointers:IF(.NOT.use_pointers)THEN SearchForFreeLocation:DO IF(hash_array(hash_value)==0)EXIT SearchForFreeLocation IF(hash_value==size_hash)hash_value=1 hash_value=hash_value+1 END DO SearchForFreeLocation ELSE FindFreeLocation:DO IF(hash_pointers(hash_value)==0)EXIT FindFreeLocation hash_value=hash_pointers(hash_value) END DO FindFreeLocation IF(hash_value/=size_hash)THEN hash_pointers(hash_value)=hash_value+1 hash_pointers(hash_guess)=hash_value+1 ELSE hash_pointers(hash_value)=1 hash_pointers(hash_guess)=1 END IF END IF UsePointers OnePass:IF(one_pass)THEN hash_array(hash_value)=permutation(index) ELSE hash_array(hash_value)=index END IF OnePass END IF HashEntry END DO TraverseArray IF(debug)THEN WRITE(message_print_unit,*)"Successfully hashed ",n_hashed," elements, out of ",size_array/n_passes IF(SIZE(hash_array)<=20)THEN WRITE(message_print_unit,"(A)")"Hash array:" WRITE(message_print_unit,"(20G5.2)")hash_array END IF IF(use_pointers)THEN IF(SIZE(hash_pointers)<=20)THEN WRITE(message_print_unit,"(A)")"Hash pointers:" WRITE(message_print_unit,"(20G5.2)")hash_pointers END IF END IF END IF IF(one_pass)THEN CompressToFront:DO hash_index=1,size_hash IF(hash_array(hash_index)/=0)THEN n_placed=n_placed+1 permutation(n_placed)=hash_array(hash_index) END IF END DO CompressToFront ELSE n_front=0 MoveToBack:DO hash_index=1,size_hash IF(hash_array(hash_index)/=0)THEN n_front=n_front+1 index=hash_array(hash_index) hash_array(n_front)=permutation(index) hash_pointers(n_front)=index FindStoreLocation:DO IF(index>=n_placed+n_front)THEN EXIT FindStoreLocation ELSE index=hash_pointers(index-n_placed) END IF END DO FindStoreLocation permutation(index)=permutation(n_placed+n_front) END IF END DO MoveToBack MoveToFront:DO hash_index=1,n_hashed n_placed=n_placed+1 permutation(n_placed)=hash_array(hash_index) END DO MoveToFront END IF IF(debug)THEN WRITE(message_print_unit,*)"Successfully placed ",n_placed," out of ",(pass*size_array)/n_passes IF(SIZE(permutation)<=20)THEN WRITE(message_print_unit,"(A)")"Permutation array:" WRITE(message_print_unit,"(20G5.2)")permutation END IF END IF END DO MakePass DEALLOCATE(hash_array) CALL RecordAllocation(n_elements=-size_hash,mold=1_i_wp) IF(ALLOCATED(hash_pointers))THEN DEALLOCATE(hash_pointers) CALL RecordAllocation(n_elements=-size_hash,mold=1_i_wp) END IF IF(ALLOCATED(hash_values))THEN DEALLOCATE(hash_values) CALL RecordAllocation(n_elements=-size_array,mold=1_i_wp) END IF IF(finish_ranking)CALL ShellInsertionRank(array=array,permutation=permutation,& partially_ranked=.TRUE.,method="Insertion") END SUBROUTINE SUBROUTINE HashRankSingle_i_dp(array,permutation,partially_ranked,mean_value,standard_deviation,number_passes,memory_factor,& &uniformity,hash_function,distribution,store_hash_values,use_hash_pointers,full_ranking) IMPLICIT NONE INTEGER(KIND=i_dp),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked INTEGER(KIND=i_dp),INTENT(IN),OPTIONAL::mean_value,standard_deviation INTEGER,INTENT(IN),OPTIONAL::number_passes REAL,INTENT(IN),OPTIONAL::memory_factor REAL,INTENT(IN),OPTIONAL::uniformity OPTIONAL::hash_function INTERFACE FUNCTION hash_function(value)RESULT(hash_factor) USE Precision INTEGER(KIND=i_dp),INTENT(IN)::value REAL(KIND=r_sp)::hash_factor END FUNCTION hash_function END INTERFACE CHARACTER(LEN=*),INTENT(IN),OPTIONAL::distribution LOGICAL,INTENT(IN),OPTIONAL::store_hash_values,use_hash_pointers,full_ranking INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::hash_array,hash_pointers,hash_values INTEGER(KIND=i_wp)::size_array,size_hash INTEGER(KIND=i_wp)::hash_value,hash_guess,hash_range(2) INTEGER(KIND=i_wp)::hash_index,n_placed,n_hashed,n_hashed_max,n_front,n_back INTEGER::n_passes,pass REAL(KIND=r_hv)::mem_factor,extra_space,hash_factor REAL(KIND=r_hv)::mean,std,inv_width,total LOGICAL::uniform_distribution,store_hash,one_pass,use_pointers,finish_ranking LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status INTEGER(KIND=i_dp)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF IF(PRESENT(number_passes).AND.(number_passes>0))THEN n_passes=number_passes ELSE n_passes=1 END IF IF(n_passes==1)one_pass=.TRUE. IF(PRESENT(memory_factor))THEN mem_factor=REAL(MAX(1.0,memory_factor),KIND=r_hv) ELSE mem_factor=2.0_r_hv END IF IF(PRESENT(uniformity))THEN extra_space=REAL(MAX(EPSILON(0.0),MIN(1.0,uniformity)),KIND=r_hv) ELSE extra_space=0.9_r_hv END IF uniform_distribution=.FALSE. IF(PRESENT(distribution))THEN IF(distribution(1:1)=='U'.OR.distribution(1:1)=='u')& uniform_distribution=.TRUE. END IF store_hash=.FALSE. IF(PRESENT(store_hash_values))THEN IF(store_hash_values.AND.(n_passes>1))store_hash=.TRUE. END IF use_pointers=.FALSE. IF(PRESENT(use_hash_pointers))use_pointers=use_hash_pointers finish_ranking=.TRUE. IF(PRESENT(full_ranking))finish_ranking=full_ranking size_array=INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) size_hash=INT(mem_factor*REAL(size_array/n_passes,KIND=r_hv)/extra_space,i_wp)+1 n_hashed_max=INT(REAL(size_array/n_passes,KIND=r_hv)/extra_space,i_wp) ALLOCATE(hash_array(size_hash),STAT=alloc_status) CALL RecordAllocation(n_elements=size_hash,mold=1_i_wp,& caller="HashRank",alloc_status=alloc_status) IF(use_pointers.OR.(.NOT.one_pass))THEN ALLOCATE(hash_pointers(size_hash),STAT=alloc_status) CALL RecordAllocation(n_elements=size_hash,mold=1_i_wp,& caller="HashRank",alloc_status=alloc_status) END IF IF(.NOT.PRESENT(hash_function))THEN IF(PRESENT(mean_value))THEN mean=REAL(mean_value,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,size_array total=total+REAL(array(element),KIND=r_hv) END DO mean=total/REAL(size_array,KIND=r_hv) END IF IF(PRESENT(standard_deviation))THEN std=REAL(standard_deviation,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,size_array total=total+(REAL(array(element),KIND=r_hv)-mean)**2 END DO std=SQRT(total/REAL(size_array,KIND=r_hv)) END IF IF(uniform_distribution)THEN inv_width=1.0_r_hv/SQRT(3.0_r_hv)/std ELSE inv_width=1.0_r_hv/SQRT(2.0_r_hv)/std END IF END IF IF(store_hash)THEN ALLOCATE(hash_values(size_array),STAT=alloc_status) CALL RecordAllocation(n_elements=size_hash,mold=1_i_wp,& caller="HashRank",alloc_status=alloc_status) DO element=1,size_array IF(PRESENT(hash_function))THEN hash_factor=hash_function(array(element)) ELSE IF(uniform_distribution)THEN hash_factor=0.5_r_hv*((REAL(array(element),KIND=r_hv)-mean)*inv_width+1.0_r_hv) ELSE hash_factor=0.5_r_hv*(Erf((REAL(array(element),KIND=r_hv)-mean)*inv_width)+1.0_r_hv) END IF END IF hash_value=INT(extra_space*hash_factor*REAL(n_passes*size_hash,KIND=r_hv),i_wp)+1 hash_values(element)=hash_value END DO END IF n_placed=0 MakePass:DO pass=1,n_passes IF(debug)THEN WRITE(message_print_unit,"(A,I2)")" Now making pass #",pass END IF hash_array=0 IF(use_pointers)hash_pointers=0 hash_range=size_hash*(/(pass-1),pass/) hash_range=INT(extra_space*REAL(hash_range,KIND=r_hv),i_wp) n_hashed=0 TraverseArray:DO index=n_placed+1,size_array element=permutation(index) IF(store_hash)THEN hash_value=hash_values(element) ELSE IF(PRESENT(hash_function))THEN hash_factor=hash_function(array(element)) ELSE IF(uniform_distribution)THEN hash_factor=0.5_r_hv*((REAL(array(element),KIND=r_hv)-mean)*inv_width+1.0_r_hv) ELSE hash_factor=0.5_r_hv*(Erf((REAL(array(element),KIND=r_hv)-mean)*inv_width)+1.0_r_hv) END IF END IF hash_value=INT(extra_space*hash_factor*REAL(n_passes*size_hash,KIND=r_hv),i_wp)+1 END IF HashEntry:IF(hash_value=n_hashed_max)EXIT TraverseArray n_hashed=n_hashed+1 hash_guess=hash_value UsePointers:IF(.NOT.use_pointers)THEN SearchForFreeLocation:DO IF(hash_array(hash_value)==0)EXIT SearchForFreeLocation IF(hash_value==size_hash)hash_value=1 hash_value=hash_value+1 END DO SearchForFreeLocation ELSE FindFreeLocation:DO IF(hash_pointers(hash_value)==0)EXIT FindFreeLocation hash_value=hash_pointers(hash_value) END DO FindFreeLocation IF(hash_value/=size_hash)THEN hash_pointers(hash_value)=hash_value+1 hash_pointers(hash_guess)=hash_value+1 ELSE hash_pointers(hash_value)=1 hash_pointers(hash_guess)=1 END IF END IF UsePointers OnePass:IF(one_pass)THEN hash_array(hash_value)=permutation(index) ELSE hash_array(hash_value)=index END IF OnePass END IF HashEntry END DO TraverseArray IF(debug)THEN WRITE(message_print_unit,*)"Successfully hashed ",n_hashed," elements, out of ",size_array/n_passes IF(SIZE(hash_array)<=20)THEN WRITE(message_print_unit,"(A)")"Hash array:" WRITE(message_print_unit,"(20G5.2)")hash_array END IF IF(use_pointers)THEN IF(SIZE(hash_pointers)<=20)THEN WRITE(message_print_unit,"(A)")"Hash pointers:" WRITE(message_print_unit,"(20G5.2)")hash_pointers END IF END IF END IF IF(one_pass)THEN CompressToFront:DO hash_index=1,size_hash IF(hash_array(hash_index)/=0)THEN n_placed=n_placed+1 permutation(n_placed)=hash_array(hash_index) END IF END DO CompressToFront ELSE n_front=0 MoveToBack:DO hash_index=1,size_hash IF(hash_array(hash_index)/=0)THEN n_front=n_front+1 index=hash_array(hash_index) hash_array(n_front)=permutation(index) hash_pointers(n_front)=index FindStoreLocation:DO IF(index>=n_placed+n_front)THEN EXIT FindStoreLocation ELSE index=hash_pointers(index-n_placed) END IF END DO FindStoreLocation permutation(index)=permutation(n_placed+n_front) END IF END DO MoveToBack MoveToFront:DO hash_index=1,n_hashed n_placed=n_placed+1 permutation(n_placed)=hash_array(hash_index) END DO MoveToFront END IF IF(debug)THEN WRITE(message_print_unit,*)"Successfully placed ",n_placed," out of ",(pass*size_array)/n_passes IF(SIZE(permutation)<=20)THEN WRITE(message_print_unit,"(A)")"Permutation array:" WRITE(message_print_unit,"(20G5.2)")permutation END IF END IF END DO MakePass DEALLOCATE(hash_array) CALL RecordAllocation(n_elements=-size_hash,mold=1_i_wp) IF(ALLOCATED(hash_pointers))THEN DEALLOCATE(hash_pointers) CALL RecordAllocation(n_elements=-size_hash,mold=1_i_wp) END IF IF(ALLOCATED(hash_values))THEN DEALLOCATE(hash_values) CALL RecordAllocation(n_elements=-size_array,mold=1_i_wp) END IF IF(finish_ranking)CALL ShellInsertionRank(array=array,permutation=permutation,& partially_ranked=.TRUE.,method="Insertion") END SUBROUTINE SUBROUTINE HashRankSingle_r_sp(array,permutation,partially_ranked,mean_value,standard_deviation,number_passes,memory_factor,& &uniformity,hash_function,distribution,store_hash_values,use_hash_pointers,full_ranking) IMPLICIT NONE REAL(KIND=r_sp),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked REAL(KIND=r_sp),INTENT(IN),OPTIONAL::mean_value,standard_deviation INTEGER,INTENT(IN),OPTIONAL::number_passes REAL,INTENT(IN),OPTIONAL::memory_factor REAL,INTENT(IN),OPTIONAL::uniformity OPTIONAL::hash_function INTERFACE FUNCTION hash_function(value)RESULT(hash_factor) USE Precision REAL(KIND=r_sp),INTENT(IN)::value REAL(KIND=r_sp)::hash_factor END FUNCTION hash_function END INTERFACE CHARACTER(LEN=*),INTENT(IN),OPTIONAL::distribution LOGICAL,INTENT(IN),OPTIONAL::store_hash_values,use_hash_pointers,full_ranking INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::hash_array,hash_pointers,hash_values INTEGER(KIND=i_wp)::size_array,size_hash INTEGER(KIND=i_wp)::hash_value,hash_guess,hash_range(2) INTEGER(KIND=i_wp)::hash_index,n_placed,n_hashed,n_hashed_max,n_front,n_back INTEGER::n_passes,pass REAL(KIND=r_hv)::mem_factor,extra_space,hash_factor REAL(KIND=r_hv)::mean,std,inv_width,total LOGICAL::uniform_distribution,store_hash,one_pass,use_pointers,finish_ranking LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status REAL(KIND=r_sp)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF IF(PRESENT(number_passes).AND.(number_passes>0))THEN n_passes=number_passes ELSE n_passes=1 END IF IF(n_passes==1)one_pass=.TRUE. IF(PRESENT(memory_factor))THEN mem_factor=REAL(MAX(1.0,memory_factor),KIND=r_hv) ELSE mem_factor=2.0_r_hv END IF IF(PRESENT(uniformity))THEN extra_space=REAL(MAX(EPSILON(0.0),MIN(1.0,uniformity)),KIND=r_hv) ELSE extra_space=0.9_r_hv END IF uniform_distribution=.FALSE. IF(PRESENT(distribution))THEN IF(distribution(1:1)=='U'.OR.distribution(1:1)=='u')& uniform_distribution=.TRUE. END IF store_hash=.FALSE. IF(PRESENT(store_hash_values))THEN IF(store_hash_values.AND.(n_passes>1))store_hash=.TRUE. END IF use_pointers=.FALSE. IF(PRESENT(use_hash_pointers))use_pointers=use_hash_pointers finish_ranking=.TRUE. IF(PRESENT(full_ranking))finish_ranking=full_ranking size_array=INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) size_hash=INT(mem_factor*REAL(size_array/n_passes,KIND=r_hv)/extra_space,i_wp)+1 n_hashed_max=INT(REAL(size_array/n_passes,KIND=r_hv)/extra_space,i_wp) ALLOCATE(hash_array(size_hash),STAT=alloc_status) CALL RecordAllocation(n_elements=size_hash,mold=1_i_wp,& caller="HashRank",alloc_status=alloc_status) IF(use_pointers.OR.(.NOT.one_pass))THEN ALLOCATE(hash_pointers(size_hash),STAT=alloc_status) CALL RecordAllocation(n_elements=size_hash,mold=1_i_wp,& caller="HashRank",alloc_status=alloc_status) END IF IF(.NOT.PRESENT(hash_function))THEN IF(PRESENT(mean_value))THEN mean=REAL(mean_value,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,size_array total=total+REAL(array(element),KIND=r_hv) END DO mean=total/REAL(size_array,KIND=r_hv) END IF IF(PRESENT(standard_deviation))THEN std=REAL(standard_deviation,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,size_array total=total+(REAL(array(element),KIND=r_hv)-mean)**2 END DO std=SQRT(total/REAL(size_array,KIND=r_hv)) END IF IF(uniform_distribution)THEN inv_width=1.0_r_hv/SQRT(3.0_r_hv)/std ELSE inv_width=1.0_r_hv/SQRT(2.0_r_hv)/std END IF END IF IF(store_hash)THEN ALLOCATE(hash_values(size_array),STAT=alloc_status) CALL RecordAllocation(n_elements=size_hash,mold=1_i_wp,& caller="HashRank",alloc_status=alloc_status) DO element=1,size_array IF(PRESENT(hash_function))THEN hash_factor=hash_function(array(element)) ELSE IF(uniform_distribution)THEN hash_factor=0.5_r_hv*((REAL(array(element),KIND=r_hv)-mean)*inv_width+1.0_r_hv) ELSE hash_factor=0.5_r_hv*(Erf((REAL(array(element),KIND=r_hv)-mean)*inv_width)+1.0_r_hv) END IF END IF hash_value=INT(extra_space*hash_factor*REAL(n_passes*size_hash,KIND=r_hv),i_wp)+1 hash_values(element)=hash_value END DO END IF n_placed=0 MakePass:DO pass=1,n_passes IF(debug)THEN WRITE(message_print_unit,"(A,I2)")" Now making pass #",pass END IF hash_array=0 IF(use_pointers)hash_pointers=0 hash_range=size_hash*(/(pass-1),pass/) hash_range=INT(extra_space*REAL(hash_range,KIND=r_hv),i_wp) n_hashed=0 TraverseArray:DO index=n_placed+1,size_array element=permutation(index) IF(store_hash)THEN hash_value=hash_values(element) ELSE IF(PRESENT(hash_function))THEN hash_factor=hash_function(array(element)) ELSE IF(uniform_distribution)THEN hash_factor=0.5_r_hv*((REAL(array(element),KIND=r_hv)-mean)*inv_width+1.0_r_hv) ELSE hash_factor=0.5_r_hv*(Erf((REAL(array(element),KIND=r_hv)-mean)*inv_width)+1.0_r_hv) END IF END IF hash_value=INT(extra_space*hash_factor*REAL(n_passes*size_hash,KIND=r_hv),i_wp)+1 END IF HashEntry:IF(hash_value=n_hashed_max)EXIT TraverseArray n_hashed=n_hashed+1 hash_guess=hash_value UsePointers:IF(.NOT.use_pointers)THEN SearchForFreeLocation:DO IF(hash_array(hash_value)==0)EXIT SearchForFreeLocation IF(hash_value==size_hash)hash_value=1 hash_value=hash_value+1 END DO SearchForFreeLocation ELSE FindFreeLocation:DO IF(hash_pointers(hash_value)==0)EXIT FindFreeLocation hash_value=hash_pointers(hash_value) END DO FindFreeLocation IF(hash_value/=size_hash)THEN hash_pointers(hash_value)=hash_value+1 hash_pointers(hash_guess)=hash_value+1 ELSE hash_pointers(hash_value)=1 hash_pointers(hash_guess)=1 END IF END IF UsePointers OnePass:IF(one_pass)THEN hash_array(hash_value)=permutation(index) ELSE hash_array(hash_value)=index END IF OnePass END IF HashEntry END DO TraverseArray IF(debug)THEN WRITE(message_print_unit,*)"Successfully hashed ",n_hashed," elements, out of ",size_array/n_passes IF(SIZE(hash_array)<=20)THEN WRITE(message_print_unit,"(A)")"Hash array:" WRITE(message_print_unit,"(20G5.2)")hash_array END IF IF(use_pointers)THEN IF(SIZE(hash_pointers)<=20)THEN WRITE(message_print_unit,"(A)")"Hash pointers:" WRITE(message_print_unit,"(20G5.2)")hash_pointers END IF END IF END IF IF(one_pass)THEN CompressToFront:DO hash_index=1,size_hash IF(hash_array(hash_index)/=0)THEN n_placed=n_placed+1 permutation(n_placed)=hash_array(hash_index) END IF END DO CompressToFront ELSE n_front=0 MoveToBack:DO hash_index=1,size_hash IF(hash_array(hash_index)/=0)THEN n_front=n_front+1 index=hash_array(hash_index) hash_array(n_front)=permutation(index) hash_pointers(n_front)=index FindStoreLocation:DO IF(index>=n_placed+n_front)THEN EXIT FindStoreLocation ELSE index=hash_pointers(index-n_placed) END IF END DO FindStoreLocation permutation(index)=permutation(n_placed+n_front) END IF END DO MoveToBack MoveToFront:DO hash_index=1,n_hashed n_placed=n_placed+1 permutation(n_placed)=hash_array(hash_index) END DO MoveToFront END IF IF(debug)THEN WRITE(message_print_unit,*)"Successfully placed ",n_placed," out of ",(pass*size_array)/n_passes IF(SIZE(permutation)<=20)THEN WRITE(message_print_unit,"(A)")"Permutation array:" WRITE(message_print_unit,"(20G5.2)")permutation END IF END IF END DO MakePass DEALLOCATE(hash_array) CALL RecordAllocation(n_elements=-size_hash,mold=1_i_wp) IF(ALLOCATED(hash_pointers))THEN DEALLOCATE(hash_pointers) CALL RecordAllocation(n_elements=-size_hash,mold=1_i_wp) END IF IF(ALLOCATED(hash_values))THEN DEALLOCATE(hash_values) CALL RecordAllocation(n_elements=-size_array,mold=1_i_wp) END IF IF(finish_ranking)CALL ShellInsertionRank(array=array,permutation=permutation,& partially_ranked=.TRUE.,method="Insertion") END SUBROUTINE SUBROUTINE HashRankSingle_r_dp(array,permutation,partially_ranked,mean_value,standard_deviation,number_passes,memory_factor,& &uniformity,hash_function,distribution,store_hash_values,use_hash_pointers,full_ranking) IMPLICIT NONE REAL(KIND=r_dp),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked REAL(KIND=r_dp),INTENT(IN),OPTIONAL::mean_value,standard_deviation INTEGER,INTENT(IN),OPTIONAL::number_passes REAL,INTENT(IN),OPTIONAL::memory_factor REAL,INTENT(IN),OPTIONAL::uniformity OPTIONAL::hash_function INTERFACE FUNCTION hash_function(value)RESULT(hash_factor) USE Precision REAL(KIND=r_dp),INTENT(IN)::value REAL(KIND=r_sp)::hash_factor END FUNCTION hash_function END INTERFACE CHARACTER(LEN=*),INTENT(IN),OPTIONAL::distribution LOGICAL,INTENT(IN),OPTIONAL::store_hash_values,use_hash_pointers,full_ranking INTEGER(KIND=i_wp),DIMENSION(:),ALLOCATABLE::hash_array,hash_pointers,hash_values INTEGER(KIND=i_wp)::size_array,size_hash INTEGER(KIND=i_wp)::hash_value,hash_guess,hash_range(2) INTEGER(KIND=i_wp)::hash_index,n_placed,n_hashed,n_hashed_max,n_front,n_back INTEGER::n_passes,pass REAL(KIND=r_hv)::mem_factor,extra_space,hash_factor REAL(KIND=r_hv)::mean,std,inv_width,total LOGICAL::uniform_distribution,store_hash,one_pass,use_pointers,finish_ranking LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status REAL(KIND=r_dp)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF IF(PRESENT(number_passes).AND.(number_passes>0))THEN n_passes=number_passes ELSE n_passes=1 END IF IF(n_passes==1)one_pass=.TRUE. IF(PRESENT(memory_factor))THEN mem_factor=REAL(MAX(1.0,memory_factor),KIND=r_hv) ELSE mem_factor=2.0_r_hv END IF IF(PRESENT(uniformity))THEN extra_space=REAL(MAX(EPSILON(0.0),MIN(1.0,uniformity)),KIND=r_hv) ELSE extra_space=0.9_r_hv END IF uniform_distribution=.FALSE. IF(PRESENT(distribution))THEN IF(distribution(1:1)=='U'.OR.distribution(1:1)=='u')& uniform_distribution=.TRUE. END IF store_hash=.FALSE. IF(PRESENT(store_hash_values))THEN IF(store_hash_values.AND.(n_passes>1))store_hash=.TRUE. END IF use_pointers=.FALSE. IF(PRESENT(use_hash_pointers))use_pointers=use_hash_pointers finish_ranking=.TRUE. IF(PRESENT(full_ranking))finish_ranking=full_ranking size_array=INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) size_hash=INT(mem_factor*REAL(size_array/n_passes,KIND=r_hv)/extra_space,i_wp)+1 n_hashed_max=INT(REAL(size_array/n_passes,KIND=r_hv)/extra_space,i_wp) ALLOCATE(hash_array(size_hash),STAT=alloc_status) CALL RecordAllocation(n_elements=size_hash,mold=1_i_wp,& caller="HashRank",alloc_status=alloc_status) IF(use_pointers.OR.(.NOT.one_pass))THEN ALLOCATE(hash_pointers(size_hash),STAT=alloc_status) CALL RecordAllocation(n_elements=size_hash,mold=1_i_wp,& caller="HashRank",alloc_status=alloc_status) END IF IF(.NOT.PRESENT(hash_function))THEN IF(PRESENT(mean_value))THEN mean=REAL(mean_value,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,size_array total=total+REAL(array(element),KIND=r_hv) END DO mean=total/REAL(size_array,KIND=r_hv) END IF IF(PRESENT(standard_deviation))THEN std=REAL(standard_deviation,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,size_array total=total+(REAL(array(element),KIND=r_hv)-mean)**2 END DO std=SQRT(total/REAL(size_array,KIND=r_hv)) END IF IF(uniform_distribution)THEN inv_width=1.0_r_hv/SQRT(3.0_r_hv)/std ELSE inv_width=1.0_r_hv/SQRT(2.0_r_hv)/std END IF END IF IF(store_hash)THEN ALLOCATE(hash_values(size_array),STAT=alloc_status) CALL RecordAllocation(n_elements=size_hash,mold=1_i_wp,& caller="HashRank",alloc_status=alloc_status) DO element=1,size_array IF(PRESENT(hash_function))THEN hash_factor=hash_function(array(element)) ELSE IF(uniform_distribution)THEN hash_factor=0.5_r_hv*((REAL(array(element),KIND=r_hv)-mean)*inv_width+1.0_r_hv) ELSE hash_factor=0.5_r_hv*(Erf((REAL(array(element),KIND=r_hv)-mean)*inv_width)+1.0_r_hv) END IF END IF hash_value=INT(extra_space*hash_factor*REAL(n_passes*size_hash,KIND=r_hv),i_wp)+1 hash_values(element)=hash_value END DO END IF n_placed=0 MakePass:DO pass=1,n_passes IF(debug)THEN WRITE(message_print_unit,"(A,I2)")" Now making pass #",pass END IF hash_array=0 IF(use_pointers)hash_pointers=0 hash_range=size_hash*(/(pass-1),pass/) hash_range=INT(extra_space*REAL(hash_range,KIND=r_hv),i_wp) n_hashed=0 TraverseArray:DO index=n_placed+1,size_array element=permutation(index) IF(store_hash)THEN hash_value=hash_values(element) ELSE IF(PRESENT(hash_function))THEN hash_factor=hash_function(array(element)) ELSE IF(uniform_distribution)THEN hash_factor=0.5_r_hv*((REAL(array(element),KIND=r_hv)-mean)*inv_width+1.0_r_hv) ELSE hash_factor=0.5_r_hv*(Erf((REAL(array(element),KIND=r_hv)-mean)*inv_width)+1.0_r_hv) END IF END IF hash_value=INT(extra_space*hash_factor*REAL(n_passes*size_hash,KIND=r_hv),i_wp)+1 END IF HashEntry:IF(hash_value=n_hashed_max)EXIT TraverseArray n_hashed=n_hashed+1 hash_guess=hash_value UsePointers:IF(.NOT.use_pointers)THEN SearchForFreeLocation:DO IF(hash_array(hash_value)==0)EXIT SearchForFreeLocation IF(hash_value==size_hash)hash_value=1 hash_value=hash_value+1 END DO SearchForFreeLocation ELSE FindFreeLocation:DO IF(hash_pointers(hash_value)==0)EXIT FindFreeLocation hash_value=hash_pointers(hash_value) END DO FindFreeLocation IF(hash_value/=size_hash)THEN hash_pointers(hash_value)=hash_value+1 hash_pointers(hash_guess)=hash_value+1 ELSE hash_pointers(hash_value)=1 hash_pointers(hash_guess)=1 END IF END IF UsePointers OnePass:IF(one_pass)THEN hash_array(hash_value)=permutation(index) ELSE hash_array(hash_value)=index END IF OnePass END IF HashEntry END DO TraverseArray IF(debug)THEN WRITE(message_print_unit,*)"Successfully hashed ",n_hashed," elements, out of ",size_array/n_passes IF(SIZE(hash_array)<=20)THEN WRITE(message_print_unit,"(A)")"Hash array:" WRITE(message_print_unit,"(20G5.2)")hash_array END IF IF(use_pointers)THEN IF(SIZE(hash_pointers)<=20)THEN WRITE(message_print_unit,"(A)")"Hash pointers:" WRITE(message_print_unit,"(20G5.2)")hash_pointers END IF END IF END IF IF(one_pass)THEN CompressToFront:DO hash_index=1,size_hash IF(hash_array(hash_index)/=0)THEN n_placed=n_placed+1 permutation(n_placed)=hash_array(hash_index) END IF END DO CompressToFront ELSE n_front=0 MoveToBack:DO hash_index=1,size_hash IF(hash_array(hash_index)/=0)THEN n_front=n_front+1 index=hash_array(hash_index) hash_array(n_front)=permutation(index) hash_pointers(n_front)=index FindStoreLocation:DO IF(index>=n_placed+n_front)THEN EXIT FindStoreLocation ELSE index=hash_pointers(index-n_placed) END IF END DO FindStoreLocation permutation(index)=permutation(n_placed+n_front) END IF END DO MoveToBack MoveToFront:DO hash_index=1,n_hashed n_placed=n_placed+1 permutation(n_placed)=hash_array(hash_index) END DO MoveToFront END IF IF(debug)THEN WRITE(message_print_unit,*)"Successfully placed ",n_placed," out of ",(pass*size_array)/n_passes IF(SIZE(permutation)<=20)THEN WRITE(message_print_unit,"(A)")"Permutation array:" WRITE(message_print_unit,"(20G5.2)")permutation END IF END IF END DO MakePass DEALLOCATE(hash_array) CALL RecordAllocation(n_elements=-size_hash,mold=1_i_wp) IF(ALLOCATED(hash_pointers))THEN DEALLOCATE(hash_pointers) CALL RecordAllocation(n_elements=-size_hash,mold=1_i_wp) END IF IF(ALLOCATED(hash_values))THEN DEALLOCATE(hash_values) CALL RecordAllocation(n_elements=-size_array,mold=1_i_wp) END IF IF(finish_ranking)CALL ShellInsertionRank(array=array,permutation=permutation,& partially_ranked=.TRUE.,method="Insertion") END SUBROUTINE SUBROUTINE QuickRankSingle_i_byte(array,permutation,partially_ranked,mean_value,standard_deviation,pivot_selection,cutoff_si& &ze,pivot_function) IMPLICIT NONE INTEGER(KIND=i_byte),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked INTEGER(KIND=i_byte),INTENT(IN),OPTIONAL::mean_value,standard_deviation CHARACTER(LEN=*),INTENT(IN),OPTIONAL::pivot_selection INTEGER(KIND=i_wp),INTENT(IN),OPTIONAL::cutoff_size OPTIONAL::pivot_function INTERFACE FUNCTION pivot_function(lower_bound,upper_bound)RESULT(pivot_value) USE Precision REAL(KIND=r_sp),INTENT(IN)::lower_bound,upper_bound INTEGER(KIND=i_byte)::pivot_value END FUNCTION pivot_function END INTERFACE INTEGER(KIND=i_wp)::cutoff INTEGER::recursion_depth,max_recursion_depth CHARACTER(LEN=1)::pivot_method REAL(KIND=r_hv)::mean,std,total INTEGER(KIND=i_byte)::min_val,max_val LOGICAL::uniform_distribution LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status INTEGER(KIND=i_byte)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF pivot_method='R' IF(PRESENT(pivot_selection))pivot_method=pivot_selection IF(PRESENT(cutoff_size))THEN cutoff=MAX(1_i_wp,cutoff_size) ELSE cutoff=25 END IF recursion_depth=0; max_recursion_depth=0 SELECT CASE(pivot_method) CASE('U','u','G','g','N','n') IF(.NOT.PRESENT(pivot_function))THEN IF(PRESENT(mean_value))THEN mean=REAL(mean_value,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) total=total+REAL(array(element),KIND=r_hv) END DO mean=total/REAL(INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp),KIND=r_hv) END IF IF(PRESENT(standard_deviation))THEN std=REAL(standard_deviation,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) total=total+(REAL(array(element),KIND=r_hv)-mean)**2 END DO std=SQRT(total/REAL(INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp),KIND=r_hv)) END IF END IF IF(pivot_method=='U'.OR.pivot_method=='u')THEN CALL QuickSortRecursive_RD_i_byte(1_i_wp,INT(SIZE(array),KIND=i_wp),& range=(/mean-std*SQRT(3.0_r_hv),mean+std*SQRT(3.0_r_hv)/),distribution='U') ELSE CALL QuickSortRecursive_RD_i_byte(1_i_wp,INT(SIZE(array),KIND=i_wp),& range=(/mean-std*10.0_r_hv,mean+std*10.0/),distribution='N') END IF CASE DEFAULT CALL QuickSortRecursive_PE_i_byte(1_i_wp,INT(SIZE(array),KIND=i_wp),(INT(SIZE(array),KIND=i_wp)+1)/2) ENDSELECT IF(debug)WRITE(message_print_unit,*)"Recursion depth=",max_recursion_depth IF(.NOT.PRESENT(cutoff_size))THEN CALL ShellInsertionRank(array=array,permutation=permutation,& partially_ranked=.TRUE.,method="Insertion") END IF CONTAINS RECURSIVE SUBROUTINE QuickSortRecursive_RD_i_byte(starting,ending,range,distribution) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::starting,ending REAL(KIND=r_hv),DIMENSION(2),INTENT(IN)::range CHARACTER(LEN=1),INTENT(IN)::distribution INTEGER(KIND=i_wp)::low,high INTEGER(KIND=i_byte)::pivot_key recursion_depth=recursion_depth+1 IF((ending-starting)=range(2))THEN pivot_key=(range(1)+range(2))/2 END IF END IF ELSE pivot_key=pivot_function(range(1),range(2)) END IF END IF NibbleArray_RD:DO NibbleLeft_RD:DO IF(low>high)EXIT NibbleArray_RD comparison_test=(array(permutation(low))<=pivot_key) IF(.NOT.comparison_test)EXIT NibbleLeft_RD low=low+1 END DO NibbleLeft_RD NibbleRight_RD:DO IF(highending).OR.(high=MAX(25_i_wp,cutoff))THEN IF(debug)WRITE(message_print_unit,*)"Recalculating range for section of size: ",ending-starting+1 min_val=HUGE(array); max_val=-HUGE(array) DO index=starting,ending temp_compare=array(permutation(index)) IF(temp_comparemax_val)max_val=temp_compare END DO CALL QuickSortRecursive_RD_i_byte(starting,ending,& range=REAL((/min_val,max_val/),KIND=r_hv),distribution=distribution) ELSE CALL QuickSortRecursive_PE_i_byte(starting,ending,(starting+ending)/2) END IF ELSE IF(low-starting>=cutoff)CALL QuickSortRecursive_RD_i_byte(starting,low-1,& range=(/range(1),REAL(pivot_key,KIND=r_hv)/),distribution=distribution) IF(ending-high>=cutoff)CALL QuickSortRecursive_RD_i_byte(high+1,ending,& range=(/REAL(pivot_key,KIND=r_hv),range(2)/),distribution=distribution) END IF max_recursion_depth=MAX(max_recursion_depth,recursion_depth) recursion_depth=recursion_depth-1 RETURN END SUBROUTINE RECURSIVE SUBROUTINE QuickSortRecursive_PE_i_byte(starting,ending,pivot_index) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::starting,ending INTEGER(KIND=i_wp),INTENT(IN)::pivot_index INTEGER(KIND=i_wp)::new_pivot_index INTEGER(KIND=i_wp)::low,high INTEGER(KIND=i_byte)::pivot_key recursion_depth=recursion_depth+1 IF((ending-starting)high)EXIT NibbleArray_PE comparison_test=(array(permutation(low))<=pivot_key) IF(.NOT.comparison_test)EXIT NibbleLeft_PE low=low+1 END DO NibbleLeft_PE NibbleRight_PE:DO IF(high=cutoff)THEN CALL RandomUniform(new_pivot_index,range=(/starting,low-1/)) CALL QuickSortRecursive_PE_i_byte(starting,low-1,new_pivot_index) END IF IF(ending-high>=cutoff)THEN CALL RandomUniform(new_pivot_index,range=(/high+1,ending/)) CALL QuickSortRecursive_PE_i_byte(high+1,ending,new_pivot_index) END IF max_recursion_depth=MAX(max_recursion_depth,recursion_depth) recursion_depth=recursion_depth-1 RETURN END SUBROUTINE END SUBROUTINE SUBROUTINE QuickRankSingle_i_sp(array,permutation,partially_ranked,mean_value,standard_deviation,pivot_selection,cutoff_size& &,pivot_function) IMPLICIT NONE INTEGER(KIND=i_sp),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked INTEGER(KIND=i_sp),INTENT(IN),OPTIONAL::mean_value,standard_deviation CHARACTER(LEN=*),INTENT(IN),OPTIONAL::pivot_selection INTEGER(KIND=i_wp),INTENT(IN),OPTIONAL::cutoff_size OPTIONAL::pivot_function INTERFACE FUNCTION pivot_function(lower_bound,upper_bound)RESULT(pivot_value) USE Precision REAL(KIND=r_sp),INTENT(IN)::lower_bound,upper_bound INTEGER(KIND=i_sp)::pivot_value END FUNCTION pivot_function END INTERFACE INTEGER(KIND=i_wp)::cutoff INTEGER::recursion_depth,max_recursion_depth CHARACTER(LEN=1)::pivot_method REAL(KIND=r_hv)::mean,std,total INTEGER(KIND=i_sp)::min_val,max_val LOGICAL::uniform_distribution LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status INTEGER(KIND=i_sp)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF pivot_method='R' IF(PRESENT(pivot_selection))pivot_method=pivot_selection IF(PRESENT(cutoff_size))THEN cutoff=MAX(1_i_wp,cutoff_size) ELSE cutoff=25 END IF recursion_depth=0; max_recursion_depth=0 SELECT CASE(pivot_method) CASE('U','u','G','g','N','n') IF(.NOT.PRESENT(pivot_function))THEN IF(PRESENT(mean_value))THEN mean=REAL(mean_value,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) total=total+REAL(array(element),KIND=r_hv) END DO mean=total/REAL(INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp),KIND=r_hv) END IF IF(PRESENT(standard_deviation))THEN std=REAL(standard_deviation,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) total=total+(REAL(array(element),KIND=r_hv)-mean)**2 END DO std=SQRT(total/REAL(INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp),KIND=r_hv)) END IF END IF IF(pivot_method=='U'.OR.pivot_method=='u')THEN CALL QuickSortRecursive_RD_i_sp(1_i_wp,INT(SIZE(array),KIND=i_wp),& range=(/mean-std*SQRT(3.0_r_hv),mean+std*SQRT(3.0_r_hv)/),distribution='U') ELSE CALL QuickSortRecursive_RD_i_sp(1_i_wp,INT(SIZE(array),KIND=i_wp),& range=(/mean-std*10.0_r_hv,mean+std*10.0/),distribution='N') END IF CASE DEFAULT CALL QuickSortRecursive_PE_i_sp(1_i_wp,INT(SIZE(array),KIND=i_wp),(INT(SIZE(array),KIND=i_wp)+1)/2) ENDSELECT IF(debug)WRITE(message_print_unit,*)"Recursion depth=",max_recursion_depth IF(.NOT.PRESENT(cutoff_size))THEN CALL ShellInsertionRank(array=array,permutation=permutation,& partially_ranked=.TRUE.,method="Insertion") END IF CONTAINS RECURSIVE SUBROUTINE QuickSortRecursive_RD_i_sp(starting,ending,range,distribution) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::starting,ending REAL(KIND=r_hv),DIMENSION(2),INTENT(IN)::range CHARACTER(LEN=1),INTENT(IN)::distribution INTEGER(KIND=i_wp)::low,high INTEGER(KIND=i_sp)::pivot_key recursion_depth=recursion_depth+1 IF((ending-starting)=range(2))THEN pivot_key=(range(1)+range(2))/2 END IF END IF ELSE pivot_key=pivot_function(range(1),range(2)) END IF END IF NibbleArray_RD:DO NibbleLeft_RD:DO IF(low>high)EXIT NibbleArray_RD comparison_test=(array(permutation(low))<=pivot_key) IF(.NOT.comparison_test)EXIT NibbleLeft_RD low=low+1 END DO NibbleLeft_RD NibbleRight_RD:DO IF(highending).OR.(high=MAX(25_i_wp,cutoff))THEN IF(debug)WRITE(message_print_unit,*)"Recalculating range for section of size: ",ending-starting+1 min_val=HUGE(array); max_val=-HUGE(array) DO index=starting,ending temp_compare=array(permutation(index)) IF(temp_comparemax_val)max_val=temp_compare END DO CALL QuickSortRecursive_RD_i_sp(starting,ending,& range=REAL((/min_val,max_val/),KIND=r_hv),distribution=distribution) ELSE CALL QuickSortRecursive_PE_i_sp(starting,ending,(starting+ending)/2) END IF ELSE IF(low-starting>=cutoff)CALL QuickSortRecursive_RD_i_sp(starting,low-1,& range=(/range(1),REAL(pivot_key,KIND=r_hv)/),distribution=distribution) IF(ending-high>=cutoff)CALL QuickSortRecursive_RD_i_sp(high+1,ending,& range=(/REAL(pivot_key,KIND=r_hv),range(2)/),distribution=distribution) END IF max_recursion_depth=MAX(max_recursion_depth,recursion_depth) recursion_depth=recursion_depth-1 RETURN END SUBROUTINE RECURSIVE SUBROUTINE QuickSortRecursive_PE_i_sp(starting,ending,pivot_index) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::starting,ending INTEGER(KIND=i_wp),INTENT(IN)::pivot_index INTEGER(KIND=i_wp)::new_pivot_index INTEGER(KIND=i_wp)::low,high INTEGER(KIND=i_sp)::pivot_key recursion_depth=recursion_depth+1 IF((ending-starting)high)EXIT NibbleArray_PE comparison_test=(array(permutation(low))<=pivot_key) IF(.NOT.comparison_test)EXIT NibbleLeft_PE low=low+1 END DO NibbleLeft_PE NibbleRight_PE:DO IF(high=cutoff)THEN CALL RandomUniform(new_pivot_index,range=(/starting,low-1/)) CALL QuickSortRecursive_PE_i_sp(starting,low-1,new_pivot_index) END IF IF(ending-high>=cutoff)THEN CALL RandomUniform(new_pivot_index,range=(/high+1,ending/)) CALL QuickSortRecursive_PE_i_sp(high+1,ending,new_pivot_index) END IF max_recursion_depth=MAX(max_recursion_depth,recursion_depth) recursion_depth=recursion_depth-1 RETURN END SUBROUTINE END SUBROUTINE SUBROUTINE QuickRankSingle_i_dp(array,permutation,partially_ranked,mean_value,standard_deviation,pivot_selection,cutoff_size& &,pivot_function) IMPLICIT NONE INTEGER(KIND=i_dp),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked INTEGER(KIND=i_dp),INTENT(IN),OPTIONAL::mean_value,standard_deviation CHARACTER(LEN=*),INTENT(IN),OPTIONAL::pivot_selection INTEGER(KIND=i_wp),INTENT(IN),OPTIONAL::cutoff_size OPTIONAL::pivot_function INTERFACE FUNCTION pivot_function(lower_bound,upper_bound)RESULT(pivot_value) USE Precision REAL(KIND=r_sp),INTENT(IN)::lower_bound,upper_bound INTEGER(KIND=i_dp)::pivot_value END FUNCTION pivot_function END INTERFACE INTEGER(KIND=i_wp)::cutoff INTEGER::recursion_depth,max_recursion_depth CHARACTER(LEN=1)::pivot_method REAL(KIND=r_hv)::mean,std,total INTEGER(KIND=i_dp)::min_val,max_val LOGICAL::uniform_distribution LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status INTEGER(KIND=i_dp)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF pivot_method='R' IF(PRESENT(pivot_selection))pivot_method=pivot_selection IF(PRESENT(cutoff_size))THEN cutoff=MAX(1_i_wp,cutoff_size) ELSE cutoff=25 END IF recursion_depth=0; max_recursion_depth=0 SELECT CASE(pivot_method) CASE('U','u','G','g','N','n') IF(.NOT.PRESENT(pivot_function))THEN IF(PRESENT(mean_value))THEN mean=REAL(mean_value,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) total=total+REAL(array(element),KIND=r_hv) END DO mean=total/REAL(INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp),KIND=r_hv) END IF IF(PRESENT(standard_deviation))THEN std=REAL(standard_deviation,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) total=total+(REAL(array(element),KIND=r_hv)-mean)**2 END DO std=SQRT(total/REAL(INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp),KIND=r_hv)) END IF END IF IF(pivot_method=='U'.OR.pivot_method=='u')THEN CALL QuickSortRecursive_RD_i_dp(1_i_wp,INT(SIZE(array),KIND=i_wp),& range=(/mean-std*SQRT(3.0_r_hv),mean+std*SQRT(3.0_r_hv)/),distribution='U') ELSE CALL QuickSortRecursive_RD_i_dp(1_i_wp,INT(SIZE(array),KIND=i_wp),& range=(/mean-std*10.0_r_hv,mean+std*10.0/),distribution='N') END IF CASE DEFAULT CALL QuickSortRecursive_PE_i_dp(1_i_wp,INT(SIZE(array),KIND=i_wp),(INT(SIZE(array),KIND=i_wp)+1)/2) ENDSELECT IF(debug)WRITE(message_print_unit,*)"Recursion depth=",max_recursion_depth IF(.NOT.PRESENT(cutoff_size))THEN CALL ShellInsertionRank(array=array,permutation=permutation,& partially_ranked=.TRUE.,method="Insertion") END IF CONTAINS RECURSIVE SUBROUTINE QuickSortRecursive_RD_i_dp(starting,ending,range,distribution) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::starting,ending REAL(KIND=r_hv),DIMENSION(2),INTENT(IN)::range CHARACTER(LEN=1),INTENT(IN)::distribution INTEGER(KIND=i_wp)::low,high INTEGER(KIND=i_dp)::pivot_key recursion_depth=recursion_depth+1 IF((ending-starting)=range(2))THEN pivot_key=(range(1)+range(2))/2 END IF END IF ELSE pivot_key=pivot_function(range(1),range(2)) END IF END IF NibbleArray_RD:DO NibbleLeft_RD:DO IF(low>high)EXIT NibbleArray_RD comparison_test=(array(permutation(low))<=pivot_key) IF(.NOT.comparison_test)EXIT NibbleLeft_RD low=low+1 END DO NibbleLeft_RD NibbleRight_RD:DO IF(highending).OR.(high=MAX(25_i_wp,cutoff))THEN IF(debug)WRITE(message_print_unit,*)"Recalculating range for section of size: ",ending-starting+1 min_val=HUGE(array); max_val=-HUGE(array) DO index=starting,ending temp_compare=array(permutation(index)) IF(temp_comparemax_val)max_val=temp_compare END DO CALL QuickSortRecursive_RD_i_dp(starting,ending,& range=REAL((/min_val,max_val/),KIND=r_hv),distribution=distribution) ELSE CALL QuickSortRecursive_PE_i_dp(starting,ending,(starting+ending)/2) END IF ELSE IF(low-starting>=cutoff)CALL QuickSortRecursive_RD_i_dp(starting,low-1,& range=(/range(1),REAL(pivot_key,KIND=r_hv)/),distribution=distribution) IF(ending-high>=cutoff)CALL QuickSortRecursive_RD_i_dp(high+1,ending,& range=(/REAL(pivot_key,KIND=r_hv),range(2)/),distribution=distribution) END IF max_recursion_depth=MAX(max_recursion_depth,recursion_depth) recursion_depth=recursion_depth-1 RETURN END SUBROUTINE RECURSIVE SUBROUTINE QuickSortRecursive_PE_i_dp(starting,ending,pivot_index) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::starting,ending INTEGER(KIND=i_wp),INTENT(IN)::pivot_index INTEGER(KIND=i_wp)::new_pivot_index INTEGER(KIND=i_wp)::low,high INTEGER(KIND=i_dp)::pivot_key recursion_depth=recursion_depth+1 IF((ending-starting)high)EXIT NibbleArray_PE comparison_test=(array(permutation(low))<=pivot_key) IF(.NOT.comparison_test)EXIT NibbleLeft_PE low=low+1 END DO NibbleLeft_PE NibbleRight_PE:DO IF(high=cutoff)THEN CALL RandomUniform(new_pivot_index,range=(/starting,low-1/)) CALL QuickSortRecursive_PE_i_dp(starting,low-1,new_pivot_index) END IF IF(ending-high>=cutoff)THEN CALL RandomUniform(new_pivot_index,range=(/high+1,ending/)) CALL QuickSortRecursive_PE_i_dp(high+1,ending,new_pivot_index) END IF max_recursion_depth=MAX(max_recursion_depth,recursion_depth) recursion_depth=recursion_depth-1 RETURN END SUBROUTINE END SUBROUTINE SUBROUTINE QuickRankMultiple_i_sp(array,permutation,partially_ranked,pivot_selection,cutoff_size) IMPLICIT NONE INTEGER(KIND=i_sp),DIMENSION(:,:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked CHARACTER(LEN=1),INTENT(IN),OPTIONAL::pivot_selection INTEGER(KIND=i_wp),INTENT(IN),OPTIONAL::cutoff_size INTEGER(KIND=i_wp)::cutoff INTEGER::recursion_depth,max_recursion_depth CHARACTER(LEN=1)::pivot_method LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element,key INTEGER::alloc_status INTEGER(KIND=i_sp)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array,DIM=2),KIND=i_wp) permutation(index)=index END DO END IF pivot_method='R' IF(PRESENT(pivot_selection))pivot_method=pivot_selection IF(PRESENT(cutoff_size))THEN cutoff=MAX(1_i_wp,cutoff_size) ELSE cutoff=25 END IF recursion_depth=0; max_recursion_depth=0 IF((pivot_method/='R').AND.(pivot_method/='r'))THEN IF(debug)WRITE(message_print_unit,*)"Requested pivot selection method not available for multiple keys! Switching to random p& &ivot selection..." END IF CALL QuickSortRecursive_PE_i_sp(1_i_wp,INT(SIZE(array,DIM=2),KIND=i_wp),(INT(SIZE(array,DIM=2),KIND=i_wp)+1)/2) IF(debug)WRITE(message_print_unit,*)"Recursion depth=",max_recursion_depth IF(.NOT.PRESENT(cutoff_size))THEN CALL ShellInsertionRank(array=array,permutation=permutation,& partially_ranked=.TRUE.,method="Insertion") END IF CONTAINS RECURSIVE SUBROUTINE QuickSortRecursive_PE_i_sp(starting,ending,pivot_index) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::starting,ending INTEGER(KIND=i_wp),INTENT(IN)::pivot_index INTEGER(KIND=i_wp)::new_pivot_index INTEGER(KIND=i_wp)::low,high INTEGER(KIND=i_sp)::pivot_key(INT(SIZE(array,1),KIND=i_wp)) recursion_depth=recursion_depth+1 IF((ending-starting)high)EXIT NibbleArray_PE SELECT CASE(SIZE(array,1)) CASE(1) comparison_test=(array(1,permutation(low))<=pivot_key(1)) CASE(2) temp_compare=(array(1,permutation(low))-pivot_key(1)) IF(temp_compare==0)THEN comparison_test=(array(2,permutation(low))<=pivot_key(2)) ELSE comparison_test=(temp_compare<0) END IF CASE DEFAULT comparison_test=.TRUE. CompareKey_PE_NL:DO key=1,SIZE(array,1) temp_compare=(array(key,permutation(low))-pivot_key(key)) IF(temp_compare==0)THEN CYCLE CompareKey_PE_NL ELSE comparison_test=(temp_compare<0) EXIT CompareKey_PE_NL END IF END DO CompareKey_PE_NL ENDSELECT IF(.NOT.comparison_test)EXIT NibbleLeft_PE low=low+1 END DO NibbleLeft_PE NibbleRight_PE:DO IF(high=cutoff)THEN CALL RandomUniform(new_pivot_index,range=(/starting,low-1/)) CALL QuickSortRecursive_PE_i_sp(starting,low-1,new_pivot_index) END IF IF(ending-high>=cutoff)THEN CALL RandomUniform(new_pivot_index,range=(/high+1,ending/)) CALL QuickSortRecursive_PE_i_sp(high+1,ending,new_pivot_index) END IF max_recursion_depth=MAX(max_recursion_depth,recursion_depth) recursion_depth=recursion_depth-1 RETURN END SUBROUTINE END SUBROUTINE SUBROUTINE QuickRankSingle_r_sp(array,permutation,partially_ranked,mean_value,standard_deviation,pivot_selection,cutoff_size& &,pivot_function) IMPLICIT NONE REAL(KIND=r_sp),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked REAL(KIND=r_sp),INTENT(IN),OPTIONAL::mean_value,standard_deviation CHARACTER(LEN=*),INTENT(IN),OPTIONAL::pivot_selection INTEGER(KIND=i_wp),INTENT(IN),OPTIONAL::cutoff_size OPTIONAL::pivot_function INTERFACE FUNCTION pivot_function(lower_bound,upper_bound)RESULT(pivot_value) USE Precision REAL(KIND=r_sp),INTENT(IN)::lower_bound,upper_bound REAL(KIND=r_sp)::pivot_value END FUNCTION pivot_function END INTERFACE INTEGER(KIND=i_wp)::cutoff INTEGER::recursion_depth,max_recursion_depth CHARACTER(LEN=1)::pivot_method REAL(KIND=r_hv)::mean,std,total REAL(KIND=r_sp)::min_val,max_val LOGICAL::uniform_distribution LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status REAL(KIND=r_sp)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF pivot_method='R' IF(PRESENT(pivot_selection))pivot_method=pivot_selection IF(PRESENT(cutoff_size))THEN cutoff=MAX(1_i_wp,cutoff_size) ELSE cutoff=25 END IF recursion_depth=0; max_recursion_depth=0 SELECT CASE(pivot_method) CASE('U','u','G','g','N','n') IF(.NOT.PRESENT(pivot_function))THEN IF(PRESENT(mean_value))THEN mean=REAL(mean_value,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) total=total+REAL(array(element),KIND=r_hv) END DO mean=total/REAL(INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp),KIND=r_hv) END IF IF(PRESENT(standard_deviation))THEN std=REAL(standard_deviation,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) total=total+(REAL(array(element),KIND=r_hv)-mean)**2 END DO std=SQRT(total/REAL(INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp),KIND=r_hv)) END IF END IF IF(pivot_method=='U'.OR.pivot_method=='u')THEN CALL QuickSortRecursive_RD_r_sp(1_i_wp,INT(SIZE(array),KIND=i_wp),& range=(/mean-std*SQRT(3.0_r_hv),mean+std*SQRT(3.0_r_hv)/),distribution='U') ELSE CALL QuickSortRecursive_RD_r_sp(1_i_wp,INT(SIZE(array),KIND=i_wp),& range=(/mean-std*10.0_r_hv,mean+std*10.0/),distribution='N') END IF CASE DEFAULT CALL QuickSortRecursive_PE_r_sp(1_i_wp,INT(SIZE(array),KIND=i_wp),(INT(SIZE(array),KIND=i_wp)+1)/2) ENDSELECT IF(debug)WRITE(message_print_unit,*)"Recursion depth=",max_recursion_depth IF(.NOT.PRESENT(cutoff_size))THEN CALL ShellInsertionRank(array=array,permutation=permutation,& partially_ranked=.TRUE.,method="Insertion") END IF CONTAINS RECURSIVE SUBROUTINE QuickSortRecursive_RD_r_sp(starting,ending,range,distribution) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::starting,ending REAL(KIND=r_hv),DIMENSION(2),INTENT(IN)::range CHARACTER(LEN=1),INTENT(IN)::distribution INTEGER(KIND=i_wp)::low,high REAL(KIND=r_sp)::pivot_key recursion_depth=recursion_depth+1 IF((ending-starting)=range(2))THEN pivot_key=(range(1)+range(2))/2 END IF END IF ELSE pivot_key=pivot_function(range(1),range(2)) END IF END IF NibbleArray_RD:DO NibbleLeft_RD:DO IF(low>high)EXIT NibbleArray_RD comparison_test=(array(permutation(low))<=pivot_key) IF(.NOT.comparison_test)EXIT NibbleLeft_RD low=low+1 END DO NibbleLeft_RD NibbleRight_RD:DO IF(highending).OR.(high=MAX(25_i_wp,cutoff))THEN IF(debug)WRITE(message_print_unit,*)"Recalculating range for section of size: ",ending-starting+1 min_val=HUGE(array); max_val=-HUGE(array) DO index=starting,ending temp_compare=array(permutation(index)) IF(temp_comparemax_val)max_val=temp_compare END DO CALL QuickSortRecursive_RD_r_sp(starting,ending,& range=REAL((/min_val,max_val/),KIND=r_hv),distribution=distribution) ELSE CALL QuickSortRecursive_PE_r_sp(starting,ending,(starting+ending)/2) END IF ELSE IF(low-starting>=cutoff)CALL QuickSortRecursive_RD_r_sp(starting,low-1,& range=(/range(1),REAL(pivot_key,KIND=r_hv)/),distribution=distribution) IF(ending-high>=cutoff)CALL QuickSortRecursive_RD_r_sp(high+1,ending,& range=(/REAL(pivot_key,KIND=r_hv),range(2)/),distribution=distribution) END IF max_recursion_depth=MAX(max_recursion_depth,recursion_depth) recursion_depth=recursion_depth-1 RETURN END SUBROUTINE RECURSIVE SUBROUTINE QuickSortRecursive_PE_r_sp(starting,ending,pivot_index) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::starting,ending INTEGER(KIND=i_wp),INTENT(IN)::pivot_index INTEGER(KIND=i_wp)::new_pivot_index INTEGER(KIND=i_wp)::low,high REAL(KIND=r_sp)::pivot_key recursion_depth=recursion_depth+1 IF((ending-starting)high)EXIT NibbleArray_PE comparison_test=(array(permutation(low))<=pivot_key) IF(.NOT.comparison_test)EXIT NibbleLeft_PE low=low+1 END DO NibbleLeft_PE NibbleRight_PE:DO IF(high=cutoff)THEN CALL RandomUniform(new_pivot_index,range=(/starting,low-1/)) CALL QuickSortRecursive_PE_r_sp(starting,low-1,new_pivot_index) END IF IF(ending-high>=cutoff)THEN CALL RandomUniform(new_pivot_index,range=(/high+1,ending/)) CALL QuickSortRecursive_PE_r_sp(high+1,ending,new_pivot_index) END IF max_recursion_depth=MAX(max_recursion_depth,recursion_depth) recursion_depth=recursion_depth-1 RETURN END SUBROUTINE END SUBROUTINE SUBROUTINE QuickRankSingle_r_dp(array,permutation,partially_ranked,mean_value,standard_deviation,pivot_selection,cutoff_size& &,pivot_function) IMPLICIT NONE REAL(KIND=r_dp),DIMENSION(:),INTENT(IN)::array INTEGER(KIND=i_wp),DIMENSION(:),INTENT(INOUT)::permutation LOGICAL,INTENT(IN),OPTIONAL::partially_ranked REAL(KIND=r_dp),INTENT(IN),OPTIONAL::mean_value,standard_deviation CHARACTER(LEN=*),INTENT(IN),OPTIONAL::pivot_selection INTEGER(KIND=i_wp),INTENT(IN),OPTIONAL::cutoff_size OPTIONAL::pivot_function INTERFACE FUNCTION pivot_function(lower_bound,upper_bound)RESULT(pivot_value) USE Precision REAL(KIND=r_sp),INTENT(IN)::lower_bound,upper_bound REAL(KIND=r_dp)::pivot_value END FUNCTION pivot_function END INTERFACE INTEGER(KIND=i_wp)::cutoff INTEGER::recursion_depth,max_recursion_depth CHARACTER(LEN=1)::pivot_method REAL(KIND=r_hv)::mean,std,total REAL(KIND=r_dp)::min_val,max_val LOGICAL::uniform_distribution LOGICAL::comparison_test,use_permutation INTEGER(KIND=i_wp)::index,element,temp_index,temp_element INTEGER::alloc_status REAL(KIND=r_dp)::temp_compare IF(PRESENT(partially_ranked))THEN use_permutation=partially_ranked ELSE use_permutation=.FALSE. END IF IF(.NOT.use_permutation)THEN DO index=1,INT(SIZE(array),KIND=i_wp) permutation(index)=index END DO END IF pivot_method='R' IF(PRESENT(pivot_selection))pivot_method=pivot_selection IF(PRESENT(cutoff_size))THEN cutoff=MAX(1_i_wp,cutoff_size) ELSE cutoff=25 END IF recursion_depth=0; max_recursion_depth=0 SELECT CASE(pivot_method) CASE('U','u','G','g','N','n') IF(.NOT.PRESENT(pivot_function))THEN IF(PRESENT(mean_value))THEN mean=REAL(mean_value,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) total=total+REAL(array(element),KIND=r_hv) END DO mean=total/REAL(INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp),KIND=r_hv) END IF IF(PRESENT(standard_deviation))THEN std=REAL(standard_deviation,KIND=r_hv) ELSE total=0.0_r_hv DO element=1,INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp) total=total+(REAL(array(element),KIND=r_hv)-mean)**2 END DO std=SQRT(total/REAL(INT(INT(SIZE(array),KIND=i_wp),KIND=i_wp),KIND=r_hv)) END IF END IF IF(pivot_method=='U'.OR.pivot_method=='u')THEN CALL QuickSortRecursive_RD_r_dp(1_i_wp,INT(SIZE(array),KIND=i_wp),& range=(/mean-std*SQRT(3.0_r_hv),mean+std*SQRT(3.0_r_hv)/),distribution='U') ELSE CALL QuickSortRecursive_RD_r_dp(1_i_wp,INT(SIZE(array),KIND=i_wp),& range=(/mean-std*10.0_r_hv,mean+std*10.0/),distribution='N') END IF CASE DEFAULT CALL QuickSortRecursive_PE_r_dp(1_i_wp,INT(SIZE(array),KIND=i_wp),(INT(SIZE(array),KIND=i_wp)+1)/2) ENDSELECT IF(debug)WRITE(message_print_unit,*)"Recursion depth=",max_recursion_depth IF(.NOT.PRESENT(cutoff_size))THEN CALL ShellInsertionRank(array=array,permutation=permutation,& partially_ranked=.TRUE.,method="Insertion") END IF CONTAINS RECURSIVE SUBROUTINE QuickSortRecursive_RD_r_dp(starting,ending,range,distribution) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::starting,ending REAL(KIND=r_hv),DIMENSION(2),INTENT(IN)::range CHARACTER(LEN=1),INTENT(IN)::distribution INTEGER(KIND=i_wp)::low,high REAL(KIND=r_dp)::pivot_key recursion_depth=recursion_depth+1 IF((ending-starting)=range(2))THEN pivot_key=(range(1)+range(2))/2 END IF END IF ELSE pivot_key=pivot_function(range(1),range(2)) END IF END IF NibbleArray_RD:DO NibbleLeft_RD:DO IF(low>high)EXIT NibbleArray_RD comparison_test=(array(permutation(low))<=pivot_key) IF(.NOT.comparison_test)EXIT NibbleLeft_RD low=low+1 END DO NibbleLeft_RD NibbleRight_RD:DO IF(highending).OR.(high=MAX(25_i_wp,cutoff))THEN IF(debug)WRITE(message_print_unit,*)"Recalculating range for section of size: ",ending-starting+1 min_val=HUGE(array); max_val=-HUGE(array) DO index=starting,ending temp_compare=array(permutation(index)) IF(temp_comparemax_val)max_val=temp_compare END DO CALL QuickSortRecursive_RD_r_dp(starting,ending,& range=REAL((/min_val,max_val/),KIND=r_hv),distribution=distribution) ELSE CALL QuickSortRecursive_PE_r_dp(starting,ending,(starting+ending)/2) END IF ELSE IF(low-starting>=cutoff)CALL QuickSortRecursive_RD_r_dp(starting,low-1,& range=(/range(1),REAL(pivot_key,KIND=r_hv)/),distribution=distribution) IF(ending-high>=cutoff)CALL QuickSortRecursive_RD_r_dp(high+1,ending,& range=(/REAL(pivot_key,KIND=r_hv),range(2)/),distribution=distribution) END IF max_recursion_depth=MAX(max_recursion_depth,recursion_depth) recursion_depth=recursion_depth-1 RETURN END SUBROUTINE RECURSIVE SUBROUTINE QuickSortRecursive_PE_r_dp(starting,ending,pivot_index) IMPLICIT NONE INTEGER(KIND=i_wp),INTENT(IN)::starting,ending INTEGER(KIND=i_wp),INTENT(IN)::pivot_index INTEGER(KIND=i_wp)::new_pivot_index INTEGER(KIND=i_wp)::low,high REAL(KIND=r_dp)::pivot_key recursion_depth=recursion_depth+1 IF((ending-starting)high)EXIT NibbleArray_PE comparison_test=(array(permutation(low))<=pivot_key) IF(.NOT.comparison_test)EXIT NibbleLeft_PE low=low+1 END DO NibbleLeft_PE NibbleRight_PE:DO IF(high=cutoff)THEN CALL RandomUniform(new_pivot_index,range=(/starting,low-1/)) CALL QuickSortRecursive_PE_r_dp(starting,low-1,new_pivot_index) END IF IF(ending-high>=cutoff)THEN CALL RandomUniform(new_pivot_index,range=(/high+1,ending/)) CALL QuickSortRecursive_PE_r_dp(high+1,ending,new_pivot_index) END IF max_recursion_depth=MAX(max_recursion_depth,recursion_depth) recursion_depth=recursion_depth-1 RETURN END SUBROUTINE END SUBROUTINE Subroutine MergeRank(XVALT,IRNGT) Real,Dimension(:),Intent(In)::XVALT Integer,Dimension(:),Intent(Out)::IRNGT Integer,Dimension(SIZE(IRNGT))::JWRKT Integer::LMTNA,LMTNC,IRNG1,IRNG2 Integer::NVAL,IIND,IWRKD,IWRK,IWRKF,JINDA,IINDA,IINDB Real(Kind(XVALT))::XVALA,XVALB NVAL=Min(SIZE(XVALT),SIZE(IRNGT)) Select Case(NVAL) Case(:0) Return Case(1) IRNGT(1)=1 Return Case Default Continue End Select Do IIND=2,NVAL,2 If(XVALT(IIND-1)<=XVALT(IIND))Then IRNGT(IIND-1)=IIND-1 IRNGT(IIND)=IIND Else IRNGT(IIND-1)=IIND IRNGT(IIND)=IIND-1 End If End Do If(Mod(NVAL,2)/=0)Then IRNGT(NVAL)=NVAL End If LMTNA=2 LMTNC=4 Do If(NVAL<=2)Exit Do IWRKD=0,NVAL-1,4 If((IWRKD+4)>NVAL)Then If((IWRKD+2)>=NVAL)Exit If(XVALT(IRNGT(IWRKD+2))<=XVALT(IRNGT(IWRKD+3)))Exit If(XVALT(IRNGT(IWRKD+1))<=XVALT(IRNGT(IWRKD+3)))Then IRNG2=IRNGT(IWRKD+2) IRNGT(IWRKD+2)=IRNGT(IWRKD+3) IRNGT(IWRKD+3)=IRNG2 Else IRNG1=IRNGT(IWRKD+1) IRNGT(IWRKD+1)=IRNGT(IWRKD+3) IRNGT(IWRKD+3)=IRNGT(IWRKD+2) IRNGT(IWRKD+2)=IRNG1 End If Exit End If If(XVALT(IRNGT(IWRKD+2))<=XVALT(IRNGT(IWRKD+3)))Cycle If(XVALT(IRNGT(IWRKD+1))<=XVALT(IRNGT(IWRKD+3)))Then IRNG2=IRNGT(IWRKD+2) IRNGT(IWRKD+2)=IRNGT(IWRKD+3) If(XVALT(IRNG2)<=XVALT(IRNGT(IWRKD+4)))Then IRNGT(IWRKD+3)=IRNG2 Else IRNGT(IWRKD+3)=IRNGT(IWRKD+4) IRNGT(IWRKD+4)=IRNG2 End If Else IRNG1=IRNGT(IWRKD+1) IRNG2=IRNGT(IWRKD+2) IRNGT(IWRKD+1)=IRNGT(IWRKD+3) If(XVALT(IRNG1)<=XVALT(IRNGT(IWRKD+4)))Then IRNGT(IWRKD+2)=IRNG1 If(XVALT(IRNG2)<=XVALT(IRNGT(IWRKD+4)))Then IRNGT(IWRKD+3)=IRNG2 Else IRNGT(IWRKD+3)=IRNGT(IWRKD+4) IRNGT(IWRKD+4)=IRNG2 End If Else IRNGT(IWRKD+2)=IRNGT(IWRKD+4) IRNGT(IWRKD+3)=IRNG1 IRNGT(IWRKD+4)=IRNG2 End If End If End Do LMTNA=4 Exit End Do Do If(LMTNA>=NVAL)Exit IWRKF=0 LMTNC=2*LMTNC Do IWRK=IWRKF IWRKD=IWRKF+1 JINDA=IWRKF+LMTNA IWRKF=IWRKF+LMTNC If(IWRKF>=NVAL)Then If(JINDA>=NVAL)Exit IWRKF=NVAL End If IINDA=1 IINDB=JINDA+1 JWRKT(1:LMTNA)=IRNGT(IWRKD:JINDA) XVALA=XVALT(JWRKT(IINDA)) XVALB=XVALT(IRNGT(IINDB)) Do IWRK=IWRK+1 If(XVALA>XVALB)Then IRNGT(IWRK)=IRNGT(IINDB) IINDB=IINDB+1 If(IINDB>IWRKF)Then IRNGT(IWRK+1:IWRKF)=JWRKT(IINDA:LMTNA) Exit End If XVALB=XVALT(IRNGT(IINDB)) Else IRNGT(IWRK)=JWRKT(IINDA) IINDA=IINDA+1 If(IINDA>LMTNA)Exit XVALA=XVALT(JWRKT(IINDA)) End If End Do End Do LMTNA=2*LMTNA End Do Return End Subroutine MergeRank END MODULE Sorting_Ranking