Mercurial > octave

comparison libcruft/arpack/src/ssortc.f @ 12274:9f5d2ef078e8 release-3-4-x

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
import ARPACK sources to libcruft from Debian package libarpack2 2.1+parpack96.dfsg-3+b1
author John W. Eaton <jwe@octave.org>
date Fri, 28 Jan 2011 14:04:33 -0500
parents
children
comparison
equal deleted inserted replaced
12273:83133b5bf392 12274:9f5d2ef078e8
1 c-----------------------------------------------------------------------
2 c\BeginDoc
3 c
4 c\Name: ssortc
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 ssortc
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, Real array of length N. (INPUT/OUTPUT)
37 c XIMAG Real and imaginary part of the array to be sorted.
38 c
39 c Y Real 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 ssortc (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 Real
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 Real
89 & temp, temp1, temp2
90 c
91 c %--------------------%
92 c | External Functions |
93 c %--------------------%
94 c
95 Real
96 & slapy2
97 external slapy2
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 = slapy2(xreal(j),ximag(j))
121 temp2 = slapy2(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 = slapy2(xreal(j),ximag(j))
162 temp2 = slapy2(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 ssortc |
342 c %---------------%
343 c
344 end