Télécharger la plaquette Cast3M


fr en

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
  98.       external   dlapy2
  99. c
  100. c     %-----------------------%
  101. c     | Executable Statements |
  102. c     %-----------------------%
  103. c
  104.       igap = n / 2
  105. c
  106.       if (which .eq. 'LM') then
  107. c
  108. c        %------------------------------------------------------%
  109. c        | Sort XREAL,XIMAG into increasing order of magnitude. |
  110. c        %------------------------------------------------------%
  111. c
  112.    10    continue
  113.          if (igap .eq. 0) go to 9000
  114. c
  115.          do 30 i = igap, n-1
  116.             j = i-igap
  117.    20       continue
  118. c
  119.             if (j.lt.0) go to 30
  120. c
  121.             temp1 = dlapy2(xreal(j),ximag(j))
  122.             temp2 = dlapy2(xreal(j+igap),ximag(j+igap))
  123. c
  124.             if (temp1.gt.temp2) then
  125.                 temp = xreal(j)
  126.                 xreal(j) = xreal(j+igap)
  127.                 xreal(j+igap) = temp
  128. c
  129.                 temp = ximag(j)
  130.                 ximag(j) = ximag(j+igap)
  131.                 ximag(j+igap) = temp
  132. c
  133.                 if (apply) then
  134.                     temp = y(j)
  135.                     y(j) = y(j+igap)
  136.                     y(j+igap) = temp
  137.                 end if
  138.             else
  139.                 go to 30
  140.             end if
  141.             j = j-igap
  142.             go to 20
  143.    30    continue
  144.          igap = igap / 2
  145.          go to 10
  146. c
  147.       else if (which .eq. 'SM') then
  148. c
  149. c        %------------------------------------------------------%
  150. c        | Sort XREAL,XIMAG into decreasing order of magnitude. |
  151. c        %------------------------------------------------------%
  152. c
  153.    40    continue
  154.          if (igap .eq. 0) go to 9000
  155. c
  156.          do 60 i = igap, n-1
  157.             j = i-igap
  158.    50       continue
  159. c
  160.             if (j .lt. 0) go to 60
  161. c
  162.             temp1 = dlapy2(xreal(j),ximag(j))
  163.             temp2 = dlapy2(xreal(j+igap),ximag(j+igap))
  164. c
  165.             if (temp1.lt.temp2) then
  166.                temp = xreal(j)
  167.                xreal(j) = xreal(j+igap)
  168.                xreal(j+igap) = temp
  169. c
  170.                temp = ximag(j)
  171.                ximag(j) = ximag(j+igap)
  172.                ximag(j+igap) = temp
  173. c
  174.                if (apply) then
  175.                   temp = y(j)
  176.                   y(j) = y(j+igap)
  177.                   y(j+igap) = temp
  178.                end if
  179.             else
  180.                go to 60
  181.             endif
  182.             j = j-igap
  183.             go to 50
  184.    60    continue
  185.          igap = igap / 2
  186.          go to 40
  187. c
  188.       else if (which .eq. 'LR') then
  189. c
  190. c        %------------------------------------------------%
  191. c        | Sort XREAL into increasing order of algebraic. |
  192. c        %------------------------------------------------%
  193. c
  194.    70    continue
  195.          if (igap .eq. 0) go to 9000
  196. c
  197.          do 90 i = igap, n-1
  198.             j = i-igap
  199.    80       continue
  200. c
  201.             if (j.lt.0) go to 90
  202. c
  203.             if (xreal(j).gt.xreal(j+igap)) then
  204.                temp = xreal(j)
  205.                xreal(j) = xreal(j+igap)
  206.                xreal(j+igap) = temp
  207. c
  208.                temp = ximag(j)
  209.                ximag(j) = ximag(j+igap)
  210.                ximag(j+igap) = temp
  211. c
  212.                if (apply) then
  213.                   temp = y(j)
  214.                   y(j) = y(j+igap)
  215.                   y(j+igap) = temp
  216.                end if
  217.             else
  218.                go to 90
  219.             endif
  220.             j = j-igap
  221.             go to 80
  222.    90    continue
  223.          igap = igap / 2
  224.          go to 70
  225. c
  226.       else if (which .eq. 'SR') then
  227. c
  228. c        %------------------------------------------------%
  229. c        | Sort XREAL into decreasing order of algebraic. |
  230. c        %------------------------------------------------%
  231. c
  232.   100    continue
  233.          if (igap .eq. 0) go to 9000
  234.          do 120 i = igap, n-1
  235.             j = i-igap
  236.   110       continue
  237. c
  238.             if (j.lt.0) go to 120
  239. c
  240.             if (xreal(j).lt.xreal(j+igap)) then
  241.                temp = xreal(j)
  242.                xreal(j) = xreal(j+igap)
  243.                xreal(j+igap) = temp
  244. c
  245.                temp = ximag(j)
  246.                ximag(j) = ximag(j+igap)
  247.                ximag(j+igap) = temp
  248. c
  249.                if (apply) then
  250.                   temp = y(j)
  251.                   y(j) = y(j+igap)
  252.                   y(j+igap) = temp
  253.                end if
  254.             else
  255.                go to 120
  256.             endif
  257.             j = j-igap
  258.             go to 110
  259.   120    continue
  260.          igap = igap / 2
  261.          go to 100
  262. c
  263.       else if (which .eq. 'LI') then
  264. c
  265. c        %------------------------------------------------%
  266. c        | Sort XIMAG into increasing order of magnitude. |
  267. c        %------------------------------------------------%
  268. c
  269.   130    continue
  270.          if (igap .eq. 0) go to 9000
  271.          do 150 i = igap, n-1
  272.             j = i-igap
  273.   140       continue
  274. c
  275.             if (j.lt.0) go to 150
  276. c
  277.             if (abs(ximag(j)).gt.abs(ximag(j+igap))) then
  278.                temp = xreal(j)
  279.                xreal(j) = xreal(j+igap)
  280.                xreal(j+igap) = temp
  281. c
  282.                temp = ximag(j)
  283.                ximag(j) = ximag(j+igap)
  284.                ximag(j+igap) = temp
  285. c
  286.                if (apply) then
  287.                   temp = y(j)
  288.                   y(j) = y(j+igap)
  289.                   y(j+igap) = temp
  290.                end if
  291.             else
  292.                go to 150
  293.             endif
  294.             j = j-igap
  295.             go to 140
  296.   150    continue
  297.          igap = igap / 2
  298.          go to 130
  299. c
  300.       else if (which .eq. 'SI') then
  301. c
  302. c        %------------------------------------------------%
  303. c        | Sort XIMAG into decreasing order of magnitude. |
  304. c        %------------------------------------------------%
  305. c
  306.   160    continue
  307.          if (igap .eq. 0) go to 9000
  308.          do 180 i = igap, n-1
  309.             j = i-igap
  310.   170       continue
  311. c
  312.             if (j.lt.0) go to 180
  313. c
  314.             if (abs(ximag(j)).lt.abs(ximag(j+igap))) then
  315.                temp = xreal(j)
  316.                xreal(j) = xreal(j+igap)
  317.                xreal(j+igap) = temp
  318. c
  319.                temp = ximag(j)
  320.                ximag(j) = ximag(j+igap)
  321.                ximag(j+igap) = temp
  322. c
  323.                if (apply) then
  324.                   temp = y(j)
  325.                   y(j) = y(j+igap)
  326.                   y(j+igap) = temp
  327.                end if
  328.             else
  329.                go to 180
  330.             endif
  331.             j = j-igap
  332.             go to 170
  333.   180    continue
  334.          igap = igap / 2
  335.          go to 160
  336.       end if
  337. c
  338.  9000 continue
  339.       return
  340. c
  341. c     %---------------%
  342. c     | End of dsortc |
  343. c     %---------------%
  344. c
  345.       end
  346.  
  347.  

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