1 SUBROUTINE CSTEQR2( COMPZ, N, D, E, Z, LDZ, NR, WORK, INFO ) 2* 3* -- ScaLAPACK routine (version 1.7) -- 4* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5* Courant Institute, Argonne National Lab, and Rice University 6* November 15, 1997 7* 8* .. Scalar Arguments .. 9 CHARACTER COMPZ 10 INTEGER INFO, LDZ, N, NR 11* .. 12* .. Array Arguments .. 13 REAL D( * ), E( * ), WORK( * ) 14 COMPLEX Z( LDZ, * ) 15* .. 16* 17* Purpose 18* ======= 19* 20* CSTEQR2 is a modified version of LAPACK routine CSTEQR. 21* CSTEQR2 computes all eigenvalues and, optionally, eigenvectors of a 22* symmetric tridiagonal matrix using the implicit QL or QR method. 23* CSTEQR2 is modified from CSTEQR to allow each ScaLAPACK process 24* running CSTEQR2 to perform updates on a distributed matrix Q. 25* Proper usage of CSTEQR2 can be gleaned from 26* examination of ScaLAPACK's * PCHEEV. 27* CSTEQR2 incorporates changes attributed to Greg Henry. 28* 29* Arguments 30* ========= 31* 32* COMPZ (input) CHARACTER*1 33* = 'N': Compute eigenvalues only. 34* = 'I': Compute eigenvalues and eigenvectors of the 35* tridiagonal matrix. Z must be initialized to the 36* identity matrix by PCLASET or CLASET prior 37* to entering this subroutine. 38* 39* N (input) INTEGER 40* The order of the matrix. N >= 0. 41* 42* D (input/output) REAL array, dimension (N) 43* On entry, the diagonal elements of the tridiagonal matrix. 44* On exit, if INFO = 0, the eigenvalues in ascending order. 45* 46* E (input/output) REAL array, dimension (N-1) 47* On entry, the (n-1) subdiagonal elements of the tridiagonal 48* matrix. 49* On exit, E has been destroyed. 50* 51* Z (local input/local output) COMPLEX array, global 52* dimension (N, N), local dimension (LDZ, NR). 53* On entry, if COMPZ = 'V', then Z contains the orthogonal 54* matrix used in the reduction to tridiagonal form. 55* On exit, if INFO = 0, then if COMPZ = 'V', Z contains the 56* orthonormal eigenvectors of the original symmetric matrix, 57* and if COMPZ = 'I', Z contains the orthonormal eigenvectors 58* of the symmetric tridiagonal matrix. 59* If COMPZ = 'N', then Z is not referenced. 60* 61* LDZ (input) INTEGER 62* The leading dimension of the array Z. LDZ >= 1, and if 63* eigenvectors are desired, then LDZ >= max(1,N). 64* 65* NR (input) INTEGER 66* NR = MAX(1, NUMROC( N, NB, MYPROW, 0, NPROCS ) ). 67* If COMPZ = 'N', then NR is not referenced. 68* 69* WORK (workspace) REAL array, dimension (max(1,2*N-2)) 70* If COMPZ = 'N', then WORK is not referenced. 71* 72* INFO (output) INTEGER 73* = 0: successful exit 74* < 0: if INFO = -i, the i-th argument had an illegal value 75* > 0: the algorithm has failed to find all the eigenvalues in 76* a total of 30*N iterations; if INFO = i, then i 77* elements of E have not converged to zero; on exit, D 78* and E contain the elements of a symmetric tridiagonal 79* matrix which is orthogonally similar to the original 80* matrix. 81* 82* ===================================================================== 83* 84* .. Parameters .. 85 REAL ZERO, ONE, TWO, THREE, HALF 86 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, 87 $ THREE = 3.0E0, HALF = 0.5E0 ) 88 COMPLEX CONE 89 PARAMETER ( CONE = ( 1.0E0, 1.0E0 ) ) 90 INTEGER MAXIT, NMAXLOOK 91 PARAMETER ( MAXIT = 30, NMAXLOOK = 15 ) 92* .. 93* .. Local Scalars .. 94 INTEGER I, ICOMPZ, II, ILAST, ISCALE, J, JTOT, K, L, 95 $ L1, LEND, LENDM1, LENDP1, LENDSV, LM1, LSV, M, 96 $ MM, MM1, NLOOK, NM1, NMAXIT 97 REAL ANORM, B, C, EPS, EPS2, F, G, GP, OLDEL, OLDGP, 98 $ OLDRP, P, R, RP, RT1, RT2, S, SAFMAX, SAFMIN, 99 $ SSFMAX, SSFMIN, TST, TST1 100* .. 101* .. External Functions .. 102 LOGICAL LSAME 103 REAL SLAMCH, SLANST, SLAPY2 104 EXTERNAL LSAME, SLAMCH, SLANST, SLAPY2 105* .. 106* .. External Subroutines .. 107 EXTERNAL CLASR, CSWAP, SLAEV2, SLARTG, SLASCL, SSTERF, 108 $ XERBLA 109* .. 110* .. Intrinsic Functions .. 111 INTRINSIC ABS, MAX, SIGN, SQRT 112* .. 113* .. Executable Statements .. 114* 115* Test the input parameters. 116* 117 ILAST = 0 118 INFO = 0 119* 120 IF( LSAME( COMPZ, 'N' ) ) THEN 121 ICOMPZ = 0 122 ELSEIF( LSAME( COMPZ, 'I' ) ) THEN 123 ICOMPZ = 1 124 ELSE 125 ICOMPZ = -1 126 ENDIF 127 IF( ICOMPZ.LT.0 ) THEN 128 INFO = -1 129 ELSEIF( N.LT.0 ) THEN 130 INFO = -2 131 ELSEIF( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, NR ) ) THEN 132 INFO = -6 133 ENDIF 134 IF( INFO.NE.0 ) THEN 135 CALL XERBLA( 'CSTEQR2', -INFO ) 136 RETURN 137 ENDIF 138* 139* Quick return if possible 140* 141 IF( N.EQ.0 ) 142 $ RETURN 143* 144* If eigenvectors aren't not desired, this is faster 145* 146 IF( ICOMPZ.EQ.0 ) THEN 147 CALL SSTERF( N, D, E, INFO ) 148 RETURN 149 ENDIF 150* 151 IF( N.EQ.1 ) THEN 152 Z( 1, 1 ) = CONE 153 RETURN 154 ENDIF 155* 156* Determine the unit roundoff and over/underflow thresholds. 157* 158 EPS = SLAMCH( 'E' ) 159 EPS2 = EPS**2 160 SAFMIN = SLAMCH( 'S' ) 161 SAFMAX = ONE / SAFMIN 162 SSFMAX = SQRT( SAFMAX ) / THREE 163 SSFMIN = SQRT( SAFMIN ) / EPS2 164* 165* Compute the eigenvalues and eigenvectors of the tridiagonal 166* matrix. 167* 168 NMAXIT = N*MAXIT 169 JTOT = 0 170* 171* Determine where the matrix splits and choose QL or QR iteration 172* for each block, according to whether top or bottom diagonal 173* element is smaller. 174* 175 L1 = 1 176 NM1 = N - 1 177* 178 10 CONTINUE 179 IF( L1.GT.N ) 180 $ GOTO 220 181 IF( L1.GT.1 ) 182 $ E( L1-1 ) = ZERO 183 IF( L1.LE.NM1 ) THEN 184 DO 20 M = L1, NM1 185 TST = ABS( E( M ) ) 186 IF( TST.EQ.ZERO ) 187 $ GOTO 30 188 IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ 189 $ 1 ) ) ) )*EPS ) THEN 190 E( M ) = ZERO 191 GOTO 30 192 ENDIF 193 20 CONTINUE 194 ENDIF 195 M = N 196* 197 30 CONTINUE 198 L = L1 199 LSV = L 200 LEND = M 201 LENDSV = LEND 202 L1 = M + 1 203 IF( LEND.EQ.L ) 204 $ GOTO 10 205* 206* Scale submatrix in rows and columns L to LEND 207* 208 ANORM = SLANST( 'I', LEND-L+1, D( L ), E( L ) ) 209 ISCALE = 0 210 IF( ANORM.EQ.ZERO ) 211 $ GOTO 10 212 IF( ANORM.GT.SSFMAX ) THEN 213 ISCALE = 1 214 CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, 215 $ INFO ) 216 CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, 217 $ INFO ) 218 ELSEIF( ANORM.LT.SSFMIN ) THEN 219 ISCALE = 2 220 CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, 221 $ INFO ) 222 CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, 223 $ INFO ) 224 ENDIF 225* 226* Choose between QL and QR iteration 227* 228 IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN 229 LEND = LSV 230 L = LENDSV 231 ENDIF 232* 233 IF( LEND.GT.L ) THEN 234* 235* QL Iteration 236* 237* Look for small subdiagonal element. 238* 239 40 CONTINUE 240 IF( L.NE.LEND ) THEN 241 LENDM1 = LEND - 1 242 DO 50 M = L, LENDM1 243 TST = ABS( E( M ) )**2 244 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+ 245 $ SAFMIN )GOTO 60 246 50 CONTINUE 247 ENDIF 248* 249 M = LEND 250* 251 60 CONTINUE 252 IF( M.LT.LEND ) 253 $ E( M ) = ZERO 254 P = D( L ) 255 IF( M.EQ.L ) 256 $ GOTO 110 257* 258* If remaining matrix is 2-by-2, use SLAE2 or SLAEV2 259* to compute its eigensystem. 260* 261 IF( M.EQ.L+1 ) THEN 262 CALL SLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S ) 263 WORK( L ) = C 264 WORK( N-1+L ) = S 265 CALL CLASR( 'R', 'V', 'B', NR, 2, WORK( L ), WORK( N-1+L ), 266 $ Z( 1, L ), LDZ ) 267 D( L ) = RT1 268 D( L+1 ) = RT2 269 E( L ) = ZERO 270 L = L + 2 271 IF( L.LE.LEND ) 272 $ GOTO 40 273 GOTO 200 274 ENDIF 275* 276 IF( JTOT.EQ.NMAXIT ) 277 $ GOTO 200 278 JTOT = JTOT + 1 279* 280* Form shift. 281* 282 G = ( D( L+1 )-P ) / ( TWO*E( L ) ) 283 R = SLAPY2( G, ONE ) 284 G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) ) 285* 286 IF( ICOMPZ.EQ.0 ) THEN 287* Do not do a lookahead! 288 GOTO 90 289 ENDIF 290* 291 OLDEL = ABS( E( L ) ) 292 GP = G 293 RP = R 294 TST = ABS( E( L ) )**2 295 TST = TST / ( ( EPS2*ABS( D( L ) ) )*ABS( D( L+1 ) )+SAFMIN ) 296* 297 NLOOK = 1 298 IF( ( TST.GT.ONE ) .AND. ( NLOOK.LE.NMAXLOOK ) ) THEN 299 70 CONTINUE 300* 301* This is the lookahead loop, going until we have 302* convergence or too many steps have been taken. 303* 304 S = ONE 305 C = ONE 306 P = ZERO 307 MM1 = M - 1 308 DO 80 I = MM1, L, -1 309 F = S*E( I ) 310 B = C*E( I ) 311 CALL SLARTG( GP, F, C, S, RP ) 312 GP = D( I+1 ) - P 313 RP = ( D( I )-GP )*S + TWO*C*B 314 P = S*RP 315 IF( I.NE.L ) 316 $ GP = C*RP - B 317 80 CONTINUE 318 OLDGP = GP 319 OLDRP = RP 320* Find GP & RP for the next iteration 321 IF( ABS( C*OLDRP-B ).GT.SAFMIN ) THEN 322 GP = ( ( OLDGP+P )-( D( L )-P ) ) / ( TWO*( C*OLDRP-B ) ) 323 ELSE 324* 325* Goto put in by G. Henry to fix ALPHA problem 326* 327 GOTO 90 328* GP = ( ( OLDGP+P )-( D( L )-P ) ) / 329* $ ( TWO*( C*OLDRP-B )+SAFMIN ) 330 ENDIF 331 RP = SLAPY2( GP, ONE ) 332 GP = D( M ) - ( D( L )-P ) + 333 $ ( ( C*OLDRP-B ) / ( GP+SIGN( RP, GP ) ) ) 334 TST1 = TST 335 TST = ABS( C*OLDRP-B )**2 336 TST = TST / ( ( EPS2*ABS( D( L )-P ) )*ABS( OLDGP+P )+ 337 $ SAFMIN ) 338* Make sure that we are making progress 339 IF( ABS( C*OLDRP-B ).GT.0.9E0*OLDEL ) THEN 340 IF( ABS( C*OLDRP-B ).GT.OLDEL ) THEN 341 GP = G 342 RP = R 343 ENDIF 344 TST = HALF 345 ELSE 346 OLDEL = ABS( C*OLDRP-B ) 347 ENDIF 348 NLOOK = NLOOK + 1 349 IF( ( TST.GT.ONE ) .AND. ( NLOOK.LE.NMAXLOOK ) ) 350 $ GOTO 70 351 ENDIF 352* 353 IF( ( TST.LE.ONE ) .AND. ( TST.NE.HALF ) .AND. 354 $ ( ABS( P ).LT.EPS*ABS( D( L ) ) ) .AND. 355 $ ( ILAST.EQ.L ) .AND. ( ABS( E( L ) )**2.LE.10000.0E0* 356 $ ( ( EPS2*ABS( D( L ) ) )*ABS( D( L+1 ) )+SAFMIN ) ) ) THEN 357* 358* Skip the current step: the subdiagonal info is just noise. 359* 360 M = L 361 E( M ) = ZERO 362 P = D( L ) 363 JTOT = JTOT - 1 364 GOTO 110 365 ENDIF 366 G = GP 367 R = RP 368* 369* Lookahead over 370* 371 90 CONTINUE 372* 373 S = ONE 374 C = ONE 375 P = ZERO 376* 377* Inner loop 378* 379 MM1 = M - 1 380 DO 100 I = MM1, L, -1 381 F = S*E( I ) 382 B = C*E( I ) 383 CALL SLARTG( G, F, C, S, R ) 384 IF( I.NE.M-1 ) 385 $ E( I+1 ) = R 386 G = D( I+1 ) - P 387 R = ( D( I )-G )*S + TWO*C*B 388 P = S*R 389 D( I+1 ) = G + P 390 G = C*R - B 391* 392* If eigenvectors are desired, then save rotations. 393* 394 WORK( I ) = C 395 WORK( N-1+I ) = -S 396* 397 100 CONTINUE 398* 399* If eigenvectors are desired, then apply saved rotations. 400* 401 MM = M - L + 1 402 CALL CLASR( 'R', 'V', 'B', NR, MM, WORK( L ), WORK( N-1+L ), 403 $ Z( 1, L ), LDZ ) 404* 405 D( L ) = D( L ) - P 406 E( L ) = G 407 ILAST = L 408 GOTO 40 409* 410* Eigenvalue found. 411* 412 110 CONTINUE 413 D( L ) = P 414* 415 L = L + 1 416 IF( L.LE.LEND ) 417 $ GOTO 40 418 GOTO 200 419* 420 ELSE 421* 422* QR Iteration 423* 424* Look for small superdiagonal element. 425* 426 120 CONTINUE 427 IF( L.NE.LEND ) THEN 428 LENDP1 = LEND + 1 429 DO 130 M = L, LENDP1, -1 430 TST = ABS( E( M-1 ) )**2 431 IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+ 432 $ SAFMIN )GOTO 140 433 130 CONTINUE 434 ENDIF 435* 436 M = LEND 437* 438 140 CONTINUE 439 IF( M.GT.LEND ) 440 $ E( M-1 ) = ZERO 441 P = D( L ) 442 IF( M.EQ.L ) 443 $ GOTO 190 444* 445* If remaining matrix is 2-by-2, use SLAE2 or SLAEV2 446* to compute its eigensystem. 447* 448 IF( M.EQ.L-1 ) THEN 449 CALL SLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S ) 450 WORK( M ) = C 451 WORK( N-1+M ) = S 452 CALL CLASR( 'R', 'V', 'F', NR, 2, WORK( M ), WORK( N-1+M ), 453 $ Z( 1, L-1 ), LDZ ) 454 D( L-1 ) = RT1 455 D( L ) = RT2 456 E( L-1 ) = ZERO 457 L = L - 2 458 IF( L.GE.LEND ) 459 $ GOTO 120 460 GOTO 200 461 ENDIF 462* 463 IF( JTOT.EQ.NMAXIT ) 464 $ GOTO 200 465 JTOT = JTOT + 1 466* 467* Form shift. 468* 469 G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) ) 470 R = SLAPY2( G, ONE ) 471 G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) ) 472* 473 IF( ICOMPZ.EQ.0 ) THEN 474* Do not do a lookahead! 475 GOTO 170 476 ENDIF 477* 478 OLDEL = ABS( E( L-1 ) ) 479 GP = G 480 RP = R 481 TST = ABS( E( L-1 ) )**2 482 TST = TST / ( ( EPS2*ABS( D( L ) ) )*ABS( D( L-1 ) )+SAFMIN ) 483 NLOOK = 1 484 IF( ( TST.GT.ONE ) .AND. ( NLOOK.LE.NMAXLOOK ) ) THEN 485 150 CONTINUE 486* 487* This is the lookahead loop, going until we have 488* convergence or too many steps have been taken. 489* 490 S = ONE 491 C = ONE 492 P = ZERO 493* 494* Inner loop 495* 496 LM1 = L - 1 497 DO 160 I = M, LM1 498 F = S*E( I ) 499 B = C*E( I ) 500 CALL SLARTG( GP, F, C, S, RP ) 501 GP = D( I ) - P 502 RP = ( D( I+1 )-GP )*S + TWO*C*B 503 P = S*RP 504 IF( I.LT.LM1 ) 505 $ GP = C*RP - B 506 160 CONTINUE 507 OLDGP = GP 508 OLDRP = RP 509* Find GP & RP for the next iteration 510 IF( ABS( C*OLDRP-B ).GT.SAFMIN ) THEN 511 GP = ( ( OLDGP+P )-( D( L )-P ) ) / ( TWO*( C*OLDRP-B ) ) 512 ELSE 513* 514* Goto put in by G. Henry to fix ALPHA problem 515* 516 GOTO 170 517* GP = ( ( OLDGP+P )-( D( L )-P ) ) / 518* $ ( TWO*( C*OLDRP-B )+SAFMIN ) 519 ENDIF 520 RP = SLAPY2( GP, ONE ) 521 GP = D( M ) - ( D( L )-P ) + 522 $ ( ( C*OLDRP-B ) / ( GP+SIGN( RP, GP ) ) ) 523 TST1 = TST 524 TST = ABS( ( C*OLDRP-B ) )**2 525 TST = TST / ( ( EPS2*ABS( D( L )-P ) )*ABS( OLDGP+P )+ 526 $ SAFMIN ) 527* Make sure that we are making progress 528 IF( ABS( C*OLDRP-B ).GT.0.9E0*OLDEL ) THEN 529 IF( ABS( C*OLDRP-B ).GT.OLDEL ) THEN 530 GP = G 531 RP = R 532 ENDIF 533 TST = HALF 534 ELSE 535 OLDEL = ABS( C*OLDRP-B ) 536 ENDIF 537 NLOOK = NLOOK + 1 538 IF( ( TST.GT.ONE ) .AND. ( NLOOK.LE.NMAXLOOK ) ) 539 $ GOTO 150 540 ENDIF 541 IF( ( TST.LE.ONE ) .AND. ( TST.NE.HALF ) .AND. 542 $ ( ABS( P ).LT.EPS*ABS( D( L ) ) ) .AND. 543 $ ( ILAST.EQ.L ) .AND. ( ABS( E( L-1 ) )**2.LE.10000.0E0* 544 $ ( ( EPS2*ABS( D( L-1 ) ) )*ABS( D( L ) )+SAFMIN ) ) ) THEN 545* 546* Skip the current step: the subdiagonal info is just noise. 547* 548 M = L 549 E( M-1 ) = ZERO 550 P = D( L ) 551 JTOT = JTOT - 1 552 GOTO 190 553 ENDIF 554* 555 G = GP 556 R = RP 557* 558* Lookahead over 559* 560 170 CONTINUE 561* 562 S = ONE 563 C = ONE 564 P = ZERO 565 DO 180 I = M, LM1 566 F = S*E( I ) 567 B = C*E( I ) 568 CALL SLARTG( G, F, C, S, R ) 569 IF( I.NE.M ) 570 $ E( I-1 ) = R 571 G = D( I ) - P 572 R = ( D( I+1 )-G )*S + TWO*C*B 573 P = S*R 574 D( I ) = G + P 575 G = C*R - B 576* 577* If eigenvectors are desired, then save rotations. 578* 579 WORK( I ) = C 580 WORK( N-1+I ) = S 581* 582 180 CONTINUE 583* 584* If eigenvectors are desired, then apply saved rotations. 585* 586 MM = L - M + 1 587 CALL CLASR( 'R', 'V', 'F', NR, MM, WORK( M ), WORK( N-1+M ), 588 $ Z( 1, M ), LDZ ) 589* 590 D( L ) = D( L ) - P 591 E( LM1 ) = G 592 ILAST = L 593 GOTO 120 594* 595* Eigenvalue found. 596* 597 190 CONTINUE 598 D( L ) = P 599* 600 L = L - 1 601 IF( L.GE.LEND ) 602 $ GOTO 120 603 GOTO 200 604* 605 ENDIF 606* 607* Undo scaling if necessary 608* 609 200 CONTINUE 610 IF( ISCALE.EQ.1 ) THEN 611 CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, 612 $ D( LSV ), N, INFO ) 613 CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ), 614 $ N, INFO ) 615 ELSEIF( ISCALE.EQ.2 ) THEN 616 CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, 617 $ D( LSV ), N, INFO ) 618 CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ), 619 $ N, INFO ) 620 ENDIF 621* 622* Check for no convergence to an eigenvalue after a total 623* of N*MAXIT iterations. 624* 625 IF( JTOT.LT.NMAXIT ) 626 $ GOTO 10 627 DO 210 I = 1, N - 1 628 IF( E( I ).NE.ZERO ) 629 $ INFO = INFO + 1 630 210 CONTINUE 631 GOTO 250 632* 633* Order eigenvalues and eigenvectors. 634* 635 220 CONTINUE 636* 637* Use Selection Sort to minimize swaps of eigenvectors 638* 639 DO 240 II = 2, N 640 I = II - 1 641 K = I 642 P = D( I ) 643 DO 230 J = II, N 644 IF( D( J ).LT.P ) THEN 645 K = J 646 P = D( J ) 647 ENDIF 648 230 CONTINUE 649 IF( K.NE.I ) THEN 650 D( K ) = D( I ) 651 D( I ) = P 652 CALL CSWAP( NR, Z( 1, I ), 1, Z( 1, K ), 1 ) 653 ENDIF 654 240 CONTINUE 655* 656 250 CONTINUE 657* WRITE( *, FMT = * )'JTOT', JTOT 658 RETURN 659* 660* End of SSTEQR2 661* 662 END 663