## "Fossies" - the Fresh Open Source Software Archive

### Member "getdp-3.1.0-source/contrib/Arpack/dsortc.f" (31 Jul 2018, 9254 Bytes) of package /linux/privat/getdp-3.1.0-source.tgz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Fortran 77 source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "dsortc.f" see the Fossies "Dox" file reference documentation.

1 c-----------------------------------------------------------------------
2 c\BeginDoc
3 c
4 c\Name: dsortc
5 c
6 c\Description:
7 c  Sorts the complex array in XREAL and XIMAG into the order
8 c  specified by WHICH and optionally applies the permutation to the
9 c  real array Y. It is assumed that if an element of XIMAG is
10 c  nonzero, then its negative is also an element. In other words,
11 c  both members of a complex conjugate pair are to be sorted and the
12 c  pairs are kept adjacent to each other.
13 c
14 c\Usage:
15 c  call dsortc
16 c     ( WHICH, APPLY, N, XREAL, XIMAG, Y )
17 c
18 c\Arguments
19 c  WHICH   Character*2.  (Input)
20 c          'LM' -> sort XREAL,XIMAG into increasing order of magnitude.
21 c          'SM' -> sort XREAL,XIMAG into decreasing order of magnitude.
22 c          'LR' -> sort XREAL into increasing order of algebraic.
23 c          'SR' -> sort XREAL into decreasing order of algebraic.
24 c          'LI' -> sort XIMAG into increasing order of magnitude.
25 c          'SI' -> sort XIMAG into decreasing order of magnitude.
26 c          NOTE: If an element of XIMAG is non-zero, then its negative
27 c                is also an element.
28 c
29 c  APPLY   Logical.  (Input)
30 c          APPLY = .TRUE.  -> apply the sorted order to array Y.
31 c          APPLY = .FALSE. -> do not apply the sorted order to array Y.
32 c
33 c  N       Integer.  (INPUT)
34 c          Size of the arrays.
35 c
36 c  XREAL,  Double precision array of length N.  (INPUT/OUTPUT)
37 c  XIMAG   Real and imaginary part of the array to be sorted.
38 c
39 c  Y       Double precision array of length N.  (INPUT/OUTPUT)
40 c
41 c\EndDoc
42 c
43 c-----------------------------------------------------------------------
44 c
45 c\BeginLib
46 c
47 c\Author
48 c     Danny Sorensen               Phuong Vu
49 c     Richard Lehoucq              CRPC / Rice University
50 c     Dept. of Computational &     Houston, Texas
51 c     Applied Mathematics
52 c     Rice University
53 c     Houston, Texas
54 c
55 c\Revision history:
56 c     xx/xx/92: Version ' 2.1'
57 c               Adapted from the sort routine in LANSO.
58 c
59 c\SCCS Information: @(#)
60 c FILE: sortc.F   SID: 2.3   DATE OF SID: 4/20/96   RELEASE: 2
61 c
62 c\EndLib
63 c
64 c-----------------------------------------------------------------------
65 c
66       subroutine dsortc (which, apply, n, xreal, ximag, y)
67 c
68 c     %------------------%
69 c     | Scalar Arguments |
70 c     %------------------%
71 c
72       character*2 which
73       logical    apply
74       integer    n
75 c
76 c     %-----------------%
77 c     | Array Arguments |
78 c     %-----------------%
79 c
80       Double precision
81      &           xreal(0:n-1), ximag(0:n-1), y(0:n-1)
82 c
83 c     %---------------%
84 c     | Local Scalars |
85 c     %---------------%
86 c
87       integer    i, igap, j
88       Double precision
89      &           temp, temp1, temp2
90 c
91 c     %--------------------%
92 c     | External Functions |
93 c     %--------------------%
94 c
95       Double precision
96      &           dlapy2
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