Mercurial > octave-nkf

comparison libcruft/lapack/csteqr.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 SUBROUTINE CSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
2 *
3 * -- LAPACK routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 CHARACTER COMPZ
9 INTEGER INFO, LDZ, N
10 * ..
11 * .. Array Arguments ..
12 REAL D( * ), E( * ), WORK( * )
13 COMPLEX Z( LDZ, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * CSTEQR computes all eigenvalues and, optionally, eigenvectors of a
20 * symmetric tridiagonal matrix using the implicit QL or QR method.
21 * The eigenvectors of a full or band complex Hermitian matrix can also
22 * be found if CHETRD or CHPTRD or CHBTRD has been used to reduce this
23 * matrix to tridiagonal form.
24 *
25 * Arguments
26 * =========
27 *
28 * COMPZ (input) CHARACTER*1
29 * = 'N': Compute eigenvalues only.
30 * = 'V': Compute eigenvalues and eigenvectors of the original
31 * Hermitian matrix. On entry, Z must contain the
32 * unitary matrix used to reduce the original matrix
33 * to tridiagonal form.
34 * = 'I': Compute eigenvalues and eigenvectors of the
35 * tridiagonal matrix. Z is initialized to the identity
36 * matrix.
37 *
38 * N (input) INTEGER
39 * The order of the matrix. N >= 0.
40 *
41 * D (input/output) REAL array, dimension (N)
42 * On entry, the diagonal elements of the tridiagonal matrix.
43 * On exit, if INFO = 0, the eigenvalues in ascending order.
44 *
45 * E (input/output) REAL array, dimension (N-1)
46 * On entry, the (n-1) subdiagonal elements of the tridiagonal
47 * matrix.
48 * On exit, E has been destroyed.
49 *
50 * Z (input/output) COMPLEX array, dimension (LDZ, N)
51 * On entry, if COMPZ = 'V', then Z contains the unitary
52 * matrix used in the reduction to tridiagonal form.
53 * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
54 * orthonormal eigenvectors of the original Hermitian matrix,
55 * and if COMPZ = 'I', Z contains the orthonormal eigenvectors
56 * of the symmetric tridiagonal matrix.
57 * If COMPZ = 'N', then Z is not referenced.
58 *
59 * LDZ (input) INTEGER
60 * The leading dimension of the array Z. LDZ >= 1, and if
61 * eigenvectors are desired, then LDZ >= max(1,N).
62 *
63 * WORK (workspace) REAL array, dimension (max(1,2*N-2))
64 * If COMPZ = 'N', then WORK is not referenced.
65 *
66 * INFO (output) INTEGER
67 * = 0: successful exit
68 * < 0: if INFO = -i, the i-th argument had an illegal value
69 * > 0: the algorithm has failed to find all the eigenvalues in
70 * a total of 30*N iterations; if INFO = i, then i
71 * elements of E have not converged to zero; on exit, D
72 * and E contain the elements of a symmetric tridiagonal
73 * matrix which is unitarily similar to the original
74 * matrix.
75 *
76 * =====================================================================
77 *
78 * .. Parameters ..
79 REAL ZERO, ONE, TWO, THREE
80 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
81 $ THREE = 3.0E0 )
82 COMPLEX CZERO, CONE
83 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ),
84 $ CONE = ( 1.0E0, 0.0E0 ) )
85 INTEGER MAXIT
86 PARAMETER ( MAXIT = 30 )
87 * ..
88 * .. Local Scalars ..
89 INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
90 $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
91 $ NM1, NMAXIT
92 REAL ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
93 $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
94 * ..
95 * .. External Functions ..
96 LOGICAL LSAME
97 REAL SLAMCH, SLANST, SLAPY2
98 EXTERNAL LSAME, SLAMCH, SLANST, SLAPY2
99 * ..
100 * .. External Subroutines ..
101 EXTERNAL CLASET, CLASR, CSWAP, SLAE2, SLAEV2, SLARTG,
102 $ SLASCL, SLASRT, XERBLA
103 * ..
104 * .. Intrinsic Functions ..
105 INTRINSIC ABS, MAX, SIGN, SQRT
106 * ..
107 * .. Executable Statements ..
108 *
109 * Test the input parameters.
110 *
111 INFO = 0
112 *
113 IF( LSAME( COMPZ, 'N' ) ) THEN
114 ICOMPZ = 0
115 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
116 ICOMPZ = 1
117 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
118 ICOMPZ = 2
119 ELSE
120 ICOMPZ = -1
121 END IF
122 IF( ICOMPZ.LT.0 ) THEN
123 INFO = -1
124 ELSE IF( N.LT.0 ) THEN
125 INFO = -2
126 ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
127 $ N ) ) ) THEN
128 INFO = -6
129 END IF
130 IF( INFO.NE.0 ) THEN
131 CALL XERBLA( 'CSTEQR', -INFO )
132 RETURN
133 END IF
134 *
135 * Quick return if possible
136 *
137 IF( N.EQ.0 )
138 $ RETURN
139 *
140 IF( N.EQ.1 ) THEN
141 IF( ICOMPZ.EQ.2 )
142 $ Z( 1, 1 ) = CONE
143 RETURN
144 END IF
145 *
146 * Determine the unit roundoff and over/underflow thresholds.
147 *
148 EPS = SLAMCH( 'E' )
149 EPS2 = EPS**2
150 SAFMIN = SLAMCH( 'S' )
151 SAFMAX = ONE / SAFMIN
152 SSFMAX = SQRT( SAFMAX ) / THREE
153 SSFMIN = SQRT( SAFMIN ) / EPS2
154 *
155 * Compute the eigenvalues and eigenvectors of the tridiagonal
156 * matrix.
157 *
158 IF( ICOMPZ.EQ.2 )
159 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
160 *
161 NMAXIT = N*MAXIT
162 JTOT = 0
163 *
164 * Determine where the matrix splits and choose QL or QR iteration
165 * for each block, according to whether top or bottom diagonal
166 * element is smaller.
167 *
168 L1 = 1
169 NM1 = N - 1
170 *
171 10 CONTINUE
172 IF( L1.GT.N )
173 $ GO TO 160
174 IF( L1.GT.1 )
175 $ E( L1-1 ) = ZERO
176 IF( L1.LE.NM1 ) THEN
177 DO 20 M = L1, NM1
178 TST = ABS( E( M ) )
179 IF( TST.EQ.ZERO )
180 $ GO TO 30
181 IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
182 $ 1 ) ) ) )*EPS ) THEN
183 E( M ) = ZERO
184 GO TO 30
185 END IF
186 20 CONTINUE
187 END IF
188 M = N
189 *
190 30 CONTINUE
191 L = L1
192 LSV = L
193 LEND = M
194 LENDSV = LEND
195 L1 = M + 1
196 IF( LEND.EQ.L )
197 $ GO TO 10
198 *
199 * Scale submatrix in rows and columns L to LEND
200 *
201 ANORM = SLANST( 'I', LEND-L+1, D( L ), E( L ) )
202 ISCALE = 0
203 IF( ANORM.EQ.ZERO )
204 $ GO TO 10
205 IF( ANORM.GT.SSFMAX ) THEN
206 ISCALE = 1
207 CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
208 $ INFO )
209 CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
210 $ INFO )
211 ELSE IF( ANORM.LT.SSFMIN ) THEN
212 ISCALE = 2
213 CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
214 $ INFO )
215 CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
216 $ INFO )
217 END IF
218 *
219 * Choose between QL and QR iteration
220 *
221 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
222 LEND = LSV
223 L = LENDSV
224 END IF
225 *
226 IF( LEND.GT.L ) THEN
227 *
228 * QL Iteration
229 *
230 * Look for small subdiagonal element.
231 *
232 40 CONTINUE
233 IF( L.NE.LEND ) THEN
234 LENDM1 = LEND - 1
235 DO 50 M = L, LENDM1
236 TST = ABS( E( M ) )**2
237 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
238 $ SAFMIN )GO TO 60
239 50 CONTINUE
240 END IF
241 *
242 M = LEND
243 *
244 60 CONTINUE
245 IF( M.LT.LEND )
246 $ E( M ) = ZERO
247 P = D( L )
248 IF( M.EQ.L )
249 $ GO TO 80
250 *
251 * If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
252 * to compute its eigensystem.
253 *
254 IF( M.EQ.L+1 ) THEN
255 IF( ICOMPZ.GT.0 ) THEN
256 CALL SLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
257 WORK( L ) = C
258 WORK( N-1+L ) = S
259 CALL CLASR( 'R', 'V', 'B', N, 2, WORK( L ),
260 $ WORK( N-1+L ), Z( 1, L ), LDZ )
261 ELSE
262 CALL SLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 )
263 END IF
264 D( L ) = RT1
265 D( L+1 ) = RT2
266 E( L ) = ZERO
267 L = L + 2
268 IF( L.LE.LEND )
269 $ GO TO 40
270 GO TO 140
271 END IF
272 *
273 IF( JTOT.EQ.NMAXIT )
274 $ GO TO 140
275 JTOT = JTOT + 1
276 *
277 * Form shift.
278 *
279 G = ( D( L+1 )-P ) / ( TWO*E( L ) )
280 R = SLAPY2( G, ONE )
281 G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
282 *
283 S = ONE
284 C = ONE
285 P = ZERO
286 *
287 * Inner loop
288 *
289 MM1 = M - 1
290 DO 70 I = MM1, L, -1
291 F = S*E( I )
292 B = C*E( I )
293 CALL SLARTG( G, F, C, S, R )
294 IF( I.NE.M-1 )
295 $ E( I+1 ) = R
296 G = D( I+1 ) - P
297 R = ( D( I )-G )*S + TWO*C*B
298 P = S*R
299 D( I+1 ) = G + P
300 G = C*R - B
301 *
302 * If eigenvectors are desired, then save rotations.
303 *
304 IF( ICOMPZ.GT.0 ) THEN
305 WORK( I ) = C
306 WORK( N-1+I ) = -S
307 END IF
308 *
309 70 CONTINUE
310 *
311 * If eigenvectors are desired, then apply saved rotations.
312 *
313 IF( ICOMPZ.GT.0 ) THEN
314 MM = M - L + 1
315 CALL CLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ),
316 $ Z( 1, L ), LDZ )
317 END IF
318 *
319 D( L ) = D( L ) - P
320 E( L ) = G
321 GO TO 40
322 *
323 * Eigenvalue found.
324 *
325 80 CONTINUE
326 D( L ) = P
327 *
328 L = L + 1
329 IF( L.LE.LEND )
330 $ GO TO 40
331 GO TO 140
332 *
333 ELSE
334 *
335 * QR Iteration
336 *
337 * Look for small superdiagonal element.
338 *
339 90 CONTINUE
340 IF( L.NE.LEND ) THEN
341 LENDP1 = LEND + 1
342 DO 100 M = L, LENDP1, -1
343 TST = ABS( E( M-1 ) )**2
344 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
345 $ SAFMIN )GO TO 110
346 100 CONTINUE
347 END IF
348 *
349 M = LEND
350 *
351 110 CONTINUE
352 IF( M.GT.LEND )
353 $ E( M-1 ) = ZERO
354 P = D( L )
355 IF( M.EQ.L )
356 $ GO TO 130
357 *
358 * If remaining matrix is 2-by-2, use SLAE2 or SLAEV2
359 * to compute its eigensystem.
360 *
361 IF( M.EQ.L-1 ) THEN
362 IF( ICOMPZ.GT.0 ) THEN
363 CALL SLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
364 WORK( M ) = C
365 WORK( N-1+M ) = S
366 CALL CLASR( 'R', 'V', 'F', N, 2, WORK( M ),
367 $ WORK( N-1+M ), Z( 1, L-1 ), LDZ )
368 ELSE
369 CALL SLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
370 END IF
371 D( L-1 ) = RT1
372 D( L ) = RT2
373 E( L-1 ) = ZERO
374 L = L - 2
375 IF( L.GE.LEND )
376 $ GO TO 90
377 GO TO 140
378 END IF
379 *
380 IF( JTOT.EQ.NMAXIT )
381 $ GO TO 140
382 JTOT = JTOT + 1
383 *
384 * Form shift.
385 *
386 G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
387 R = SLAPY2( G, ONE )
388 G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
389 *
390 S = ONE
391 C = ONE
392 P = ZERO
393 *
394 * Inner loop
395 *
396 LM1 = L - 1
397 DO 120 I = M, LM1
398 F = S*E( I )
399 B = C*E( I )
400 CALL SLARTG( G, F, C, S, R )
401 IF( I.NE.M )
402 $ E( I-1 ) = R
403 G = D( I ) - P
404 R = ( D( I+1 )-G )*S + TWO*C*B
405 P = S*R
406 D( I ) = G + P
407 G = C*R - B
408 *
409 * If eigenvectors are desired, then save rotations.
410 *
411 IF( ICOMPZ.GT.0 ) THEN
412 WORK( I ) = C
413 WORK( N-1+I ) = S
414 END IF
415 *
416 120 CONTINUE
417 *
418 * If eigenvectors are desired, then apply saved rotations.
419 *
420 IF( ICOMPZ.GT.0 ) THEN
421 MM = L - M + 1
422 CALL CLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ),
423 $ Z( 1, M ), LDZ )
424 END IF
425 *
426 D( L ) = D( L ) - P
427 E( LM1 ) = G
428 GO TO 90
429 *
430 * Eigenvalue found.
431 *
432 130 CONTINUE
433 D( L ) = P
434 *
435 L = L - 1
436 IF( L.GE.LEND )
437 $ GO TO 90
438 GO TO 140
439 *
440 END IF
441 *
442 * Undo scaling if necessary
443 *
444 140 CONTINUE
445 IF( ISCALE.EQ.1 ) THEN
446 CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
447 $ D( LSV ), N, INFO )
448 CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
449 $ N, INFO )
450 ELSE IF( ISCALE.EQ.2 ) THEN
451 CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
452 $ D( LSV ), N, INFO )
453 CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
454 $ N, INFO )
455 END IF
456 *
457 * Check for no convergence to an eigenvalue after a total
458 * of N*MAXIT iterations.
459 *
460 IF( JTOT.EQ.NMAXIT ) THEN
461 DO 150 I = 1, N - 1
462 IF( E( I ).NE.ZERO )
463 $ INFO = INFO + 1
464 150 CONTINUE
465 RETURN
466 END IF
467 GO TO 10
468 *
469 * Order eigenvalues and eigenvectors.
470 *
471 160 CONTINUE
472 IF( ICOMPZ.EQ.0 ) THEN
473 *
474 * Use Quick Sort
475 *
476 CALL SLASRT( 'I', N, D, INFO )
477 *
478 ELSE
479 *
480 * Use Selection Sort to minimize swaps of eigenvectors
481 *
482 DO 180 II = 2, N
483 I = II - 1
484 K = I
485 P = D( I )
486 DO 170 J = II, N
487 IF( D( J ).LT.P ) THEN
488 K = J
489 P = D( J )
490 END IF
491 170 CONTINUE
492 IF( K.NE.I ) THEN
493 D( K ) = D( I )
494 D( I ) = P
495 CALL CSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
496 END IF
497 180 CONTINUE
498 END IF
499 RETURN
500 *
501 * End of CSTEQR
502 *
503 END