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