1*> \brief \b ZTREVC3 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZTREVC3 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrevc3.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrevc3.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrevc3.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 22* $ LDVR, MM, M, WORK, LWORK, RWORK, LRWORK, INFO) 23* 24* .. Scalar Arguments .. 25* CHARACTER HOWMNY, SIDE 26* INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N 27* .. 28* .. Array Arguments .. 29* LOGICAL SELECT( * ) 30* DOUBLE PRECISION RWORK( * ) 31* COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 32* $ WORK( * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> ZTREVC3 computes some or all of the right and/or left eigenvectors of 42*> a complex upper triangular matrix T. 43*> Matrices of this type are produced by the Schur factorization of 44*> a complex general matrix: A = Q*T*Q**H, as computed by ZHSEQR. 45*> 46*> The right eigenvector x and the left eigenvector y of T corresponding 47*> to an eigenvalue w are defined by: 48*> 49*> T*x = w*x, (y**H)*T = w*(y**H) 50*> 51*> where y**H denotes the conjugate transpose of the vector y. 52*> The eigenvalues are not input to this routine, but are read directly 53*> from the diagonal of T. 54*> 55*> This routine returns the matrices X and/or Y of right and left 56*> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an 57*> input matrix. If Q is the unitary factor that reduces a matrix A to 58*> Schur form T, then Q*X and Q*Y are the matrices of right and left 59*> eigenvectors of A. 60*> 61*> This uses a Level 3 BLAS version of the back transformation. 62*> \endverbatim 63* 64* Arguments: 65* ========== 66* 67*> \param[in] SIDE 68*> \verbatim 69*> SIDE is CHARACTER*1 70*> = 'R': compute right eigenvectors only; 71*> = 'L': compute left eigenvectors only; 72*> = 'B': compute both right and left eigenvectors. 73*> \endverbatim 74*> 75*> \param[in] HOWMNY 76*> \verbatim 77*> HOWMNY is CHARACTER*1 78*> = 'A': compute all right and/or left eigenvectors; 79*> = 'B': compute all right and/or left eigenvectors, 80*> backtransformed using the matrices supplied in 81*> VR and/or VL; 82*> = 'S': compute selected right and/or left eigenvectors, 83*> as indicated by the logical array SELECT. 84*> \endverbatim 85*> 86*> \param[in] SELECT 87*> \verbatim 88*> SELECT is LOGICAL array, dimension (N) 89*> If HOWMNY = 'S', SELECT specifies the eigenvectors to be 90*> computed. 91*> The eigenvector corresponding to the j-th eigenvalue is 92*> computed if SELECT(j) = .TRUE.. 93*> Not referenced if HOWMNY = 'A' or 'B'. 94*> \endverbatim 95*> 96*> \param[in] N 97*> \verbatim 98*> N is INTEGER 99*> The order of the matrix T. N >= 0. 100*> \endverbatim 101*> 102*> \param[in,out] T 103*> \verbatim 104*> T is COMPLEX*16 array, dimension (LDT,N) 105*> The upper triangular matrix T. T is modified, but restored 106*> on exit. 107*> \endverbatim 108*> 109*> \param[in] LDT 110*> \verbatim 111*> LDT is INTEGER 112*> The leading dimension of the array T. LDT >= max(1,N). 113*> \endverbatim 114*> 115*> \param[in,out] VL 116*> \verbatim 117*> VL is COMPLEX*16 array, dimension (LDVL,MM) 118*> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 119*> contain an N-by-N matrix Q (usually the unitary matrix Q of 120*> Schur vectors returned by ZHSEQR). 121*> On exit, if SIDE = 'L' or 'B', VL contains: 122*> if HOWMNY = 'A', the matrix Y of left eigenvectors of T; 123*> if HOWMNY = 'B', the matrix Q*Y; 124*> if HOWMNY = 'S', the left eigenvectors of T specified by 125*> SELECT, stored consecutively in the columns 126*> of VL, in the same order as their 127*> eigenvalues. 128*> Not referenced if SIDE = 'R'. 129*> \endverbatim 130*> 131*> \param[in] LDVL 132*> \verbatim 133*> LDVL is INTEGER 134*> The leading dimension of the array VL. 135*> LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N. 136*> \endverbatim 137*> 138*> \param[in,out] VR 139*> \verbatim 140*> VR is COMPLEX*16 array, dimension (LDVR,MM) 141*> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 142*> contain an N-by-N matrix Q (usually the unitary matrix Q of 143*> Schur vectors returned by ZHSEQR). 144*> On exit, if SIDE = 'R' or 'B', VR contains: 145*> if HOWMNY = 'A', the matrix X of right eigenvectors of T; 146*> if HOWMNY = 'B', the matrix Q*X; 147*> if HOWMNY = 'S', the right eigenvectors of T specified by 148*> SELECT, stored consecutively in the columns 149*> of VR, in the same order as their 150*> eigenvalues. 151*> Not referenced if SIDE = 'L'. 152*> \endverbatim 153*> 154*> \param[in] LDVR 155*> \verbatim 156*> LDVR is INTEGER 157*> The leading dimension of the array VR. 158*> LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N. 159*> \endverbatim 160*> 161*> \param[in] MM 162*> \verbatim 163*> MM is INTEGER 164*> The number of columns in the arrays VL and/or VR. MM >= M. 165*> \endverbatim 166*> 167*> \param[out] M 168*> \verbatim 169*> M is INTEGER 170*> The number of columns in the arrays VL and/or VR actually 171*> used to store the eigenvectors. 172*> If HOWMNY = 'A' or 'B', M is set to N. 173*> Each selected eigenvector occupies one column. 174*> \endverbatim 175*> 176*> \param[out] WORK 177*> \verbatim 178*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 179*> \endverbatim 180*> 181*> \param[in] LWORK 182*> \verbatim 183*> LWORK is INTEGER 184*> The dimension of array WORK. LWORK >= max(1,2*N). 185*> For optimum performance, LWORK >= N + 2*N*NB, where NB is 186*> the optimal blocksize. 187*> 188*> If LWORK = -1, then a workspace query is assumed; the routine 189*> only calculates the optimal size of the WORK array, returns 190*> this value as the first entry of the WORK array, and no error 191*> message related to LWORK is issued by XERBLA. 192*> \endverbatim 193*> 194*> \param[out] RWORK 195*> \verbatim 196*> RWORK is DOUBLE PRECISION array, dimension (LRWORK) 197*> \endverbatim 198*> 199*> \param[in] LRWORK 200*> \verbatim 201*> LRWORK is INTEGER 202*> The dimension of array RWORK. LRWORK >= max(1,N). 203*> 204*> If LRWORK = -1, then a workspace query is assumed; the routine 205*> only calculates the optimal size of the RWORK array, returns 206*> this value as the first entry of the RWORK array, and no error 207*> message related to LRWORK is issued by XERBLA. 208*> \endverbatim 209*> 210*> \param[out] INFO 211*> \verbatim 212*> INFO is INTEGER 213*> = 0: successful exit 214*> < 0: if INFO = -i, the i-th argument had an illegal value 215*> \endverbatim 216* 217* Authors: 218* ======== 219* 220*> \author Univ. of Tennessee 221*> \author Univ. of California Berkeley 222*> \author Univ. of Colorado Denver 223*> \author NAG Ltd. 224* 225*> \ingroup complex16OTHERcomputational 226* 227*> \par Further Details: 228* ===================== 229*> 230*> \verbatim 231*> 232*> The algorithm used in this program is basically backward (forward) 233*> substitution, with scaling to make the the code robust against 234*> possible overflow. 235*> 236*> Each eigenvector is normalized so that the element of largest 237*> magnitude has magnitude 1; here the magnitude of a complex number 238*> (x,y) is taken to be |x| + |y|. 239*> \endverbatim 240*> 241* ===================================================================== 242 SUBROUTINE ZTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 243 $ LDVR, MM, M, WORK, LWORK, RWORK, LRWORK, INFO) 244 IMPLICIT NONE 245* 246* -- LAPACK computational routine -- 247* -- LAPACK is a software package provided by Univ. of Tennessee, -- 248* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 249* 250* .. Scalar Arguments .. 251 CHARACTER HOWMNY, SIDE 252 INTEGER INFO, LDT, LDVL, LDVR, LWORK, LRWORK, M, MM, N 253* .. 254* .. Array Arguments .. 255 LOGICAL SELECT( * ) 256 DOUBLE PRECISION RWORK( * ) 257 COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 258 $ WORK( * ) 259* .. 260* 261* ===================================================================== 262* 263* .. Parameters .. 264 DOUBLE PRECISION ZERO, ONE 265 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 266 COMPLEX*16 CZERO, CONE 267 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 268 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 269 INTEGER NBMIN, NBMAX 270 PARAMETER ( NBMIN = 8, NBMAX = 128 ) 271* .. 272* .. Local Scalars .. 273 LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, RIGHTV, SOMEV 274 INTEGER I, II, IS, J, K, KI, IV, MAXWRK, NB 275 DOUBLE PRECISION OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL 276 COMPLEX*16 CDUM 277* .. 278* .. External Functions .. 279 LOGICAL LSAME 280 INTEGER ILAENV, IZAMAX 281 DOUBLE PRECISION DLAMCH, DZASUM 282 EXTERNAL LSAME, ILAENV, IZAMAX, DLAMCH, DZASUM 283* .. 284* .. External Subroutines .. 285 EXTERNAL XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLATRS, 286 $ ZGEMM, DLABAD, ZLASET, ZLACPY 287* .. 288* .. Intrinsic Functions .. 289 INTRINSIC ABS, DBLE, DCMPLX, CONJG, DIMAG, MAX 290* .. 291* .. Statement Functions .. 292 DOUBLE PRECISION CABS1 293* .. 294* .. Statement Function definitions .. 295 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 296* .. 297* .. Executable Statements .. 298* 299* Decode and test the input parameters 300* 301 BOTHV = LSAME( SIDE, 'B' ) 302 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV 303 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV 304* 305 ALLV = LSAME( HOWMNY, 'A' ) 306 OVER = LSAME( HOWMNY, 'B' ) 307 SOMEV = LSAME( HOWMNY, 'S' ) 308* 309* Set M to the number of columns required to store the selected 310* eigenvectors. 311* 312 IF( SOMEV ) THEN 313 M = 0 314 DO 10 J = 1, N 315 IF( SELECT( J ) ) 316 $ M = M + 1 317 10 CONTINUE 318 ELSE 319 M = N 320 END IF 321* 322 INFO = 0 323 NB = ILAENV( 1, 'ZTREVC', SIDE // HOWMNY, N, -1, -1, -1 ) 324 MAXWRK = N + 2*N*NB 325 WORK(1) = MAXWRK 326 RWORK(1) = N 327 LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 ) 328 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN 329 INFO = -1 330 ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN 331 INFO = -2 332 ELSE IF( N.LT.0 ) THEN 333 INFO = -4 334 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 335 INFO = -6 336 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN 337 INFO = -8 338 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN 339 INFO = -10 340 ELSE IF( MM.LT.M ) THEN 341 INFO = -11 342 ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN 343 INFO = -14 344 ELSE IF ( LRWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 345 INFO = -16 346 END IF 347 IF( INFO.NE.0 ) THEN 348 CALL XERBLA( 'ZTREVC3', -INFO ) 349 RETURN 350 ELSE IF( LQUERY ) THEN 351 RETURN 352 END IF 353* 354* Quick return if possible. 355* 356 IF( N.EQ.0 ) 357 $ RETURN 358* 359* Use blocked version of back-transformation if sufficient workspace. 360* Zero-out the workspace to avoid potential NaN propagation. 361* 362 IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN 363 NB = (LWORK - N) / (2*N) 364 NB = MIN( NB, NBMAX ) 365 CALL ZLASET( 'F', N, 1+2*NB, CZERO, CZERO, WORK, N ) 366 ELSE 367 NB = 1 368 END IF 369* 370* Set the constants to control overflow. 371* 372 UNFL = DLAMCH( 'Safe minimum' ) 373 OVFL = ONE / UNFL 374 CALL DLABAD( UNFL, OVFL ) 375 ULP = DLAMCH( 'Precision' ) 376 SMLNUM = UNFL*( N / ULP ) 377* 378* Store the diagonal elements of T in working array WORK. 379* 380 DO 20 I = 1, N 381 WORK( I ) = T( I, I ) 382 20 CONTINUE 383* 384* Compute 1-norm of each column of strictly upper triangular 385* part of T to control overflow in triangular solver. 386* 387 RWORK( 1 ) = ZERO 388 DO 30 J = 2, N 389 RWORK( J ) = DZASUM( J-1, T( 1, J ), 1 ) 390 30 CONTINUE 391* 392 IF( RIGHTV ) THEN 393* 394* ============================================================ 395* Compute right eigenvectors. 396* 397* IV is index of column in current block. 398* Non-blocked version always uses IV=NB=1; 399* blocked version starts with IV=NB, goes down to 1. 400* (Note the "0-th" column is used to store the original diagonal.) 401 IV = NB 402 IS = M 403 DO 80 KI = N, 1, -1 404 IF( SOMEV ) THEN 405 IF( .NOT.SELECT( KI ) ) 406 $ GO TO 80 407 END IF 408 SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM ) 409* 410* -------------------------------------------------------- 411* Complex right eigenvector 412* 413 WORK( KI + IV*N ) = CONE 414* 415* Form right-hand side. 416* 417 DO 40 K = 1, KI - 1 418 WORK( K + IV*N ) = -T( K, KI ) 419 40 CONTINUE 420* 421* Solve upper triangular system: 422* [ T(1:KI-1,1:KI-1) - T(KI,KI) ]*X = SCALE*WORK. 423* 424 DO 50 K = 1, KI - 1 425 T( K, K ) = T( K, K ) - T( KI, KI ) 426 IF( CABS1( T( K, K ) ).LT.SMIN ) 427 $ T( K, K ) = SMIN 428 50 CONTINUE 429* 430 IF( KI.GT.1 ) THEN 431 CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y', 432 $ KI-1, T, LDT, WORK( 1 + IV*N ), SCALE, 433 $ RWORK, INFO ) 434 WORK( KI + IV*N ) = SCALE 435 END IF 436* 437* Copy the vector x or Q*x to VR and normalize. 438* 439 IF( .NOT.OVER ) THEN 440* ------------------------------ 441* no back-transform: copy x to VR and normalize. 442 CALL ZCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 ) 443* 444 II = IZAMAX( KI, VR( 1, IS ), 1 ) 445 REMAX = ONE / CABS1( VR( II, IS ) ) 446 CALL ZDSCAL( KI, REMAX, VR( 1, IS ), 1 ) 447* 448 DO 60 K = KI + 1, N 449 VR( K, IS ) = CZERO 450 60 CONTINUE 451* 452 ELSE IF( NB.EQ.1 ) THEN 453* ------------------------------ 454* version 1: back-transform each vector with GEMV, Q*x. 455 IF( KI.GT.1 ) 456 $ CALL ZGEMV( 'N', N, KI-1, CONE, VR, LDVR, 457 $ WORK( 1 + IV*N ), 1, DCMPLX( SCALE ), 458 $ VR( 1, KI ), 1 ) 459* 460 II = IZAMAX( N, VR( 1, KI ), 1 ) 461 REMAX = ONE / CABS1( VR( II, KI ) ) 462 CALL ZDSCAL( N, REMAX, VR( 1, KI ), 1 ) 463* 464 ELSE 465* ------------------------------ 466* version 2: back-transform block of vectors with GEMM 467* zero out below vector 468 DO K = KI + 1, N 469 WORK( K + IV*N ) = CZERO 470 END DO 471* 472* Columns IV:NB of work are valid vectors. 473* When the number of vectors stored reaches NB, 474* or if this was last vector, do the GEMM 475 IF( (IV.EQ.1) .OR. (KI.EQ.1) ) THEN 476 CALL ZGEMM( 'N', 'N', N, NB-IV+1, KI+NB-IV, CONE, 477 $ VR, LDVR, 478 $ WORK( 1 + (IV)*N ), N, 479 $ CZERO, 480 $ WORK( 1 + (NB+IV)*N ), N ) 481* normalize vectors 482 DO K = IV, NB 483 II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 ) 484 REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) ) 485 CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 ) 486 END DO 487 CALL ZLACPY( 'F', N, NB-IV+1, 488 $ WORK( 1 + (NB+IV)*N ), N, 489 $ VR( 1, KI ), LDVR ) 490 IV = NB 491 ELSE 492 IV = IV - 1 493 END IF 494 END IF 495* 496* Restore the original diagonal elements of T. 497* 498 DO 70 K = 1, KI - 1 499 T( K, K ) = WORK( K ) 500 70 CONTINUE 501* 502 IS = IS - 1 503 80 CONTINUE 504 END IF 505* 506 IF( LEFTV ) THEN 507* 508* ============================================================ 509* Compute left eigenvectors. 510* 511* IV is index of column in current block. 512* Non-blocked version always uses IV=1; 513* blocked version starts with IV=1, goes up to NB. 514* (Note the "0-th" column is used to store the original diagonal.) 515 IV = 1 516 IS = 1 517 DO 130 KI = 1, N 518* 519 IF( SOMEV ) THEN 520 IF( .NOT.SELECT( KI ) ) 521 $ GO TO 130 522 END IF 523 SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM ) 524* 525* -------------------------------------------------------- 526* Complex left eigenvector 527* 528 WORK( KI + IV*N ) = CONE 529* 530* Form right-hand side. 531* 532 DO 90 K = KI + 1, N 533 WORK( K + IV*N ) = -CONJG( T( KI, K ) ) 534 90 CONTINUE 535* 536* Solve conjugate-transposed triangular system: 537* [ T(KI+1:N,KI+1:N) - T(KI,KI) ]**H * X = SCALE*WORK. 538* 539 DO 100 K = KI + 1, N 540 T( K, K ) = T( K, K ) - T( KI, KI ) 541 IF( CABS1( T( K, K ) ).LT.SMIN ) 542 $ T( K, K ) = SMIN 543 100 CONTINUE 544* 545 IF( KI.LT.N ) THEN 546 CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', 547 $ 'Y', N-KI, T( KI+1, KI+1 ), LDT, 548 $ WORK( KI+1 + IV*N ), SCALE, RWORK, INFO ) 549 WORK( KI + IV*N ) = SCALE 550 END IF 551* 552* Copy the vector x or Q*x to VL and normalize. 553* 554 IF( .NOT.OVER ) THEN 555* ------------------------------ 556* no back-transform: copy x to VL and normalize. 557 CALL ZCOPY( N-KI+1, WORK( KI + IV*N ), 1, VL(KI,IS), 1 ) 558* 559 II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1 560 REMAX = ONE / CABS1( VL( II, IS ) ) 561 CALL ZDSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 562* 563 DO 110 K = 1, KI - 1 564 VL( K, IS ) = CZERO 565 110 CONTINUE 566* 567 ELSE IF( NB.EQ.1 ) THEN 568* ------------------------------ 569* version 1: back-transform each vector with GEMV, Q*x. 570 IF( KI.LT.N ) 571 $ CALL ZGEMV( 'N', N, N-KI, CONE, VL( 1, KI+1 ), LDVL, 572 $ WORK( KI+1 + IV*N ), 1, DCMPLX( SCALE ), 573 $ VL( 1, KI ), 1 ) 574* 575 II = IZAMAX( N, VL( 1, KI ), 1 ) 576 REMAX = ONE / CABS1( VL( II, KI ) ) 577 CALL ZDSCAL( N, REMAX, VL( 1, KI ), 1 ) 578* 579 ELSE 580* ------------------------------ 581* version 2: back-transform block of vectors with GEMM 582* zero out above vector 583* could go from KI-NV+1 to KI-1 584 DO K = 1, KI - 1 585 WORK( K + IV*N ) = CZERO 586 END DO 587* 588* Columns 1:IV of work are valid vectors. 589* When the number of vectors stored reaches NB, 590* or if this was last vector, do the GEMM 591 IF( (IV.EQ.NB) .OR. (KI.EQ.N) ) THEN 592 CALL ZGEMM( 'N', 'N', N, IV, N-KI+IV, CONE, 593 $ VL( 1, KI-IV+1 ), LDVL, 594 $ WORK( KI-IV+1 + (1)*N ), N, 595 $ CZERO, 596 $ WORK( 1 + (NB+1)*N ), N ) 597* normalize vectors 598 DO K = 1, IV 599 II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 ) 600 REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) ) 601 CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 ) 602 END DO 603 CALL ZLACPY( 'F', N, IV, 604 $ WORK( 1 + (NB+1)*N ), N, 605 $ VL( 1, KI-IV+1 ), LDVL ) 606 IV = 1 607 ELSE 608 IV = IV + 1 609 END IF 610 END IF 611* 612* Restore the original diagonal elements of T. 613* 614 DO 120 K = KI + 1, N 615 T( K, K ) = WORK( K ) 616 120 CONTINUE 617* 618 IS = IS + 1 619 130 CONTINUE 620 END IF 621* 622 RETURN 623* 624* End of ZTREVC3 625* 626 END 627