Multiple Key Sort in Fortran
in [Fortran]
|
From: Joe Krahn on 3 Mar 2006 15:37 Archie wrote: > I am looking for Fortran 77/95 source code for a multiple key sort. > > Input records can contain a variable number of (integer) fields where > the number of fields is typically between 4 and 40. For a file of 40 > field records, the records need to be sorted into ascending order of > field 40 within field 39, withing field 38 .... within field 1. > > File sizes are usually relatively small with between, say, 1,000 and > 20,000 records. > > It would generally be better to sort within memory, for subsequent > processing, rather than writing to an output file. > > The sort might optionally be repeated based on subsets of the input > data file. > > Most sort code I have been able to find is for a single sort key. > > Is anyone able to recommend suitable Fortran source code? > I have put together an customizable derived-type sort function, based on GLIBC's QSORT, which is actually QSORT down to a threshold, where it switches to INSERTION sorting, which is faster for smaller ranges. The ideal threshold depends on the data size and the CPU. The compare function is user-defined, typically an internal function, which has a good chance of being inlined. You can do any comparisons you want there. In my tests, it is almost as fast as an inline QSORT in C, but both are considerably faster than standard library functions because there's no callback overhead. It's small enough to post in this message, so here it is. There are some long comment lines which may get wrapped. Let me know Joe Krahn !/* ! Fast inline QSORT for Fortran. -- Joe Krahn, 2006 ! Filename: qsort_inline.inc ! ! Sorts a one-dimensional array; normally of a derived-type ! containing a sort key. ! ! Uses a optimized combination of QSORT and INSERTION sorting. ! The algorithm is based on code used in GLIBC. ! ! Using a multidimensional allocatable array to hold ! multiple data columns is VERY slow. ! It is better to sort a derived type of key/index pairs, ! then reorder data using the ! sorted indices. ! ! You must define the function QSORT_COMPARE(), ! and the integer QSORT_THRESHOLD ! ! logical function QSORT_COMPARE(a,b) ! Must be an internal function -- it accesses array() from the parent. ! Should return the logical result of (array(a)%key < array(b)%key) ! assuming an ascending order sort. ! ! NOTE: For speed, the compare function should be inlined. ! If in doubt, use a CPP macro instead of a proper Fortran function. ! ! QSORT_THRESHOLD: ! The QSORT is used down to the QSORT_THRESHOLD size sorted blocks. Then ! insertion sort is used for the remainder, because it is faster for small sort ranges. ! ! Large sorting elements or small cache may make a smaller threshold more useful. ! You can also set this to a run-time argument/variable with no performance loss. ! ! EXAMPLE: ! ! subroutine QSORT_DATA(array) ! implicit none ! type(DATA_TYPE) :: hold, array(:) ! integer, parameter :: QSORT_THRESHOLD = 96 ! include "qsort_inline.inc" ! contains ! logical function QSORT_COMPARE(a,b) ! QSORT_COMPARE = ( array(a)%key < array(b)%key ) ! end function QSORT_COMPARE ! end subroutine QSORT_DATA !*/ integer :: stack_top, right_size,left_size integer :: mid, left, right, low, high ! A stack of 32 can handle the entire extent of a 32-bit index. ! Use 64 for 64-bit indexed arrays, if they might contain more ! than 2^32 elements. integer, parameter :: QSORT_STACK_SIZE = 32 type qsort_stack; integer :: low, high; end type type(qsort_stack) :: stack(QSORT_STACK_SIZE) if (size(array,1) > QSORT_THRESHOLD) then low = lbound(array,1) high = ubound(array,1) stack_top = 0 QSORT_LOOP: do mid = (low + high)/2 if (QSORT_COMPARE (mid, low)) then hold=array(mid);array(mid)=array(low);array(low)=hold; end if if (QSORT_COMPARE (high, mid)) then hold=array(high);array(high)=array(mid);array(mid)=hold; if (QSORT_COMPARE (mid, low)) then hold=array(mid);array(mid)=array(low);array(low)=hold; end if end if left = low + 1 right = high - 1 COLLAPSE_WALLS: do do while (QSORT_COMPARE (left, mid)); left=left+1; end do do while (QSORT_COMPARE (mid, right)); right=right-1; end do if (left < right) then hold=array(left);array(left)=array(right);array(right)=hold; if (mid == left) then mid = right else if (mid == right) then mid = left end if left=left+1 right=right-1 else if (left == right) then left=left+1 right=right-1 end if exit COLLAPSE_WALLS end if end do COLLAPSE_WALLS ! Set up indices for the next iteration. ! Determine left and right partition sizes. ! Defer partitions smaller than the QSORT_THRESHOLD. ! If both partitions are significant, push the larger one onto the stack. right_size = right - low left_size = high - left if (right_size <= QSORT_THRESHOLD) then if (left_size <= QSORT_THRESHOLD) then ! Ignore both small partitions: Pop a partition or exit. if (stack_top<1) exit QSORT_LOOP low=stack(stack_top)%low; high=stack(stack_top)%high stack_top=stack_top-1 else ! Ignore small left partition. low = left end if else if (left_size <= QSORT_THRESHOLD) then ! Ignore small right partition. high = right else if (right_size > left_size) then ! Push larger left partition indices. stack_top=stack_top+1 stack(stack_top)=qsort_stack(low,right) low = left else ! Push larger right partition indices. stack_top=stack_top+1 stack(stack_top)=qsort_stack(left,high) high = right end if end do QSORT_LOOP end if ! Sort the remaining small partitions using insertion sort, which ! should be faster for partitions smaller than the appropriate ! QSORT_THRESHOLD. ! First, find smallest element in first QSORT_THRESHOLD and place ! it at the array's beginning. This is the smallest array element, ! and the operation speeds up insertion sort's inner loop. low = lbound(array,1) high = ubound(array,1) left=low do right = low+1, MIN(low+QSORT_THRESHOLD,high) if (QSORT_COMPARE(right,left)) left=right end do if (left/=low) then hold=array(left);array(left)=array(low);array(low)=hold; end if ! Insertion sort, starting from left-hand-side, ! moving to right-hand-side. INSERTION_SORT: do right = low+2,high left=right-1 if (QSORT_COMPARE(right,left)) then do while (QSORT_COMPARE (right, left-1)); left=left-1; end do hold=array(right) array(left+1:right)=array(left:right-1) array(left)=hold end if end do INSERTION_SORT
From: Brooks Moses on 3 Mar 2006 17:11 Joe Krahn wrote: > I have put together an customizable derived-type sort function, based on > GLIBC's QSORT, which is actually QSORT down to a threshold, where it > switches to INSERTION sorting, which is faster for smaller ranges. The > ideal threshold depends on the data size and the CPU. > > The compare function is user-defined, typically an internal function, > which has a good chance of being inlined. You can do any comparisons you > want there. > > In my tests, it is almost as fast as an inline QSORT in C, but both are > considerably faster than standard library functions because there's no > callback overhead. > > It's small enough to post in this message, so here it is. Thanks; this looks useful! Since you didn't explicitly state this: What your preference with regard to redistributing this? In particular, would it be ok with you if I were to include it (with proper attribution, of course!) as part of a program that's distributed under an open-source license? (I don't have a particular program in mind at the moment, but if I do find it useful in the future, it's likely to be in something like that.) - Brooks -- The "bmoses-nospam" address is valid; no unmunging needed.
From: Joe Krahn on 3 Mar 2006 18:52
Brooks Moses wrote: > Joe Krahn wrote: > >> I have put together an customizable derived-type sort function, based >> on GLIBC's QSORT, which is actually QSORT down to a threshold, where >> it switches to INSERTION sorting, which is faster for smaller ranges. >> The ideal threshold depends on the data size and the CPU. >> >> The compare function is user-defined, typically an internal function, >> which has a good chance of being inlined. You can do any comparisons >> you want there. >> >> In my tests, it is almost as fast as an inline QSORT in C, but both >> are considerably faster than standard library functions because >> there's no callback overhead. >> >> It's small enough to post in this message, so here it is. > > > Thanks; this looks useful! > > Since you didn't explicitly state this: What your preference with regard > to redistributing this? In particular, would it be ok with you if I > were to include it (with proper attribution, of course!) as part of a > program that's distributed under an open-source license? > > (I don't have a particular program in mind at the moment, but if I do > find it useful in the future, it's likely to be in something like that.) > > - Brooks > > I was hoping to get a bit of feedback and then make a 'proper' version available with a bit more info in the header for redistribution purposes. So, let me know if it needs more documentation or features. I'm hoping to put together a number of general-purpose NON-math utility sources together. The lack of good, non-math utilities is one of the biggest (IMHO) deficiencies in Fortran. With better strings and C binding on the horizon, many things will be possible in Fortran that were either difficult or best done outside of Fortran. Joe |