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