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