Line data Source code
1 :
2 : ! http://fortranwiki.org/fortran/show/qsort_inline
3 :
4 :
5 : !======================================================================
6 : ! Fast in-line QSORT+INSERTION SORT for Fortran.
7 : ! Author: Joseph M. Krahn
8 : ! FILE: qsort_inline.inc
9 : ! PURPOSE:
10 : ! Generate a custom array sort procedure for a specific type,
11 : ! without the comparison-callback overhead of a generic sort procedure.
12 : ! This is essentially the same as an in-line optimization, which generally
13 : ! is not feasible for a library-based generic sort procedure.
14 : !
15 : ! This implementation is as generic as possible, while avoiding the need
16 : ! for a code pre-processor. The success of this approach assumes that
17 : ! internal procedures are always in-lined with optimized Fortran compilation.
18 : !
19 : ! USAGE:
20 : ! This file contains the sort subroutine body. You must supply
21 : ! an integer parameter QSORT_THRESHOLD, and internal procedures:
22 : ! subroutine INIT()
23 : ! logical function LESS_THAN(a,b)
24 : ! subroutine SWAP(a,b)
25 : ! subroutine RSHIFT(left,right)
26 : !
27 : ! Descriptions:
28 : !
29 : ! SUBROUTINE INIT()
30 : ! Any user initialization code. This is needed because executable
31 : ! statements cannot precede this code, which begins with declarations.
32 : ! In many cases, this is just an empty procedure.
33 : !
34 : ! LOGICAL FUNCTION LESS_THAN(a,b)
35 : ! Return TRUE if array member 'a' is less than array member 'b'.
36 : ! Only a TRUE value causes a change in sort order. This minimizes data
37 : ! manipulation, and maintains the original order for values that are
38 : ! equivalent by the sort comparison. It also avoids the need to
39 : ! distinguish equality from greater-than.
40 : !
41 : ! SUBROUTINE SWAP(A,B)
42 : ! Swap array members 'a' and 'b'
43 : !
44 : ! SUBROUTINE RSHIFT(LEFT,RIGHT)
45 : ! Perform a circular shift of array members LEFT through RIGHT,
46 : ! shifting the element at RIGHT back to the position at LEFT.
47 : !
48 : ! QSORT_THRESHOLD:
49 : ! The QSORT is used down to the QSORT_THRESHOLD size sorted blocks.
50 : ! Then insertion sort is used for the remainder, because it is faster
51 : ! for small sort ranges. The optimal size is not critical. Most of
52 : ! the benefit is in blocks of 8 or less, and values of 16 to 128
53 : ! are generally about equal speed. However, the optimal value
54 : ! depends a lot on the hardware and the data being sorted, so this
55 : ! is left as a tunable parameter for cases where there is an
56 : ! effect on performance.
57 : !
58 : !---------------------------------------------------------------------
59 : ! NOTES:
60 : ! The procedure uses a optimized combination of QSORT and INSERTION
61 : ! sorting. The algorithm is based on code used in GLIBC.
62 : ! A stack is used in place of recursive calls. The stack size must
63 : ! be at least as big as the number of bits in the largest array index.
64 : !
65 : ! Sorting vectors of a multidimensional allocatable array can be
66 : ! VERY slow. In this case, or with large derived types, it is better
67 : ! to sort a simple derived type of key/index pairs, then reorder
68 : ! the actual data using the sorted indices.
69 : !
70 : !---------------------------------------------------------------------
71 : integer :: stack_top, right_size, left_size
72 : integer :: mid, left, right, low, high
73 :
74 : ! A stack of 32 can handle the entire extent of a 32-bit
75 : ! index, so this value is fixed. If you have 64-bit indexed
76 : ! arrays, which might contain more than 2^32 elements, this
77 : ! should be set to 64.
78 : integer, parameter :: QSORT_STACK_SIZE = 32, QSORT_THRESHOLD = 10
79 : type qsort_stack; integer :: low, high; end type qsort_stack
80 : type(qsort_stack) :: stack(QSORT_STACK_SIZE)
81 :
82 22 : call init()
83 :
84 22 : if (array_size > QSORT_THRESHOLD) then
85 22 : low = 1
86 22 : high = array_size
87 22 : stack_top = 0
88 :
89 : QSORT_LOOP: &
90 : do
91 173 : mid = (low + high)/2
92 173 : if (LESS_THAN (mid, low)) then
93 63 : call SWAP(mid,low)
94 : end if
95 173 : if (LESS_THAN (high, mid)) then
96 128 : call SWAP(high,mid)
97 128 : if (LESS_THAN (mid, low)) then
98 42 : call SWAP(mid,low)
99 : end if
100 : end if
101 173 : left = low + 1
102 173 : right = high - 1
103 :
104 : COLLAPSE_WALLS: &
105 953 : do
106 2306 : do while (LESS_THAN (left, mid))
107 1180 : left=left+1
108 : end do
109 2956 : do while (LESS_THAN (mid, right))
110 1830 : right=right-1
111 : end do
112 1126 : if (left < right) then
113 953 : call SWAP(left,right)
114 953 : if (mid == left) then
115 87 : mid = right
116 866 : else if (mid == right) then
117 76 : mid = left
118 : end if
119 953 : left=left+1
120 953 : right=right-1
121 : else
122 173 : if (left == right) then
123 10 : left=left+1
124 10 : right=right-1
125 : end if
126 : exit COLLAPSE_WALLS
127 : end if
128 : end do COLLAPSE_WALLS
129 :
130 : ! Set up indices for the next iteration.
131 : ! Determine left and right partition sizes.
132 : ! Defer partitions smaller than the QSORT_THRESHOLD.
133 : ! If both partitions are significant,
134 : ! push the larger one onto the stack.
135 173 : right_size = right - low
136 173 : left_size = high - left
137 195 : if (right_size <= QSORT_THRESHOLD) then
138 104 : if (left_size <= QSORT_THRESHOLD) then
139 : ! Ignore both small partitions: Pop a partition or exit.
140 66 : if (stack_top<1) exit QSORT_LOOP
141 44 : low=stack(stack_top)%low; high=stack(stack_top)%high
142 44 : stack_top=stack_top-1
143 : else
144 : ! Ignore small left partition.
145 38 : low = left
146 : end if
147 69 : else if (left_size <= QSORT_THRESHOLD) then
148 : ! Ignore small right partition.
149 25 : high = right
150 44 : else if (right_size > left_size) then
151 : ! Push larger left partition indices.
152 19 : stack_top=stack_top+1
153 19 : stack(stack_top)=qsort_stack(low,right)
154 19 : low = left
155 : else
156 : ! Push larger right partition indices.
157 25 : stack_top=stack_top+1
158 25 : stack(stack_top)=qsort_stack(left,high)
159 25 : high = right
160 : end if
161 : end do QSORT_LOOP
162 : end if ! (array_size > QSORT_THRESHOLD)
163 :
164 : ! Sort the remaining small partitions using insertion sort,
165 : ! which should be faster for partitions smaller than the
166 : ! appropriate QSORT_THRESHOLD.
167 :
168 : ! First, find smallest element in first QSORT_THRESHOLD and
169 : ! place it at the array's beginning. This places a lower
170 : ! bound 'guard' position, and speeds up the inner loop
171 : ! below, because it will not need a lower-bound test.
172 22 : low = 1
173 22 : high = array_size
174 :
175 : ! left is the MIN_LOC index here:
176 22 : left=low
177 242 : do right = low+1, min(low+QSORT_THRESHOLD,high)
178 242 : if (LESS_THAN(right,left)) left=right
179 : end do
180 22 : if (left/=low) call SWAP(left,low)
181 :
182 : ! Insertion sort, from left to right.
183 : ! (assuming that the left is the lowest numbered index)
184 : INSERTION_SORT: &
185 1431 : do right = low+2,high
186 1409 : left=right-1
187 1431 : if (LESS_THAN(right,left)) then
188 2328 : do while (LESS_THAN(right,left-1))
189 1532 : left=left-1
190 : end do
191 796 : call RSHIFT(left,right)
192 : end if
193 : end do INSERTION_SORT
194 : !--------------------------------------------------------------
|
