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