Télécharger dsortc.eso

Retour à la liste

Numérotation des lignes :

dsortc
  1. C DSORTC SOURCE BP208322 15/10/13 21:15:54 8670
  2. c-----------------------------------------------------------------------
  3. c\BeginDoc
  4. c
  5. c\Name: dsortc
  6. c
  7. c\Description:
  8. c Sorts the complex array in XREAL and XIMAG into the order
  9. c specified by WHICH and optionally applies the permutation to the
  10. c real array Y. It is assumed that if an element of XIMAG is
  11. c nonzero, then its negative is also an element. In other words,
  12. c both members of a complex conjugate pair are to be sorted and the
  13. c pairs are kept adjacent to each other.
  14. c
  15. c\Usage:
  16. c call dsortc
  17. c ( WHICH, APPLY, N, XREAL, XIMAG, Y )
  18. c
  19. c\Arguments
  20. c WHICH Character*2. (Input)
  21. c 'LM' -> sort XREAL,XIMAG into increasing order of magnitude.
  22. c 'SM' -> sort XREAL,XIMAG into decreasing order of magnitude.
  23. c 'LR' -> sort XREAL into increasing order of algebraic.
  24. c 'SR' -> sort XREAL into decreasing order of algebraic.
  25. c 'LI' -> sort XIMAG into increasing order of magnitude.
  26. c 'SI' -> sort XIMAG into decreasing order of magnitude.
  27. c NOTE: If an element of XIMAG is non-zero, then its negative
  28. c is also an element.
  29. c
  30. c APPLY Logical. (Input)
  31. c APPLY = .TRUE. -> apply the sorted order to array Y.
  32. c APPLY = .FALSE. -> do not apply the sorted order to array Y.
  33. c
  34. c N Integer. (INPUT)
  35. c Size of the arrays.
  36. c
  37. c XREAL, Double precision array of length N. (INPUT/OUTPUT)
  38. c XIMAG Real and imaginary part of the array to be sorted.
  39. c
  40. c Y REAL*8 array of length N. (INPUT/OUTPUT)
  41. c
  42. c\EndDoc
  43. c
  44. c-----------------------------------------------------------------------
  45. c
  46. c\BeginLib
  47. c
  48. c\Author
  49. c Danny Sorensen Phuong Vu
  50. c Richard Lehoucq CRPC / Rice University
  51. c Dept. of Computational & Houston, Texas
  52. c Applied Mathematics
  53. c Rice University
  54. c Houston, Texas
  55. c
  56. c\Revision history:
  57. c xx/xx/92: Version ' 2.1'
  58. c Adapted from the sort routine in LANSO.
  59. c
  60. c\SCCS Information: @(#)
  61. c FILE: sortc.F SID: 2.3 DATE OF SID: 4/20/96 RELEASE: 2
  62. c
  63. c\EndLib
  64. c
  65. c-----------------------------------------------------------------------
  66. c
  67. subroutine dsortc (which, apply, n, xreal, ximag, y)
  68. c
  69. c %------------------%
  70. c | Scalar Arguments |
  71. c %------------------%
  72. c
  73. character*2 which
  74. logical apply
  75. integer n
  76. c
  77. c %-----------------%
  78. c | Array Arguments |
  79. c %-----------------%
  80. c
  81. REAL*8
  82. & xreal(0:n-1), ximag(0:n-1), y(0:n-1)
  83. c
  84. c %---------------%
  85. c | Local Scalars |
  86. c %---------------%
  87. c
  88. integer i, igap, j
  89. REAL*8
  90. & temp, temp1, temp2
  91. c
  92. c %--------------------%
  93. c | External Functions |
  94. c %--------------------%
  95. c
  96. REAL*8
  97. external dlapy2
  98. c
  99. c %-----------------------%
  100. c | Executable Statements |
  101. c %-----------------------%
  102. c
  103. igap = n / 2
  104. c
  105. if (which .eq. 'LM') then
  106. c
  107. c %------------------------------------------------------%
  108. c | Sort XREAL,XIMAG into increasing order of magnitude. |
  109. c %------------------------------------------------------%
  110. c
  111. 10 continue
  112. if (igap .eq. 0) go to 9000
  113. c
  114. do 30 i = igap, n-1
  115. j = i-igap
  116. 20 continue
  117. c
  118. if (j.lt.0) go to 30
  119. c
  120. temp1 = dlapy2(xreal(j),ximag(j))
  121. temp2 = dlapy2(xreal(j+igap),ximag(j+igap))
  122. c
  123. if (temp1.gt.temp2) then
  124. temp = xreal(j)
  125. xreal(j) = xreal(j+igap)
  126. xreal(j+igap) = temp
  127. c
  128. temp = ximag(j)
  129. ximag(j) = ximag(j+igap)
  130. ximag(j+igap) = temp
  131. c
  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
  145. c
  146. else if (which .eq. 'SM') then
  147. c
  148. c %------------------------------------------------------%
  149. c | Sort XREAL,XIMAG into decreasing order of magnitude. |
  150. c %------------------------------------------------------%
  151. c
  152. 40 continue
  153. if (igap .eq. 0) go to 9000
  154. c
  155. do 60 i = igap, n-1
  156. j = i-igap
  157. 50 continue
  158. c
  159. if (j .lt. 0) go to 60
  160. c
  161. temp1 = dlapy2(xreal(j),ximag(j))
  162. temp2 = dlapy2(xreal(j+igap),ximag(j+igap))
  163. c
  164. if (temp1.lt.temp2) then
  165. temp = xreal(j)
  166. xreal(j) = xreal(j+igap)
  167. xreal(j+igap) = temp
  168. c
  169. temp = ximag(j)
  170. ximag(j) = ximag(j+igap)
  171. ximag(j+igap) = temp
  172. c
  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
  186. c
  187. else if (which .eq. 'LR') then
  188. c
  189. c %------------------------------------------------%
  190. c | Sort XREAL into increasing order of algebraic. |
  191. c %------------------------------------------------%
  192. c
  193. 70 continue
  194. if (igap .eq. 0) go to 9000
  195. c
  196. do 90 i = igap, n-1
  197. j = i-igap
  198. 80 continue
  199. c
  200. if (j.lt.0) go to 90
  201. c
  202. if (xreal(j).gt.xreal(j+igap)) then
  203. temp = xreal(j)
  204. xreal(j) = xreal(j+igap)
  205. xreal(j+igap) = temp
  206. c
  207. temp = ximag(j)
  208. ximag(j) = ximag(j+igap)
  209. ximag(j+igap) = temp
  210. c
  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
  224. c
  225. else if (which .eq. 'SR') then
  226. c
  227. c %------------------------------------------------%
  228. c | Sort XREAL into decreasing order of algebraic. |
  229. c %------------------------------------------------%
  230. c
  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
  236. c
  237. if (j.lt.0) go to 120
  238. c
  239. if (xreal(j).lt.xreal(j+igap)) then
  240. temp = xreal(j)
  241. xreal(j) = xreal(j+igap)
  242. xreal(j+igap) = temp
  243. c
  244. temp = ximag(j)
  245. ximag(j) = ximag(j+igap)
  246. ximag(j+igap) = temp
  247. c
  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
  261. c
  262. else if (which .eq. 'LI') then
  263. c
  264. c %------------------------------------------------%
  265. c | Sort XIMAG into increasing order of magnitude. |
  266. c %------------------------------------------------%
  267. c
  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
  273. c
  274. if (j.lt.0) go to 150
  275. c
  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
  280. c
  281. temp = ximag(j)
  282. ximag(j) = ximag(j+igap)
  283. ximag(j+igap) = temp
  284. c
  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
  298. c
  299. else if (which .eq. 'SI') then
  300. c
  301. c %------------------------------------------------%
  302. c | Sort XIMAG into decreasing order of magnitude. |
  303. c %------------------------------------------------%
  304. c
  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
  310. c
  311. if (j.lt.0) go to 180
  312. c
  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
  317. c
  318. temp = ximag(j)
  319. ximag(j) = ximag(j+igap)
  320. ximag(j+igap) = temp
  321. c
  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
  336. c
  337. 9000 continue
  338. return
  339. c
  340. c %---------------%
  341. c | End of dsortc |
  342. c %---------------%
  343. c
  344. end
  345.  
  346.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales