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