1! 2! SLATEC: public domain 3! 4*deck isort 5 subroutine isortid (ix, dy, n, kflag) 6c 7c changed on 01.02.2001: auxiliary array is now real*8 8c 9C***BEGIN PROLOGUE ISORT 10C***PURPOSE Sort an array and optionally make the same interchanges in 11C an auxiliary array. The array may be sorted in increasing 12C or decreasing order. A slightly modified QUICKSORT 13C algorithm is used. 14C***LIBRARY SLATEC 15C***CATEGORY N6A2A 16C***TYPE INTEGER (SSORT-S, DSORT-D, ISORT-I) 17C***KEYWORDS SINGLETON QUICKSORT, SORT, SORTING 18C***AUTHOR Jones, R. E., (SNLA) 19C Kahaner, D. K., (NBS) 20C Wisniewski, J. A., (SNLA) 21C***DESCRIPTION 22C 23C ISORT sorts array IX and optionally makes the same interchanges in 24C array DY. The array IX may be sorted in increasing order or 25C decreasing order. A slightly modified quicksort algorithm is used. 26C 27C Description of Parameters 28C IX - integer array of values to be sorted 29C DY - real*8 array to be (optionally) carried along 30C N - number of values in integer array IX to be sorted 31C KFLAG - control parameter 32C = 2 means sort IX in increasing order and carry DY along. 33C = 1 means sort IX in increasing order (ignoring DY) 34C = -1 means sort IX in decreasing order (ignoring DY) 35C = -2 means sort IX in decreasing order and carry DY along. 36C 37C***REFERENCES R. C. Singleton, Algorithm 347, An efficient algorithm 38C for sorting with minimal storage, Communications of 39C the ACM, 12, 3 (1969), pp. 185-187. 40C***ROUTINES CALLED XERMSG 41C***REVISION HISTORY (YYMMDD) 42C 761118 DATE WRITTEN 43C 810801 Modified by David K. Kahaner. 44C 890531 Changed all specific intrinsics to generic. (WRB) 45C 890831 Modified array declarations. (WRB) 46C 891009 Removed unreferenced statement labels. (WRB) 47C 891009 REVISION DATE from Version 3.2 48C 891214 Prologue converted to Version 4.0 format. (BAB) 49C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 50C 901012 Declared all variables; changed X,Y to IX,DY. (M. McClain) 51C 920501 Reformatted the REFERENCES section. (DWL, WRB) 52C 920519 Clarified error messages. (DWL) 53C 920801 Declarations section rebuilt and code restructured to use 54C IF-THEN-ELSE-ENDIF. (RWC, WRB) 55! 100411 changed the dimension of IL and IU from 21 to 31. 56! 57! field IL and IU have the dimension 31. This is log2 of the largest 58! array size to be sorted. If arrays larger than 2**31 in length have 59! to be sorted, this dimension has to be modified accordingly 60! 61! 62 implicit none 63C***END PROLOGUE ISORT 64C .. Scalar Arguments .. 65 integer kflag, n 66C .. Array Arguments .. 67 integer ix(*) 68 real*8 dy(*),ty,tty 69C .. Local Scalars .. 70 real r 71 integer i, ij, j, k, kk, l, m, nn, t, tt 72C .. Local Arrays .. 73 integer il(31), iu(31) 74C .. External Subroutines .. 75! EXTERNAL XERMSG 76C .. Intrinsic Functions .. 77 intrinsic abs, int 78C***FIRST EXECUTABLE STATEMENT ISORT 79 nn = n 80 if (nn .lt. 1) then 81! CALL XERMSG ('SLATEC', 'ISORT', 82! + 'The number of values to be sorted is not positive.', 1, 1) 83 return 84 endif 85C 86 kk = abs(kflag) 87 if (kk.ne.1 .and. kk.ne.2) then 88! CALL XERMSG ('SLATEC', 'ISORT', 89! + 'The sort control parameter, K, is not 2, 1, -1, or -2.', 2, 90! + 1) 91 return 92 endif 93C 94C Alter array IX to get decreasing order if needed 95C 96 if (kflag .le. -1) then 97 do 10 i=1,nn 98 ix(i) = -ix(i) 99 10 continue 100 endif 101C 102 if (kk .eq. 2) go to 100 103C 104C Sort IX only 105C 106 m = 1 107 i = 1 108 j = nn 109 r = 0.375e0 110C 111 20 if (i .eq. j) go to 60 112 if (r .le. 0.5898437e0) then 113 r = r+3.90625e-2 114 else 115 r = r-0.21875e0 116 endif 117C 118 30 k = i 119C 120C Select a central element of the array and save it in location T 121C 122 ij = i + int((j-i)*r) 123 t = ix(ij) 124C 125C If first element of array is greater than T, interchange with T 126C 127 if (ix(i) .gt. t) then 128 ix(ij) = ix(i) 129 ix(i) = t 130 t = ix(ij) 131 endif 132 l = j 133C 134C If last element of array is less than than T, interchange with T 135C 136 if (ix(j) .lt. t) then 137 ix(ij) = ix(j) 138 ix(j) = t 139 t = ix(ij) 140C 141C If first element of array is greater than T, interchange with T 142C 143 if (ix(i) .gt. t) then 144 ix(ij) = ix(i) 145 ix(i) = t 146 t = ix(ij) 147 endif 148 endif 149C 150C Find an element in the second half of the array which is smaller 151C than T 152C 153 40 l = l-1 154 if (ix(l) .gt. t) go to 40 155C 156C Find an element in the first half of the array which is greater 157C than T 158C 159 50 k = k+1 160 if (ix(k) .lt. t) go to 50 161C 162C Interchange these elements 163C 164 if (k .le. l) then 165 tt = ix(l) 166 ix(l) = ix(k) 167 ix(k) = tt 168 go to 40 169 endif 170C 171C Save upper and lower subscripts of the array yet to be sorted 172C 173 if (l-i .gt. j-k) then 174 il(m) = i 175 iu(m) = l 176 i = k 177 m = m+1 178 else 179 il(m) = k 180 iu(m) = j 181 j = l 182 m = m+1 183 endif 184 go to 70 185C 186C Begin again on another portion of the unsorted array 187C 188 60 m = m-1 189 if (m .eq. 0) go to 190 190 i = il(m) 191 j = iu(m) 192C 193 70 if (j-i .ge. 1) go to 30 194 if (i .eq. 1) go to 20 195 i = i-1 196C 197 80 i = i+1 198 if (i .eq. j) go to 60 199 t = ix(i+1) 200 if (ix(i) .le. t) go to 80 201 k = i 202C 203 90 ix(k+1) = ix(k) 204 k = k-1 205 if (t .lt. ix(k)) go to 90 206 ix(k+1) = t 207 go to 80 208C 209C Sort IX and carry DY along 210C 211 100 m = 1 212 i = 1 213 j = nn 214 r = 0.375e0 215C 216 110 if (i .eq. j) go to 150 217 if (r .le. 0.5898437e0) then 218 r = r+3.90625e-2 219 else 220 r = r-0.21875e0 221 endif 222C 223 120 k = i 224C 225C Select a central element of the array and save it in location T 226C 227 ij = i + int((j-i)*r) 228 t = ix(ij) 229 ty = dy(ij) 230C 231C If first element of array is greater than T, interchange with T 232C 233 if (ix(i) .gt. t) then 234 ix(ij) = ix(i) 235 ix(i) = t 236 t = ix(ij) 237 dy(ij) = dy(i) 238 dy(i) = ty 239 ty = dy(ij) 240 endif 241 l = j 242C 243C If last element of array is less than T, interchange with T 244C 245 if (ix(j) .lt. t) then 246 ix(ij) = ix(j) 247 ix(j) = t 248 t = ix(ij) 249 dy(ij) = dy(j) 250 dy(j) = ty 251 ty = dy(ij) 252C 253C If first element of array is greater than T, interchange with T 254C 255 if (ix(i) .gt. t) then 256 ix(ij) = ix(i) 257 ix(i) = t 258 t = ix(ij) 259 dy(ij) = dy(i) 260 dy(i) = ty 261 ty = dy(ij) 262 endif 263 endif 264C 265C Find an element in the second half of the array which is smaller 266C than T 267C 268 130 l = l-1 269 if (ix(l) .gt. t) go to 130 270C 271C Find an element in the first half of the array which is greater 272C than T 273C 274 140 k = k+1 275 if (ix(k) .lt. t) go to 140 276C 277C Interchange these elements 278C 279 if (k .le. l) then 280 tt = ix(l) 281 ix(l) = ix(k) 282 ix(k) = tt 283 tty = dy(l) 284 dy(l) = dy(k) 285 dy(k) = tty 286 go to 130 287 endif 288C 289C Save upper and lower subscripts of the array yet to be sorted 290C 291 if (l-i .gt. j-k) then 292 il(m) = i 293 iu(m) = l 294 i = k 295 m = m+1 296 else 297 il(m) = k 298 iu(m) = j 299 j = l 300 m = m+1 301 endif 302 go to 160 303C 304C Begin again on another portion of the unsorted array 305C 306 150 m = m-1 307 if (m .eq. 0) go to 190 308 i = il(m) 309 j = iu(m) 310C 311 160 if (j-i .ge. 1) go to 120 312 if (i .eq. 1) go to 110 313 i = i-1 314C 315 170 i = i+1 316 if (i .eq. j) go to 150 317 t = ix(i+1) 318 ty = dy(i+1) 319 if (ix(i) .le. t) go to 170 320 k = i 321C 322 180 ix(k+1) = ix(k) 323 dy(k+1) = dy(k) 324 k = k-1 325 if (t .lt. ix(k)) go to 180 326 ix(k+1) = t 327 dy(k+1) = ty 328 go to 170 329C 330C Clean up 331C 332 190 if (kflag .le. -1) then 333 do 200 i=1,nn 334 ix(i) = -ix(i) 335 200 continue 336 endif 337 return 338 end 339