1c----------------------------------------------------------------------- 2c\BeginDoc 3c 4c\Name: dsortc 5c 6c\Description: 7c Sorts the complex array in XREAL and XIMAG into the order 8c specified by WHICH and optionally applies the permutation to the 9c real array Y. It is assumed that if an element of XIMAG is 10c nonzero, then its negative is also an element. In other words, 11c both members of a complex conjugate pair are to be sorted and the 12c pairs are kept adjacent to each other. 13c 14c\Usage: 15c call dsortc 16c ( WHICH, APPLY, N, XREAL, XIMAG, Y ) 17c 18c\Arguments 19c WHICH Character*2. (Input) 20c 'LM' -> sort XREAL,XIMAG into increasing order of magnitude. 21c 'SM' -> sort XREAL,XIMAG into decreasing order of magnitude. 22c 'LR' -> sort XREAL into increasing order of algebraic. 23c 'SR' -> sort XREAL into decreasing order of algebraic. 24c 'LI' -> sort XIMAG into increasing order of magnitude. 25c 'SI' -> sort XIMAG into decreasing order of magnitude. 26c NOTE: If an element of XIMAG is non-zero, then its negative 27c is also an element. 28c 29c APPLY Logical. (Input) 30c APPLY = .TRUE. -> apply the sorted order to array Y. 31c APPLY = .FALSE. -> do not apply the sorted order to array Y. 32c 33c N Integer. (INPUT) 34c Size of the arrays. 35c 36c XREAL, Double precision array of length N. (INPUT/OUTPUT) 37c XIMAG Real and imaginary part of the array to be sorted. 38c 39c Y Double precision array of length N. (INPUT/OUTPUT) 40c 41c\EndDoc 42c 43c----------------------------------------------------------------------- 44c 45c\BeginLib 46c 47c\Author 48c Danny Sorensen Phuong Vu 49c Richard Lehoucq CRPC / Rice University 50c Dept. of Computational & Houston, Texas 51c Applied Mathematics 52c Rice University 53c Houston, Texas 54c 55c\Revision history: 56c xx/xx/92: Version ' 2.1' 57c Adapted from the sort routine in LANSO. 58c 59c\SCCS Information: @(#) 60c FILE: sortc.F SID: 2.3 DATE OF SID: 4/20/96 RELEASE: 2 61c 62c\EndLib 63c 64c----------------------------------------------------------------------- 65c 66 subroutine dsortc (which, apply, n, xreal, ximag, y) 67c 68c %------------------% 69c | Scalar Arguments | 70c %------------------% 71c 72 character*2 which 73 logical apply 74 integer n 75c 76c %-----------------% 77c | Array Arguments | 78c %-----------------% 79c 80 Double precision 81 & xreal(0:n-1), ximag(0:n-1), y(0:n-1) 82c 83c %---------------% 84c | Local Scalars | 85c %---------------% 86c 87 integer i, igap, j 88 Double precision 89 & temp, temp1, temp2 90c 91c %--------------------% 92c | External Functions | 93c %--------------------% 94c 95 Double precision 96 & dlapy2 97 external dlapy2 98c 99c %-----------------------% 100c | Executable Statements | 101c %-----------------------% 102c 103 igap = n / 2 104c 105 if (which .eq. 'LM') then 106c 107c %------------------------------------------------------% 108c | Sort XREAL,XIMAG into increasing order of magnitude. | 109c %------------------------------------------------------% 110c 111 10 continue 112 if (igap .eq. 0) go to 9000 113c 114 do 30 i = igap, n-1 115 j = i-igap 116 20 continue 117c 118 if (j.lt.0) go to 30 119c 120 temp1 = dlapy2(xreal(j),ximag(j)) 121 temp2 = dlapy2(xreal(j+igap),ximag(j+igap)) 122c 123 if (temp1.gt.temp2) then 124 temp = xreal(j) 125 xreal(j) = xreal(j+igap) 126 xreal(j+igap) = temp 127c 128 temp = ximag(j) 129 ximag(j) = ximag(j+igap) 130 ximag(j+igap) = temp 131c 132 if (apply) then 133 temp = y(j) 134 y(j) = y(j+igap) 135 y(j+igap) = temp 136 end if 137 else 138 go to 30 139 end if 140 j = j-igap 141 go to 20 142 30 continue 143 igap = igap / 2 144 go to 10 145c 146 else if (which .eq. 'SM') then 147c 148c %------------------------------------------------------% 149c | Sort XREAL,XIMAG into decreasing order of magnitude. | 150c %------------------------------------------------------% 151c 152 40 continue 153 if (igap .eq. 0) go to 9000 154c 155 do 60 i = igap, n-1 156 j = i-igap 157 50 continue 158c 159 if (j .lt. 0) go to 60 160c 161 temp1 = dlapy2(xreal(j),ximag(j)) 162 temp2 = dlapy2(xreal(j+igap),ximag(j+igap)) 163c 164 if (temp1.lt.temp2) then 165 temp = xreal(j) 166 xreal(j) = xreal(j+igap) 167 xreal(j+igap) = temp 168c 169 temp = ximag(j) 170 ximag(j) = ximag(j+igap) 171 ximag(j+igap) = temp 172c 173 if (apply) then 174 temp = y(j) 175 y(j) = y(j+igap) 176 y(j+igap) = temp 177 end if 178 else 179 go to 60 180 endif 181 j = j-igap 182 go to 50 183 60 continue 184 igap = igap / 2 185 go to 40 186c 187 else if (which .eq. 'LR') then 188c 189c %------------------------------------------------% 190c | Sort XREAL into increasing order of algebraic. | 191c %------------------------------------------------% 192c 193 70 continue 194 if (igap .eq. 0) go to 9000 195c 196 do 90 i = igap, n-1 197 j = i-igap 198 80 continue 199c 200 if (j.lt.0) go to 90 201c 202 if (xreal(j).gt.xreal(j+igap)) then 203 temp = xreal(j) 204 xreal(j) = xreal(j+igap) 205 xreal(j+igap) = temp 206c 207 temp = ximag(j) 208 ximag(j) = ximag(j+igap) 209 ximag(j+igap) = temp 210c 211 if (apply) then 212 temp = y(j) 213 y(j) = y(j+igap) 214 y(j+igap) = temp 215 end if 216 else 217 go to 90 218 endif 219 j = j-igap 220 go to 80 221 90 continue 222 igap = igap / 2 223 go to 70 224c 225 else if (which .eq. 'SR') then 226c 227c %------------------------------------------------% 228c | Sort XREAL into decreasing order of algebraic. | 229c %------------------------------------------------% 230c 231 100 continue 232 if (igap .eq. 0) go to 9000 233 do 120 i = igap, n-1 234 j = i-igap 235 110 continue 236c 237 if (j.lt.0) go to 120 238c 239 if (xreal(j).lt.xreal(j+igap)) then 240 temp = xreal(j) 241 xreal(j) = xreal(j+igap) 242 xreal(j+igap) = temp 243c 244 temp = ximag(j) 245 ximag(j) = ximag(j+igap) 246 ximag(j+igap) = temp 247c 248 if (apply) then 249 temp = y(j) 250 y(j) = y(j+igap) 251 y(j+igap) = temp 252 end if 253 else 254 go to 120 255 endif 256 j = j-igap 257 go to 110 258 120 continue 259 igap = igap / 2 260 go to 100 261c 262 else if (which .eq. 'LI') then 263c 264c %------------------------------------------------% 265c | Sort XIMAG into increasing order of magnitude. | 266c %------------------------------------------------% 267c 268 130 continue 269 if (igap .eq. 0) go to 9000 270 do 150 i = igap, n-1 271 j = i-igap 272 140 continue 273c 274 if (j.lt.0) go to 150 275c 276 if (abs(ximag(j)).gt.abs(ximag(j+igap))) then 277 temp = xreal(j) 278 xreal(j) = xreal(j+igap) 279 xreal(j+igap) = temp 280c 281 temp = ximag(j) 282 ximag(j) = ximag(j+igap) 283 ximag(j+igap) = temp 284c 285 if (apply) then 286 temp = y(j) 287 y(j) = y(j+igap) 288 y(j+igap) = temp 289 end if 290 else 291 go to 150 292 endif 293 j = j-igap 294 go to 140 295 150 continue 296 igap = igap / 2 297 go to 130 298c 299 else if (which .eq. 'SI') then 300c 301c %------------------------------------------------% 302c | Sort XIMAG into decreasing order of magnitude. | 303c %------------------------------------------------% 304c 305 160 continue 306 if (igap .eq. 0) go to 9000 307 do 180 i = igap, n-1 308 j = i-igap 309 170 continue 310c 311 if (j.lt.0) go to 180 312c 313 if (abs(ximag(j)).lt.abs(ximag(j+igap))) then 314 temp = xreal(j) 315 xreal(j) = xreal(j+igap) 316 xreal(j+igap) = temp 317c 318 temp = ximag(j) 319 ximag(j) = ximag(j+igap) 320 ximag(j+igap) = temp 321c 322 if (apply) then 323 temp = y(j) 324 y(j) = y(j+igap) 325 y(j+igap) = temp 326 end if 327 else 328 go to 180 329 endif 330 j = j-igap 331 go to 170 332 180 continue 333 igap = igap / 2 334 go to 160 335 end if 336c 337 9000 continue 338 return 339c 340c %---------------% 341c | End of dsortc | 342c %---------------% 343c 344 end 345