1*> \brief \b DTGEVC 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DTGEVC + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgevc.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgevc.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgevc.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, 22* LDVL, VR, LDVR, MM, M, WORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER HOWMNY, SIDE 26* INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N 27* .. 28* .. Array Arguments .. 29* LOGICAL SELECT( * ) 30* DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ), 31* $ VR( LDVR, * ), WORK( * ) 32* .. 33* 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> DTGEVC computes some or all of the right and/or left eigenvectors of 42*> a pair of real matrices (S,P), where S is a quasi-triangular matrix 43*> and P is upper triangular. Matrix pairs of this type are produced by 44*> the generalized Schur factorization of a matrix pair (A,B): 45*> 46*> A = Q*S*Z**T, B = Q*P*Z**T 47*> 48*> as computed by DGGHRD + DHGEQZ. 49*> 50*> The right eigenvector x and the left eigenvector y of (S,P) 51*> corresponding to an eigenvalue w are defined by: 52*> 53*> S*x = w*P*x, (y**H)*S = w*(y**H)*P, 54*> 55*> where y**H denotes the conjugate tranpose of y. 56*> The eigenvalues are not input to this routine, but are computed 57*> directly from the diagonal blocks of S and P. 58*> 59*> This routine returns the matrices X and/or Y of right and left 60*> eigenvectors of (S,P), or the products Z*X and/or Q*Y, 61*> where Z and Q are input matrices. 62*> If Q and Z are the orthogonal factors from the generalized Schur 63*> factorization of a matrix pair (A,B), then Z*X and Q*Y 64*> are the matrices of right and left eigenvectors of (A,B). 65*> 66*> \endverbatim 67* 68* Arguments: 69* ========== 70* 71*> \param[in] SIDE 72*> \verbatim 73*> SIDE is CHARACTER*1 74*> = 'R': compute right eigenvectors only; 75*> = 'L': compute left eigenvectors only; 76*> = 'B': compute both right and left eigenvectors. 77*> \endverbatim 78*> 79*> \param[in] HOWMNY 80*> \verbatim 81*> HOWMNY is CHARACTER*1 82*> = 'A': compute all right and/or left eigenvectors; 83*> = 'B': compute all right and/or left eigenvectors, 84*> backtransformed by the matrices in VR and/or VL; 85*> = 'S': compute selected right and/or left eigenvectors, 86*> specified by the logical array SELECT. 87*> \endverbatim 88*> 89*> \param[in] SELECT 90*> \verbatim 91*> SELECT is LOGICAL array, dimension (N) 92*> If HOWMNY='S', SELECT specifies the eigenvectors to be 93*> computed. If w(j) is a real eigenvalue, the corresponding 94*> real eigenvector is computed if SELECT(j) is .TRUE.. 95*> If w(j) and w(j+1) are the real and imaginary parts of a 96*> complex eigenvalue, the corresponding complex eigenvector 97*> is computed if either SELECT(j) or SELECT(j+1) is .TRUE., 98*> and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is 99*> set to .FALSE.. 100*> Not referenced if HOWMNY = 'A' or 'B'. 101*> \endverbatim 102*> 103*> \param[in] N 104*> \verbatim 105*> N is INTEGER 106*> The order of the matrices S and P. N >= 0. 107*> \endverbatim 108*> 109*> \param[in] S 110*> \verbatim 111*> S is DOUBLE PRECISION array, dimension (LDS,N) 112*> The upper quasi-triangular matrix S from a generalized Schur 113*> factorization, as computed by DHGEQZ. 114*> \endverbatim 115*> 116*> \param[in] LDS 117*> \verbatim 118*> LDS is INTEGER 119*> The leading dimension of array S. LDS >= max(1,N). 120*> \endverbatim 121*> 122*> \param[in] P 123*> \verbatim 124*> P is DOUBLE PRECISION array, dimension (LDP,N) 125*> The upper triangular matrix P from a generalized Schur 126*> factorization, as computed by DHGEQZ. 127*> 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks 128*> of S must be in positive diagonal form. 129*> \endverbatim 130*> 131*> \param[in] LDP 132*> \verbatim 133*> LDP is INTEGER 134*> The leading dimension of array P. LDP >= max(1,N). 135*> \endverbatim 136*> 137*> \param[in,out] VL 138*> \verbatim 139*> VL is DOUBLE PRECISION array, dimension (LDVL,MM) 140*> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 141*> contain an N-by-N matrix Q (usually the orthogonal matrix Q 142*> of left Schur vectors returned by DHGEQZ). 143*> On exit, if SIDE = 'L' or 'B', VL contains: 144*> if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P); 145*> if HOWMNY = 'B', the matrix Q*Y; 146*> if HOWMNY = 'S', the left eigenvectors of (S,P) specified by 147*> SELECT, stored consecutively in the columns of 148*> VL, in the same order as their eigenvalues. 149*> 150*> A complex eigenvector corresponding to a complex eigenvalue 151*> is stored in two consecutive columns, the first holding the 152*> real part, and the second the imaginary part. 153*> 154*> Not referenced if SIDE = 'R'. 155*> \endverbatim 156*> 157*> \param[in] LDVL 158*> \verbatim 159*> LDVL is INTEGER 160*> The leading dimension of array VL. LDVL >= 1, and if 161*> SIDE = 'L' or 'B', LDVL >= N. 162*> \endverbatim 163*> 164*> \param[in,out] VR 165*> \verbatim 166*> VR is DOUBLE PRECISION array, dimension (LDVR,MM) 167*> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 168*> contain an N-by-N matrix Z (usually the orthogonal matrix Z 169*> of right Schur vectors returned by DHGEQZ). 170*> 171*> On exit, if SIDE = 'R' or 'B', VR contains: 172*> if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P); 173*> if HOWMNY = 'B' or 'b', the matrix Z*X; 174*> if HOWMNY = 'S' or 's', the right eigenvectors of (S,P) 175*> specified by SELECT, stored consecutively in the 176*> columns of VR, in the same order as their 177*> eigenvalues. 178*> 179*> A complex eigenvector corresponding to a complex eigenvalue 180*> is stored in two consecutive columns, the first holding the 181*> real part and the second the imaginary part. 182*> 183*> Not referenced if SIDE = 'L'. 184*> \endverbatim 185*> 186*> \param[in] LDVR 187*> \verbatim 188*> LDVR is INTEGER 189*> The leading dimension of the array VR. LDVR >= 1, and if 190*> SIDE = 'R' or 'B', LDVR >= N. 191*> \endverbatim 192*> 193*> \param[in] MM 194*> \verbatim 195*> MM is INTEGER 196*> The number of columns in the arrays VL and/or VR. MM >= M. 197*> \endverbatim 198*> 199*> \param[out] M 200*> \verbatim 201*> M is INTEGER 202*> The number of columns in the arrays VL and/or VR actually 203*> used to store the eigenvectors. If HOWMNY = 'A' or 'B', M 204*> is set to N. Each selected real eigenvector occupies one 205*> column and each selected complex eigenvector occupies two 206*> columns. 207*> \endverbatim 208*> 209*> \param[out] WORK 210*> \verbatim 211*> WORK is DOUBLE PRECISION array, dimension (6*N) 212*> \endverbatim 213*> 214*> \param[out] INFO 215*> \verbatim 216*> INFO is INTEGER 217*> = 0: successful exit. 218*> < 0: if INFO = -i, the i-th argument had an illegal value. 219*> > 0: the 2-by-2 block (INFO:INFO+1) does not have a complex 220*> eigenvalue. 221*> \endverbatim 222* 223* Authors: 224* ======== 225* 226*> \author Univ. of Tennessee 227*> \author Univ. of California Berkeley 228*> \author Univ. of Colorado Denver 229*> \author NAG Ltd. 230* 231*> \date December 2016 232* 233*> \ingroup doubleGEcomputational 234* 235*> \par Further Details: 236* ===================== 237*> 238*> \verbatim 239*> 240*> Allocation of workspace: 241*> ---------- -- --------- 242*> 243*> WORK( j ) = 1-norm of j-th column of A, above the diagonal 244*> WORK( N+j ) = 1-norm of j-th column of B, above the diagonal 245*> WORK( 2*N+1:3*N ) = real part of eigenvector 246*> WORK( 3*N+1:4*N ) = imaginary part of eigenvector 247*> WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector 248*> WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector 249*> 250*> Rowwise vs. columnwise solution methods: 251*> ------- -- ---------- -------- ------- 252*> 253*> Finding a generalized eigenvector consists basically of solving the 254*> singular triangular system 255*> 256*> (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left) 257*> 258*> Consider finding the i-th right eigenvector (assume all eigenvalues 259*> are real). The equation to be solved is: 260*> n i 261*> 0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1 262*> k=j k=j 263*> 264*> where C = (A - w B) (The components v(i+1:n) are 0.) 265*> 266*> The "rowwise" method is: 267*> 268*> (1) v(i) := 1 269*> for j = i-1,. . .,1: 270*> i 271*> (2) compute s = - sum C(j,k) v(k) and 272*> k=j+1 273*> 274*> (3) v(j) := s / C(j,j) 275*> 276*> Step 2 is sometimes called the "dot product" step, since it is an 277*> inner product between the j-th row and the portion of the eigenvector 278*> that has been computed so far. 279*> 280*> The "columnwise" method consists basically in doing the sums 281*> for all the rows in parallel. As each v(j) is computed, the 282*> contribution of v(j) times the j-th column of C is added to the 283*> partial sums. Since FORTRAN arrays are stored columnwise, this has 284*> the advantage that at each step, the elements of C that are accessed 285*> are adjacent to one another, whereas with the rowwise method, the 286*> elements accessed at a step are spaced LDS (and LDP) words apart. 287*> 288*> When finding left eigenvectors, the matrix in question is the 289*> transpose of the one in storage, so the rowwise method then 290*> actually accesses columns of A and B at each step, and so is the 291*> preferred method. 292*> \endverbatim 293*> 294* ===================================================================== 295 SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, 296 $ LDVL, VR, LDVR, MM, M, WORK, INFO ) 297* 298* -- LAPACK computational routine (version 3.7.0) -- 299* -- LAPACK is a software package provided by Univ. of Tennessee, -- 300* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 301* December 2016 302* 303* .. Scalar Arguments .. 304 CHARACTER HOWMNY, SIDE 305 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N 306* .. 307* .. Array Arguments .. 308 LOGICAL SELECT( * ) 309 DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ), 310 $ VR( LDVR, * ), WORK( * ) 311* .. 312* 313* 314* ===================================================================== 315* 316* .. Parameters .. 317 DOUBLE PRECISION ZERO, ONE, SAFETY 318 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, 319 $ SAFETY = 1.0D+2 ) 320* .. 321* .. Local Scalars .. 322 LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK, 323 $ ILBBAD, ILCOMP, ILCPLX, LSA, LSB 324 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE, 325 $ J, JA, JC, JE, JR, JW, NA, NW 326 DOUBLE PRECISION ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI, 327 $ BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A, 328 $ CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA, 329 $ CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE, 330 $ SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX, 331 $ XSCALE 332* .. 333* .. Local Arrays .. 334 DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ), 335 $ SUMP( 2, 2 ) 336* .. 337* .. External Functions .. 338 LOGICAL LSAME 339 DOUBLE PRECISION DLAMCH 340 EXTERNAL LSAME, DLAMCH 341* .. 342* .. External Subroutines .. 343 EXTERNAL DGEMV, DLABAD, DLACPY, DLAG2, DLALN2, XERBLA 344* .. 345* .. Intrinsic Functions .. 346 INTRINSIC ABS, MAX, MIN 347* .. 348* .. Executable Statements .. 349* 350* Decode and Test the input parameters 351* 352 IF( LSAME( HOWMNY, 'A' ) ) THEN 353 IHWMNY = 1 354 ILALL = .TRUE. 355 ILBACK = .FALSE. 356 ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN 357 IHWMNY = 2 358 ILALL = .FALSE. 359 ILBACK = .FALSE. 360 ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN 361 IHWMNY = 3 362 ILALL = .TRUE. 363 ILBACK = .TRUE. 364 ELSE 365 IHWMNY = -1 366 ILALL = .TRUE. 367 END IF 368* 369 IF( LSAME( SIDE, 'R' ) ) THEN 370 ISIDE = 1 371 COMPL = .FALSE. 372 COMPR = .TRUE. 373 ELSE IF( LSAME( SIDE, 'L' ) ) THEN 374 ISIDE = 2 375 COMPL = .TRUE. 376 COMPR = .FALSE. 377 ELSE IF( LSAME( SIDE, 'B' ) ) THEN 378 ISIDE = 3 379 COMPL = .TRUE. 380 COMPR = .TRUE. 381 ELSE 382 ISIDE = -1 383 END IF 384* 385 INFO = 0 386 IF( ISIDE.LT.0 ) THEN 387 INFO = -1 388 ELSE IF( IHWMNY.LT.0 ) THEN 389 INFO = -2 390 ELSE IF( N.LT.0 ) THEN 391 INFO = -4 392 ELSE IF( LDS.LT.MAX( 1, N ) ) THEN 393 INFO = -6 394 ELSE IF( LDP.LT.MAX( 1, N ) ) THEN 395 INFO = -8 396 END IF 397 IF( INFO.NE.0 ) THEN 398 CALL XERBLA( 'DTGEVC', -INFO ) 399 RETURN 400 END IF 401* 402* Count the number of eigenvectors to be computed 403* 404 IF( .NOT.ILALL ) THEN 405 IM = 0 406 ILCPLX = .FALSE. 407 DO 10 J = 1, N 408 IF( ILCPLX ) THEN 409 ILCPLX = .FALSE. 410 GO TO 10 411 END IF 412 IF( J.LT.N ) THEN 413 IF( S( J+1, J ).NE.ZERO ) 414 $ ILCPLX = .TRUE. 415 END IF 416 IF( ILCPLX ) THEN 417 IF( SELECT( J ) .OR. SELECT( J+1 ) ) 418 $ IM = IM + 2 419 ELSE 420 IF( SELECT( J ) ) 421 $ IM = IM + 1 422 END IF 423 10 CONTINUE 424 ELSE 425 IM = N 426 END IF 427* 428* Check 2-by-2 diagonal blocks of A, B 429* 430 ILABAD = .FALSE. 431 ILBBAD = .FALSE. 432 DO 20 J = 1, N - 1 433 IF( S( J+1, J ).NE.ZERO ) THEN 434 IF( P( J, J ).EQ.ZERO .OR. P( J+1, J+1 ).EQ.ZERO .OR. 435 $ P( J, J+1 ).NE.ZERO )ILBBAD = .TRUE. 436 IF( J.LT.N-1 ) THEN 437 IF( S( J+2, J+1 ).NE.ZERO ) 438 $ ILABAD = .TRUE. 439 END IF 440 END IF 441 20 CONTINUE 442* 443 IF( ILABAD ) THEN 444 INFO = -5 445 ELSE IF( ILBBAD ) THEN 446 INFO = -7 447 ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN 448 INFO = -10 449 ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN 450 INFO = -12 451 ELSE IF( MM.LT.IM ) THEN 452 INFO = -13 453 END IF 454 IF( INFO.NE.0 ) THEN 455 CALL XERBLA( 'DTGEVC', -INFO ) 456 RETURN 457 END IF 458* 459* Quick return if possible 460* 461 M = IM 462 IF( N.EQ.0 ) 463 $ RETURN 464* 465* Machine Constants 466* 467 SAFMIN = DLAMCH( 'Safe minimum' ) 468 BIG = ONE / SAFMIN 469 CALL DLABAD( SAFMIN, BIG ) 470 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 471 SMALL = SAFMIN*N / ULP 472 BIG = ONE / SMALL 473 BIGNUM = ONE / ( SAFMIN*N ) 474* 475* Compute the 1-norm of each column of the strictly upper triangular 476* part (i.e., excluding all elements belonging to the diagonal 477* blocks) of A and B to check for possible overflow in the 478* triangular solver. 479* 480 ANORM = ABS( S( 1, 1 ) ) 481 IF( N.GT.1 ) 482 $ ANORM = ANORM + ABS( S( 2, 1 ) ) 483 BNORM = ABS( P( 1, 1 ) ) 484 WORK( 1 ) = ZERO 485 WORK( N+1 ) = ZERO 486* 487 DO 50 J = 2, N 488 TEMP = ZERO 489 TEMP2 = ZERO 490 IF( S( J, J-1 ).EQ.ZERO ) THEN 491 IEND = J - 1 492 ELSE 493 IEND = J - 2 494 END IF 495 DO 30 I = 1, IEND 496 TEMP = TEMP + ABS( S( I, J ) ) 497 TEMP2 = TEMP2 + ABS( P( I, J ) ) 498 30 CONTINUE 499 WORK( J ) = TEMP 500 WORK( N+J ) = TEMP2 501 DO 40 I = IEND + 1, MIN( J+1, N ) 502 TEMP = TEMP + ABS( S( I, J ) ) 503 TEMP2 = TEMP2 + ABS( P( I, J ) ) 504 40 CONTINUE 505 ANORM = MAX( ANORM, TEMP ) 506 BNORM = MAX( BNORM, TEMP2 ) 507 50 CONTINUE 508* 509 ASCALE = ONE / MAX( ANORM, SAFMIN ) 510 BSCALE = ONE / MAX( BNORM, SAFMIN ) 511* 512* Left eigenvectors 513* 514 IF( COMPL ) THEN 515 IEIG = 0 516* 517* Main loop over eigenvalues 518* 519 ILCPLX = .FALSE. 520 DO 220 JE = 1, N 521* 522* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or 523* (b) this would be the second of a complex pair. 524* Check for complex eigenvalue, so as to be sure of which 525* entry(-ies) of SELECT to look at. 526* 527 IF( ILCPLX ) THEN 528 ILCPLX = .FALSE. 529 GO TO 220 530 END IF 531 NW = 1 532 IF( JE.LT.N ) THEN 533 IF( S( JE+1, JE ).NE.ZERO ) THEN 534 ILCPLX = .TRUE. 535 NW = 2 536 END IF 537 END IF 538 IF( ILALL ) THEN 539 ILCOMP = .TRUE. 540 ELSE IF( ILCPLX ) THEN 541 ILCOMP = SELECT( JE ) .OR. SELECT( JE+1 ) 542 ELSE 543 ILCOMP = SELECT( JE ) 544 END IF 545 IF( .NOT.ILCOMP ) 546 $ GO TO 220 547* 548* Decide if (a) singular pencil, (b) real eigenvalue, or 549* (c) complex eigenvalue. 550* 551 IF( .NOT.ILCPLX ) THEN 552 IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND. 553 $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN 554* 555* Singular matrix pencil -- return unit eigenvector 556* 557 IEIG = IEIG + 1 558 DO 60 JR = 1, N 559 VL( JR, IEIG ) = ZERO 560 60 CONTINUE 561 VL( IEIG, IEIG ) = ONE 562 GO TO 220 563 END IF 564 END IF 565* 566* Clear vector 567* 568 DO 70 JR = 1, NW*N 569 WORK( 2*N+JR ) = ZERO 570 70 CONTINUE 571* T 572* Compute coefficients in ( a A - b B ) y = 0 573* a is ACOEF 574* b is BCOEFR + i*BCOEFI 575* 576 IF( .NOT.ILCPLX ) THEN 577* 578* Real eigenvalue 579* 580 TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE, 581 $ ABS( P( JE, JE ) )*BSCALE, SAFMIN ) 582 SALFAR = ( TEMP*S( JE, JE ) )*ASCALE 583 SBETA = ( TEMP*P( JE, JE ) )*BSCALE 584 ACOEF = SBETA*ASCALE 585 BCOEFR = SALFAR*BSCALE 586 BCOEFI = ZERO 587* 588* Scale to avoid underflow 589* 590 SCALE = ONE 591 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL 592 LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT. 593 $ SMALL 594 IF( LSA ) 595 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 596 IF( LSB ) 597 $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )* 598 $ MIN( BNORM, BIG ) ) 599 IF( LSA .OR. LSB ) THEN 600 SCALE = MIN( SCALE, ONE / 601 $ ( SAFMIN*MAX( ONE, ABS( ACOEF ), 602 $ ABS( BCOEFR ) ) ) ) 603 IF( LSA ) THEN 604 ACOEF = ASCALE*( SCALE*SBETA ) 605 ELSE 606 ACOEF = SCALE*ACOEF 607 END IF 608 IF( LSB ) THEN 609 BCOEFR = BSCALE*( SCALE*SALFAR ) 610 ELSE 611 BCOEFR = SCALE*BCOEFR 612 END IF 613 END IF 614 ACOEFA = ABS( ACOEF ) 615 BCOEFA = ABS( BCOEFR ) 616* 617* First component is 1 618* 619 WORK( 2*N+JE ) = ONE 620 XMAX = ONE 621 ELSE 622* 623* Complex eigenvalue 624* 625 CALL DLAG2( S( JE, JE ), LDS, P( JE, JE ), LDP, 626 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, 627 $ BCOEFI ) 628 BCOEFI = -BCOEFI 629 IF( BCOEFI.EQ.ZERO ) THEN 630 INFO = JE 631 RETURN 632 END IF 633* 634* Scale to avoid over/underflow 635* 636 ACOEFA = ABS( ACOEF ) 637 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 638 SCALE = ONE 639 IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN ) 640 $ SCALE = ( SAFMIN / ULP ) / ACOEFA 641 IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN ) 642 $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA ) 643 IF( SAFMIN*ACOEFA.GT.ASCALE ) 644 $ SCALE = ASCALE / ( SAFMIN*ACOEFA ) 645 IF( SAFMIN*BCOEFA.GT.BSCALE ) 646 $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) ) 647 IF( SCALE.NE.ONE ) THEN 648 ACOEF = SCALE*ACOEF 649 ACOEFA = ABS( ACOEF ) 650 BCOEFR = SCALE*BCOEFR 651 BCOEFI = SCALE*BCOEFI 652 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 653 END IF 654* 655* Compute first two components of eigenvector 656* 657 TEMP = ACOEF*S( JE+1, JE ) 658 TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE ) 659 TEMP2I = -BCOEFI*P( JE, JE ) 660 IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN 661 WORK( 2*N+JE ) = ONE 662 WORK( 3*N+JE ) = ZERO 663 WORK( 2*N+JE+1 ) = -TEMP2R / TEMP 664 WORK( 3*N+JE+1 ) = -TEMP2I / TEMP 665 ELSE 666 WORK( 2*N+JE+1 ) = ONE 667 WORK( 3*N+JE+1 ) = ZERO 668 TEMP = ACOEF*S( JE, JE+1 ) 669 WORK( 2*N+JE ) = ( BCOEFR*P( JE+1, JE+1 )-ACOEF* 670 $ S( JE+1, JE+1 ) ) / TEMP 671 WORK( 3*N+JE ) = BCOEFI*P( JE+1, JE+1 ) / TEMP 672 END IF 673 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), 674 $ ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) ) 675 END IF 676* 677 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 678* 679* T 680* Triangular solve of (a A - b B) y = 0 681* 682* T 683* (rowwise in (a A - b B) , or columnwise in (a A - b B) ) 684* 685 IL2BY2 = .FALSE. 686* 687 DO 160 J = JE + NW, N 688 IF( IL2BY2 ) THEN 689 IL2BY2 = .FALSE. 690 GO TO 160 691 END IF 692* 693 NA = 1 694 BDIAG( 1 ) = P( J, J ) 695 IF( J.LT.N ) THEN 696 IF( S( J+1, J ).NE.ZERO ) THEN 697 IL2BY2 = .TRUE. 698 BDIAG( 2 ) = P( J+1, J+1 ) 699 NA = 2 700 END IF 701 END IF 702* 703* Check whether scaling is necessary for dot products 704* 705 XSCALE = ONE / MAX( ONE, XMAX ) 706 TEMP = MAX( WORK( J ), WORK( N+J ), 707 $ ACOEFA*WORK( J )+BCOEFA*WORK( N+J ) ) 708 IF( IL2BY2 ) 709 $ TEMP = MAX( TEMP, WORK( J+1 ), WORK( N+J+1 ), 710 $ ACOEFA*WORK( J+1 )+BCOEFA*WORK( N+J+1 ) ) 711 IF( TEMP.GT.BIGNUM*XSCALE ) THEN 712 DO 90 JW = 0, NW - 1 713 DO 80 JR = JE, J - 1 714 WORK( ( JW+2 )*N+JR ) = XSCALE* 715 $ WORK( ( JW+2 )*N+JR ) 716 80 CONTINUE 717 90 CONTINUE 718 XMAX = XMAX*XSCALE 719 END IF 720* 721* Compute dot products 722* 723* j-1 724* SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k) 725* k=je 726* 727* To reduce the op count, this is done as 728* 729* _ j-1 _ j-1 730* a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) ) 731* k=je k=je 732* 733* which may cause underflow problems if A or B are close 734* to underflow. (E.g., less than SMALL.) 735* 736* 737 DO 120 JW = 1, NW 738 DO 110 JA = 1, NA 739 SUMS( JA, JW ) = ZERO 740 SUMP( JA, JW ) = ZERO 741* 742 DO 100 JR = JE, J - 1 743 SUMS( JA, JW ) = SUMS( JA, JW ) + 744 $ S( JR, J+JA-1 )* 745 $ WORK( ( JW+1 )*N+JR ) 746 SUMP( JA, JW ) = SUMP( JA, JW ) + 747 $ P( JR, J+JA-1 )* 748 $ WORK( ( JW+1 )*N+JR ) 749 100 CONTINUE 750 110 CONTINUE 751 120 CONTINUE 752* 753 DO 130 JA = 1, NA 754 IF( ILCPLX ) THEN 755 SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) + 756 $ BCOEFR*SUMP( JA, 1 ) - 757 $ BCOEFI*SUMP( JA, 2 ) 758 SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) + 759 $ BCOEFR*SUMP( JA, 2 ) + 760 $ BCOEFI*SUMP( JA, 1 ) 761 ELSE 762 SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) + 763 $ BCOEFR*SUMP( JA, 1 ) 764 END IF 765 130 CONTINUE 766* 767* T 768* Solve ( a A - b B ) y = SUM(,) 769* with scaling and perturbation of the denominator 770* 771 CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS, 772 $ BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR, 773 $ BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP, 774 $ IINFO ) 775 IF( SCALE.LT.ONE ) THEN 776 DO 150 JW = 0, NW - 1 777 DO 140 JR = JE, J - 1 778 WORK( ( JW+2 )*N+JR ) = SCALE* 779 $ WORK( ( JW+2 )*N+JR ) 780 140 CONTINUE 781 150 CONTINUE 782 XMAX = SCALE*XMAX 783 END IF 784 XMAX = MAX( XMAX, TEMP ) 785 160 CONTINUE 786* 787* Copy eigenvector to VL, back transforming if 788* HOWMNY='B'. 789* 790 IEIG = IEIG + 1 791 IF( ILBACK ) THEN 792 DO 170 JW = 0, NW - 1 793 CALL DGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL, 794 $ WORK( ( JW+2 )*N+JE ), 1, ZERO, 795 $ WORK( ( JW+4 )*N+1 ), 1 ) 796 170 CONTINUE 797 CALL DLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ), 798 $ LDVL ) 799 IBEG = 1 800 ELSE 801 CALL DLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ), 802 $ LDVL ) 803 IBEG = JE 804 END IF 805* 806* Scale eigenvector 807* 808 XMAX = ZERO 809 IF( ILCPLX ) THEN 810 DO 180 J = IBEG, N 811 XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+ 812 $ ABS( VL( J, IEIG+1 ) ) ) 813 180 CONTINUE 814 ELSE 815 DO 190 J = IBEG, N 816 XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) ) 817 190 CONTINUE 818 END IF 819* 820 IF( XMAX.GT.SAFMIN ) THEN 821 XSCALE = ONE / XMAX 822* 823 DO 210 JW = 0, NW - 1 824 DO 200 JR = IBEG, N 825 VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW ) 826 200 CONTINUE 827 210 CONTINUE 828 END IF 829 IEIG = IEIG + NW - 1 830* 831 220 CONTINUE 832 END IF 833* 834* Right eigenvectors 835* 836 IF( COMPR ) THEN 837 IEIG = IM + 1 838* 839* Main loop over eigenvalues 840* 841 ILCPLX = .FALSE. 842 DO 500 JE = N, 1, -1 843* 844* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or 845* (b) this would be the second of a complex pair. 846* Check for complex eigenvalue, so as to be sure of which 847* entry(-ies) of SELECT to look at -- if complex, SELECT(JE) 848* or SELECT(JE-1). 849* If this is a complex pair, the 2-by-2 diagonal block 850* corresponding to the eigenvalue is in rows/columns JE-1:JE 851* 852 IF( ILCPLX ) THEN 853 ILCPLX = .FALSE. 854 GO TO 500 855 END IF 856 NW = 1 857 IF( JE.GT.1 ) THEN 858 IF( S( JE, JE-1 ).NE.ZERO ) THEN 859 ILCPLX = .TRUE. 860 NW = 2 861 END IF 862 END IF 863 IF( ILALL ) THEN 864 ILCOMP = .TRUE. 865 ELSE IF( ILCPLX ) THEN 866 ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 ) 867 ELSE 868 ILCOMP = SELECT( JE ) 869 END IF 870 IF( .NOT.ILCOMP ) 871 $ GO TO 500 872* 873* Decide if (a) singular pencil, (b) real eigenvalue, or 874* (c) complex eigenvalue. 875* 876 IF( .NOT.ILCPLX ) THEN 877 IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND. 878 $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN 879* 880* Singular matrix pencil -- unit eigenvector 881* 882 IEIG = IEIG - 1 883 DO 230 JR = 1, N 884 VR( JR, IEIG ) = ZERO 885 230 CONTINUE 886 VR( IEIG, IEIG ) = ONE 887 GO TO 500 888 END IF 889 END IF 890* 891* Clear vector 892* 893 DO 250 JW = 0, NW - 1 894 DO 240 JR = 1, N 895 WORK( ( JW+2 )*N+JR ) = ZERO 896 240 CONTINUE 897 250 CONTINUE 898* 899* Compute coefficients in ( a A - b B ) x = 0 900* a is ACOEF 901* b is BCOEFR + i*BCOEFI 902* 903 IF( .NOT.ILCPLX ) THEN 904* 905* Real eigenvalue 906* 907 TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE, 908 $ ABS( P( JE, JE ) )*BSCALE, SAFMIN ) 909 SALFAR = ( TEMP*S( JE, JE ) )*ASCALE 910 SBETA = ( TEMP*P( JE, JE ) )*BSCALE 911 ACOEF = SBETA*ASCALE 912 BCOEFR = SALFAR*BSCALE 913 BCOEFI = ZERO 914* 915* Scale to avoid underflow 916* 917 SCALE = ONE 918 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL 919 LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT. 920 $ SMALL 921 IF( LSA ) 922 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 923 IF( LSB ) 924 $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )* 925 $ MIN( BNORM, BIG ) ) 926 IF( LSA .OR. LSB ) THEN 927 SCALE = MIN( SCALE, ONE / 928 $ ( SAFMIN*MAX( ONE, ABS( ACOEF ), 929 $ ABS( BCOEFR ) ) ) ) 930 IF( LSA ) THEN 931 ACOEF = ASCALE*( SCALE*SBETA ) 932 ELSE 933 ACOEF = SCALE*ACOEF 934 END IF 935 IF( LSB ) THEN 936 BCOEFR = BSCALE*( SCALE*SALFAR ) 937 ELSE 938 BCOEFR = SCALE*BCOEFR 939 END IF 940 END IF 941 ACOEFA = ABS( ACOEF ) 942 BCOEFA = ABS( BCOEFR ) 943* 944* First component is 1 945* 946 WORK( 2*N+JE ) = ONE 947 XMAX = ONE 948* 949* Compute contribution from column JE of A and B to sum 950* (See "Further Details", above.) 951* 952 DO 260 JR = 1, JE - 1 953 WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) - 954 $ ACOEF*S( JR, JE ) 955 260 CONTINUE 956 ELSE 957* 958* Complex eigenvalue 959* 960 CALL DLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP, 961 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, 962 $ BCOEFI ) 963 IF( BCOEFI.EQ.ZERO ) THEN 964 INFO = JE - 1 965 RETURN 966 END IF 967* 968* Scale to avoid over/underflow 969* 970 ACOEFA = ABS( ACOEF ) 971 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 972 SCALE = ONE 973 IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN ) 974 $ SCALE = ( SAFMIN / ULP ) / ACOEFA 975 IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN ) 976 $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA ) 977 IF( SAFMIN*ACOEFA.GT.ASCALE ) 978 $ SCALE = ASCALE / ( SAFMIN*ACOEFA ) 979 IF( SAFMIN*BCOEFA.GT.BSCALE ) 980 $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) ) 981 IF( SCALE.NE.ONE ) THEN 982 ACOEF = SCALE*ACOEF 983 ACOEFA = ABS( ACOEF ) 984 BCOEFR = SCALE*BCOEFR 985 BCOEFI = SCALE*BCOEFI 986 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 987 END IF 988* 989* Compute first two components of eigenvector 990* and contribution to sums 991* 992 TEMP = ACOEF*S( JE, JE-1 ) 993 TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE ) 994 TEMP2I = -BCOEFI*P( JE, JE ) 995 IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN 996 WORK( 2*N+JE ) = ONE 997 WORK( 3*N+JE ) = ZERO 998 WORK( 2*N+JE-1 ) = -TEMP2R / TEMP 999 WORK( 3*N+JE-1 ) = -TEMP2I / TEMP 1000 ELSE 1001 WORK( 2*N+JE-1 ) = ONE 1002 WORK( 3*N+JE-1 ) = ZERO 1003 TEMP = ACOEF*S( JE-1, JE ) 1004 WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF* 1005 $ S( JE-1, JE-1 ) ) / TEMP 1006 WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP 1007 END IF 1008* 1009 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), 1010 $ ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) ) 1011* 1012* Compute contribution from columns JE and JE-1 1013* of A and B to the sums. 1014* 1015 CREALA = ACOEF*WORK( 2*N+JE-1 ) 1016 CIMAGA = ACOEF*WORK( 3*N+JE-1 ) 1017 CREALB = BCOEFR*WORK( 2*N+JE-1 ) - 1018 $ BCOEFI*WORK( 3*N+JE-1 ) 1019 CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) + 1020 $ BCOEFR*WORK( 3*N+JE-1 ) 1021 CRE2A = ACOEF*WORK( 2*N+JE ) 1022 CIM2A = ACOEF*WORK( 3*N+JE ) 1023 CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE ) 1024 CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE ) 1025 DO 270 JR = 1, JE - 2 1026 WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) + 1027 $ CREALB*P( JR, JE-1 ) - 1028 $ CRE2A*S( JR, JE ) + CRE2B*P( JR, JE ) 1029 WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) + 1030 $ CIMAGB*P( JR, JE-1 ) - 1031 $ CIM2A*S( JR, JE ) + CIM2B*P( JR, JE ) 1032 270 CONTINUE 1033 END IF 1034* 1035 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 1036* 1037* Columnwise triangular solve of (a A - b B) x = 0 1038* 1039 IL2BY2 = .FALSE. 1040 DO 370 J = JE - NW, 1, -1 1041* 1042* If a 2-by-2 block, is in position j-1:j, wait until 1043* next iteration to process it (when it will be j:j+1) 1044* 1045 IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN 1046 IF( S( J, J-1 ).NE.ZERO ) THEN 1047 IL2BY2 = .TRUE. 1048 GO TO 370 1049 END IF 1050 END IF 1051 BDIAG( 1 ) = P( J, J ) 1052 IF( IL2BY2 ) THEN 1053 NA = 2 1054 BDIAG( 2 ) = P( J+1, J+1 ) 1055 ELSE 1056 NA = 1 1057 END IF 1058* 1059* Compute x(j) (and x(j+1), if 2-by-2 block) 1060* 1061 CALL DLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ), 1062 $ LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ), 1063 $ N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP, 1064 $ IINFO ) 1065 IF( SCALE.LT.ONE ) THEN 1066* 1067 DO 290 JW = 0, NW - 1 1068 DO 280 JR = 1, JE 1069 WORK( ( JW+2 )*N+JR ) = SCALE* 1070 $ WORK( ( JW+2 )*N+JR ) 1071 280 CONTINUE 1072 290 CONTINUE 1073 END IF 1074 XMAX = MAX( SCALE*XMAX, TEMP ) 1075* 1076 DO 310 JW = 1, NW 1077 DO 300 JA = 1, NA 1078 WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW ) 1079 300 CONTINUE 1080 310 CONTINUE 1081* 1082* w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling 1083* 1084 IF( J.GT.1 ) THEN 1085* 1086* Check whether scaling is necessary for sum. 1087* 1088 XSCALE = ONE / MAX( ONE, XMAX ) 1089 TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J ) 1090 IF( IL2BY2 ) 1091 $ TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA* 1092 $ WORK( N+J+1 ) ) 1093 TEMP = MAX( TEMP, ACOEFA, BCOEFA ) 1094 IF( TEMP.GT.BIGNUM*XSCALE ) THEN 1095* 1096 DO 330 JW = 0, NW - 1 1097 DO 320 JR = 1, JE 1098 WORK( ( JW+2 )*N+JR ) = XSCALE* 1099 $ WORK( ( JW+2 )*N+JR ) 1100 320 CONTINUE 1101 330 CONTINUE 1102 XMAX = XMAX*XSCALE 1103 END IF 1104* 1105* Compute the contributions of the off-diagonals of 1106* column j (and j+1, if 2-by-2 block) of A and B to the 1107* sums. 1108* 1109* 1110 DO 360 JA = 1, NA 1111 IF( ILCPLX ) THEN 1112 CREALA = ACOEF*WORK( 2*N+J+JA-1 ) 1113 CIMAGA = ACOEF*WORK( 3*N+J+JA-1 ) 1114 CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) - 1115 $ BCOEFI*WORK( 3*N+J+JA-1 ) 1116 CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) + 1117 $ BCOEFR*WORK( 3*N+J+JA-1 ) 1118 DO 340 JR = 1, J - 1 1119 WORK( 2*N+JR ) = WORK( 2*N+JR ) - 1120 $ CREALA*S( JR, J+JA-1 ) + 1121 $ CREALB*P( JR, J+JA-1 ) 1122 WORK( 3*N+JR ) = WORK( 3*N+JR ) - 1123 $ CIMAGA*S( JR, J+JA-1 ) + 1124 $ CIMAGB*P( JR, J+JA-1 ) 1125 340 CONTINUE 1126 ELSE 1127 CREALA = ACOEF*WORK( 2*N+J+JA-1 ) 1128 CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) 1129 DO 350 JR = 1, J - 1 1130 WORK( 2*N+JR ) = WORK( 2*N+JR ) - 1131 $ CREALA*S( JR, J+JA-1 ) + 1132 $ CREALB*P( JR, J+JA-1 ) 1133 350 CONTINUE 1134 END IF 1135 360 CONTINUE 1136 END IF 1137* 1138 IL2BY2 = .FALSE. 1139 370 CONTINUE 1140* 1141* Copy eigenvector to VR, back transforming if 1142* HOWMNY='B'. 1143* 1144 IEIG = IEIG - NW 1145 IF( ILBACK ) THEN 1146* 1147 DO 410 JW = 0, NW - 1 1148 DO 380 JR = 1, N 1149 WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )* 1150 $ VR( JR, 1 ) 1151 380 CONTINUE 1152* 1153* A series of compiler directives to defeat 1154* vectorization for the next loop 1155* 1156* 1157 DO 400 JC = 2, JE 1158 DO 390 JR = 1, N 1159 WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) + 1160 $ WORK( ( JW+2 )*N+JC )*VR( JR, JC ) 1161 390 CONTINUE 1162 400 CONTINUE 1163 410 CONTINUE 1164* 1165 DO 430 JW = 0, NW - 1 1166 DO 420 JR = 1, N 1167 VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR ) 1168 420 CONTINUE 1169 430 CONTINUE 1170* 1171 IEND = N 1172 ELSE 1173 DO 450 JW = 0, NW - 1 1174 DO 440 JR = 1, N 1175 VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR ) 1176 440 CONTINUE 1177 450 CONTINUE 1178* 1179 IEND = JE 1180 END IF 1181* 1182* Scale eigenvector 1183* 1184 XMAX = ZERO 1185 IF( ILCPLX ) THEN 1186 DO 460 J = 1, IEND 1187 XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+ 1188 $ ABS( VR( J, IEIG+1 ) ) ) 1189 460 CONTINUE 1190 ELSE 1191 DO 470 J = 1, IEND 1192 XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) ) 1193 470 CONTINUE 1194 END IF 1195* 1196 IF( XMAX.GT.SAFMIN ) THEN 1197 XSCALE = ONE / XMAX 1198 DO 490 JW = 0, NW - 1 1199 DO 480 JR = 1, IEND 1200 VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW ) 1201 480 CONTINUE 1202 490 CONTINUE 1203 END IF 1204 500 CONTINUE 1205 END IF 1206* 1207 RETURN 1208* 1209* End of DTGEVC 1210* 1211 END 1212