Mercurial > octave-nkf

comparison libcruft/arpack/src/csortc.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\BeginDoc
2 c
3 c\Name: csortc
4 c
5 c\Description:
6 c Sorts the Complex array in X into the order
7 c specified by WHICH and optionally applies the permutation to the
8 c Real array Y.
9 c
10 c\Usage:
11 c call csortc
12 c ( WHICH, APPLY, N, X, Y )
13 c
14 c\Arguments
15 c WHICH Character*2. (Input)
16 c 'LM' -> sort X into increasing order of magnitude.
17 c 'SM' -> sort X into decreasing order of magnitude.
18 c 'LR' -> sort X with real(X) in increasing algebraic order
19 c 'SR' -> sort X with real(X) in decreasing algebraic order
20 c 'LI' -> sort X with imag(X) in increasing algebraic order
21 c 'SI' -> sort X with imag(X) in decreasing algebraic order
22 c
23 c APPLY Logical. (Input)
24 c APPLY = .TRUE. -> apply the sorted order to array Y.
25 c APPLY = .FALSE. -> do not apply the sorted order to array Y.
26 c
27 c N Integer. (INPUT)
28 c Size of the arrays.
29 c
30 c X Complex array of length N. (INPUT/OUTPUT)
31 c This is the array to be sorted.
32 c
33 c Y Complex array of length N. (INPUT/OUTPUT)
34 c
35 c\EndDoc
36 c
37 c-----------------------------------------------------------------------
38 c
39 c\BeginLib
40 c
41 c\Routines called:
42 c slapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully.
43 c
44 c\Author
45 c Danny Sorensen Phuong Vu
46 c Richard Lehoucq CRPC / Rice University
47 c Dept. of Computational & Houston, Texas
48 c Applied Mathematics
49 c Rice University
50 c Houston, Texas
51 c
52 c Adapted from the sort routine in LANSO.
53 c
54 c\SCCS Information: @(#)
55 c FILE: sortc.F SID: 2.2 DATE OF SID: 4/20/96 RELEASE: 2
56 c
57 c\EndLib
58 c
59 c-----------------------------------------------------------------------
60 c
61 subroutine csortc (which, apply, n, x, y)
62 c
63 c %------------------%
64 c | Scalar Arguments |
65 c %------------------%
66 c
67 character*2 which
68 logical apply
69 integer n
70 c
71 c %-----------------%
72 c | Array Arguments |
73 c %-----------------%
74 c
75 Complex
76 & x(0:n-1), y(0:n-1)
77 c
78 c %---------------%
79 c | Local Scalars |
80 c %---------------%
81 c
82 integer i, igap, j
83 Complex
84 & temp
85 Real
86 & temp1, temp2
87 c
88 c %--------------------%
89 c | External functions |
90 c %--------------------%
91 c
92 Real
93 & slapy2
94 c
95 c %--------------------%
96 c | Intrinsic Functions |
97 c %--------------------%
98 Intrinsic
99 & real, aimag
100 c
101 c %-----------------------%
102 c | Executable Statements |
103 c %-----------------------%
104 c
105 igap = n / 2
106 c
107 if (which .eq. 'LM') then
108 c
109 c %--------------------------------------------%
110 c | Sort X into increasing order of magnitude. |
111 c %--------------------------------------------%
112 c
113 10 continue
114 if (igap .eq. 0) go to 9000
115 c
116 do 30 i = igap, n-1
117 j = i-igap
118 20 continue
119 c
120 if (j.lt.0) go to 30
121 c
122 temp1 = slapy2(real(x(j)),aimag(x(j)))
123 temp2 = slapy2(real(x(j+igap)),aimag(x(j+igap)))
124 c
125 if (temp1.gt.temp2) then
126 temp = x(j)
127 x(j) = x(j+igap)
128 x(j+igap) = temp
129 c
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
143 c
144 else if (which .eq. 'SM') then
145 c
146 c %--------------------------------------------%
147 c | Sort X into decreasing order of magnitude. |
148 c %--------------------------------------------%
149 c
150 40 continue
151 if (igap .eq. 0) go to 9000
152 c
153 do 60 i = igap, n-1
154 j = i-igap
155 50 continue
156 c
157 if (j .lt. 0) go to 60
158 c
159 temp1 = slapy2(real(x(j)),aimag(x(j)))
160 temp2 = slapy2(real(x(j+igap)),aimag(x(j+igap)))
161 c
162 if (temp1.lt.temp2) then
163 temp = x(j)
164 x(j) = x(j+igap)
165 x(j+igap) = temp
166 c
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
180 c
181 else if (which .eq. 'LR') then
182 c
183 c %------------------------------------------------%
184 c | Sort XREAL into increasing order of algebraic. |
185 c %------------------------------------------------%
186 c
187 70 continue
188 if (igap .eq. 0) go to 9000
189 c
190 do 90 i = igap, n-1
191 j = i-igap
192 80 continue
193 c
194 if (j.lt.0) go to 90
195 c
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
200 c
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
214 c
215 else if (which .eq. 'SR') then
216 c
217 c %------------------------------------------------%
218 c | Sort XREAL into decreasing order of algebraic. |
219 c %------------------------------------------------%
220 c
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
226 c
227 if (j.lt.0) go to 120
228 c
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
233 c
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
247 c
248 else if (which .eq. 'LI') then
249 c
250 c %--------------------------------------------%
251 c | Sort XIMAG into increasing algebraic order |
252 c %--------------------------------------------%
253 c
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
259 c
260 if (j.lt.0) go to 150
261 c
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
266 c
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
280 c
281 else if (which .eq. 'SI') then
282 c
283 c %---------------------------------------------%
284 c | Sort XIMAG into decreasing algebraic order |
285 c %---------------------------------------------%
286 c
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
292 c
293 if (j.lt.0) go to 180
294 c
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
299 c
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
314 c
315 9000 continue
316 return
317 c
318 c %---------------%
319 c | End of csortc |
320 c %---------------%
321 c
322 end