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