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*> \date December 2016 129* 130*> \ingroup complexOTHERcomputational 131* 132* ===================================================================== 133 SUBROUTINE CSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO ) 134* 135* -- LAPACK computational routine (version 3.7.0) -- 136* -- LAPACK is a software package provided by Univ. of Tennessee, -- 137* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 138* December 2016 139* 140* .. Scalar Arguments .. 141 CHARACTER COMPZ 142 INTEGER INFO, LDZ, N 143* .. 144* .. Array Arguments .. 145 REAL D( * ), E( * ), WORK( * ) 146 COMPLEX Z( LDZ, * ) 147* .. 148* 149* ===================================================================== 150* 151* .. Parameters .. 152 REAL ZERO, ONE, TWO, THREE 153 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, 154 $ THREE = 3.0E0 ) 155 COMPLEX CZERO, CONE 156 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ), 157 $ CONE = ( 1.0E0, 0.0E0 ) ) 158 INTEGER MAXIT 159 PARAMETER ( MAXIT = 30 ) 160* .. 161* .. Local Scalars .. 162 INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND, 163 $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1, 164 $ NM1, NMAXIT 165 REAL ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2, 166 $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST 167* .. 168* .. External Functions .. 169 LOGICAL LSAME 170 REAL SLAMCH, SLANST, SLAPY2 171 EXTERNAL LSAME, SLAMCH, SLANST, SLAPY2 172* .. 173* .. External Subroutines .. 174 EXTERNAL CLASET, CLASR, CSWAP, SLAE2, SLAEV2, SLARTG, 175 $ SLASCL, SLASRT, XERBLA 176* .. 177* .. Intrinsic Functions .. 178 INTRINSIC ABS, MAX, SIGN, SQRT 179* .. 180* .. Executable Statements .. 181* 182* Test the input parameters. 183* 184 INFO = 0 185* 186 IF( LSAME( COMPZ, 'N' ) ) THEN 187 ICOMPZ = 0 188 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 189 ICOMPZ = 1 190 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 191 ICOMPZ = 2 192 ELSE 193 ICOMPZ = -1 194 END IF 195 IF( ICOMPZ.LT.0 ) THEN 196 INFO = -1 197 ELSE IF( N.LT.0 ) THEN 198 INFO = -2 199 ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, 200 $ N ) ) ) THEN 201 INFO = -6 202 END IF 203 IF( INFO.NE.0 ) THEN 204 CALL XERBLA( 'CSTEQR', -INFO ) 205 RETURN 206 END IF 207* 208* Quick return if possible 209* 210 IF( N.EQ.0 ) 211 $ RETURN 212* 213 IF( N.EQ.1 ) THEN 214 IF( ICOMPZ.EQ.2 ) 215 $ Z( 1, 1 ) = CONE 216 RETURN 217 END IF 218* 219* Determine the unit roundoff and over/underflow thresholds. 220* 221 EPS = SLAMCH( 'E' ) 222 EPS2 = EPS**2 223 SAFMIN = SLAMCH( 'S' ) 224 SAFMAX = ONE / SAFMIN 225 SSFMAX = SQRT( SAFMAX ) / THREE 226 SSFMIN = SQRT( SAFMIN ) / EPS2 227* 228* Compute the eigenvalues and eigenvectors of the tridiagonal 229* matrix. 230* 231 IF( ICOMPZ.EQ.2 ) 232 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) 233* 234 NMAXIT = N*MAXIT 235 JTOT = 0 236* 237* Determine where the matrix splits and choose QL or QR iteration 238* for each block, according to whether top or bottom diagonal 239* element is smaller. 240* 241 L1 = 1 242 NM1 = N - 1 243* 244 10 CONTINUE 245 IF( L1.GT.N ) 246 $ GO TO 160 247 IF( L1.GT.1 ) 248 $ E( L1-1 ) = ZERO 249 IF( L1.LE.NM1 ) THEN 250 DO 20 M = L1, NM1 251 TST = ABS( E( M ) ) 252 IF( TST.EQ.ZERO ) 253 $ GO TO 30 254 IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ 255 $ 1 ) ) ) )*EPS ) THEN 256 E( M ) = ZERO 257 GO TO 30 258 END IF 259 20 CONTINUE 260 END IF 261 M = N 262* 263 30 CONTINUE 264 L = L1 265 LSV = L 266 LEND = M 267 LENDSV = LEND 268 L1 = M + 1 269 IF( LEND.EQ.L ) 270 $ GO TO 10 271* 272* Scale submatrix in rows and columns L to LEND 273* 274 ANORM = SLANST( 'I', LEND-L+1, D( L ), E( L ) ) 275 ISCALE = 0 276 IF( ANORM.EQ.ZERO ) 277 $ GO TO 10 278 IF( ANORM.GT.SSFMAX ) THEN 279 ISCALE = 1 280 CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, 281 $ INFO ) 282 CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, 283 $ INFO ) 284 ELSE IF( ANORM.LT.SSFMIN ) THEN 285 ISCALE = 2 286 CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, 287 $ INFO ) 288 CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, 289 $ INFO ) 290 END IF 291* 292* Choose between QL and QR iteration 293* 294 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN 295 LEND = LSV 296 L = LENDSV 297 END IF 298* 299 IF( LEND.GT.L ) THEN 300* 301* QL Iteration 302* 303* Look for small subdiagonal element. 304* 305 40 CONTINUE 306 IF( L.NE.LEND ) THEN 307 LENDM1 = LEND - 1 308 DO 50 M = L, LENDM1 309 TST = ABS( E( M ) )**2 310 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+ 311 $ SAFMIN )GO TO 60 312 50 CONTINUE 313 END IF 314* 315 M = LEND 316* 317 60 CONTINUE 318 IF( M.LT.LEND ) 319 $ E( M ) = ZERO 320 P = D( L ) 321 IF( M.EQ.L ) 322 $ GO TO 80 323* 324* If remaining matrix is 2-by-2, use SLAE2 or SLAEV2 325* to compute its eigensystem. 326* 327 IF( M.EQ.L+1 ) THEN 328 IF( ICOMPZ.GT.0 ) THEN 329 CALL SLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S ) 330 WORK( L ) = C 331 WORK( N-1+L ) = S 332 CALL CLASR( 'R', 'V', 'B', N, 2, WORK( L ), 333 $ WORK( N-1+L ), Z( 1, L ), LDZ ) 334 ELSE 335 CALL SLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 ) 336 END IF 337 D( L ) = RT1 338 D( L+1 ) = RT2 339 E( L ) = ZERO 340 L = L + 2 341 IF( L.LE.LEND ) 342 $ GO TO 40 343 GO TO 140 344 END IF 345* 346 IF( JTOT.EQ.NMAXIT ) 347 $ GO TO 140 348 JTOT = JTOT + 1 349* 350* Form shift. 351* 352 G = ( D( L+1 )-P ) / ( TWO*E( L ) ) 353 R = SLAPY2( G, ONE ) 354 G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) ) 355* 356 S = ONE 357 C = ONE 358 P = ZERO 359* 360* Inner loop 361* 362 MM1 = M - 1 363 DO 70 I = MM1, L, -1 364 F = S*E( I ) 365 B = C*E( I ) 366 CALL SLARTG( G, F, C, S, R ) 367 IF( I.NE.M-1 ) 368 $ E( I+1 ) = R 369 G = D( I+1 ) - P 370 R = ( D( I )-G )*S + TWO*C*B 371 P = S*R 372 D( I+1 ) = G + P 373 G = C*R - B 374* 375* If eigenvectors are desired, then save rotations. 376* 377 IF( ICOMPZ.GT.0 ) THEN 378 WORK( I ) = C 379 WORK( N-1+I ) = -S 380 END IF 381* 382 70 CONTINUE 383* 384* If eigenvectors are desired, then apply saved rotations. 385* 386 IF( ICOMPZ.GT.0 ) THEN 387 MM = M - L + 1 388 CALL CLASR( 'R', 'V', 'B', N, MM, WORK( L ), WORK( N-1+L ), 389 $ Z( 1, L ), LDZ ) 390 END IF 391* 392 D( L ) = D( L ) - P 393 E( L ) = G 394 GO TO 40 395* 396* Eigenvalue found. 397* 398 80 CONTINUE 399 D( L ) = P 400* 401 L = L + 1 402 IF( L.LE.LEND ) 403 $ GO TO 40 404 GO TO 140 405* 406 ELSE 407* 408* QR Iteration 409* 410* Look for small superdiagonal element. 411* 412 90 CONTINUE 413 IF( L.NE.LEND ) THEN 414 LENDP1 = LEND + 1 415 DO 100 M = L, LENDP1, -1 416 TST = ABS( E( M-1 ) )**2 417 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+ 418 $ SAFMIN )GO TO 110 419 100 CONTINUE 420 END IF 421* 422 M = LEND 423* 424 110 CONTINUE 425 IF( M.GT.LEND ) 426 $ E( M-1 ) = ZERO 427 P = D( L ) 428 IF( M.EQ.L ) 429 $ GO TO 130 430* 431* If remaining matrix is 2-by-2, use SLAE2 or SLAEV2 432* to compute its eigensystem. 433* 434 IF( M.EQ.L-1 ) THEN 435 IF( ICOMPZ.GT.0 ) THEN 436 CALL SLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S ) 437 WORK( M ) = C 438 WORK( N-1+M ) = S 439 CALL CLASR( 'R', 'V', 'F', N, 2, WORK( M ), 440 $ WORK( N-1+M ), Z( 1, L-1 ), LDZ ) 441 ELSE 442 CALL SLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 ) 443 END IF 444 D( L-1 ) = RT1 445 D( L ) = RT2 446 E( L-1 ) = ZERO 447 L = L - 2 448 IF( L.GE.LEND ) 449 $ GO TO 90 450 GO TO 140 451 END IF 452* 453 IF( JTOT.EQ.NMAXIT ) 454 $ GO TO 140 455 JTOT = JTOT + 1 456* 457* Form shift. 458* 459 G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ) 460 R = SLAPY2( G, ONE ) 461 G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) ) 462* 463 S = ONE 464 C = ONE 465 P = ZERO 466* 467* Inner loop 468* 469 LM1 = L - 1 470 DO 120 I = M, LM1 471 F = S*E( I ) 472 B = C*E( I ) 473 CALL SLARTG( G, F, C, S, R ) 474 IF( I.NE.M ) 475 $ E( I-1 ) = R 476 G = D( I ) - P 477 R = ( D( I+1 )-G )*S + TWO*C*B 478 P = S*R 479 D( I ) = G + P 480 G = C*R - B 481* 482* If eigenvectors are desired, then save rotations. 483* 484 IF( ICOMPZ.GT.0 ) THEN 485 WORK( I ) = C 486 WORK( N-1+I ) = S 487 END IF 488* 489 120 CONTINUE 490* 491* If eigenvectors are desired, then apply saved rotations. 492* 493 IF( ICOMPZ.GT.0 ) THEN 494 MM = L - M + 1 495 CALL CLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ), 496 $ Z( 1, M ), LDZ ) 497 END IF 498* 499 D( L ) = D( L ) - P 500 E( LM1 ) = G 501 GO TO 90 502* 503* Eigenvalue found. 504* 505 130 CONTINUE 506 D( L ) = P 507* 508 L = L - 1 509 IF( L.GE.LEND ) 510 $ GO TO 90 511 GO TO 140 512* 513 END IF 514* 515* Undo scaling if necessary 516* 517 140 CONTINUE 518 IF( ISCALE.EQ.1 ) THEN 519 CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, 520 $ D( LSV ), N, INFO ) 521 CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ), 522 $ N, INFO ) 523 ELSE IF( ISCALE.EQ.2 ) THEN 524 CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, 525 $ D( LSV ), N, INFO ) 526 CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ), 527 $ N, INFO ) 528 END IF 529* 530* Check for no convergence to an eigenvalue after a total 531* of N*MAXIT iterations. 532* 533 IF( JTOT.EQ.NMAXIT ) THEN 534 DO 150 I = 1, N - 1 535 IF( E( I ).NE.ZERO ) 536 $ INFO = INFO + 1 537 150 CONTINUE 538 RETURN 539 END IF 540 GO TO 10 541* 542* Order eigenvalues and eigenvectors. 543* 544 160 CONTINUE 545 IF( ICOMPZ.EQ.0 ) THEN 546* 547* Use Quick Sort 548* 549 CALL SLASRT( 'I', N, D, INFO ) 550* 551 ELSE 552* 553* Use Selection Sort to minimize swaps of eigenvectors 554* 555 DO 180 II = 2, N 556 I = II - 1 557 K = I 558 P = D( I ) 559 DO 170 J = II, N 560 IF( D( J ).LT.P ) THEN 561 K = J 562 P = D( J ) 563 END IF 564 170 CONTINUE 565 IF( K.NE.I ) THEN 566 D( K ) = D( I ) 567 D( I ) = P 568 CALL CSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 ) 569 END IF 570 180 CONTINUE 571 END IF 572 RETURN 573* 574* End of CSTEQR 575* 576 END 577