Mercurial > octave
comparison libcruft/arpack/src/dsortc.f @ 12274:9f5d2ef078e8 release-3-4-x
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: 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 |