SRC/sggev3.f

   1 *> \brief <b> SGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SGGEV3 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggev3.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggev3.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggev3.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SGGEV3( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR,
  22 *      $                   ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK,
  23 *      $                   INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          JOBVL, JOBVR
  27 *       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
  31 *      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
  32 *      $                   VR( LDVR, * ), WORK( * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> SGGEV3 computes for a pair of N-by-N real nonsymmetric matrices (A,B)
  42 *> the generalized eigenvalues, and optionally, the left and/or right
  43 *> generalized eigenvectors.
  44 *>
  45 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
  46 *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
  47 *> singular. It is usually represented as the pair (alpha,beta), as
  48 *> there is a reasonable interpretation for beta=0, and even for both
  49 *> being zero.
  50 *>
  51 *> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
  52 *> of (A,B) satisfies
  53 *>
  54 *>                  A * v(j) = lambda(j) * B * v(j).
  55 *>
  56 *> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
  57 *> of (A,B) satisfies
  58 *>
  59 *>                  u(j)**H * A  = lambda(j) * u(j)**H * B .
  60 *>
  61 *> where u(j)**H is the conjugate-transpose of u(j).
  62 *>
  63 *> \endverbatim
  64 *
  65 *  Arguments:
  66 *  ==========
  67 *
  68 *> \param[in] JOBVL
  69 *> \verbatim
  70 *>          JOBVL is CHARACTER*1
  71 *>          = 'N':  do not compute the left generalized eigenvectors;
  72 *>          = 'V':  compute the left generalized eigenvectors.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] JOBVR
  76 *> \verbatim
  77 *>          JOBVR is CHARACTER*1
  78 *>          = 'N':  do not compute the right generalized eigenvectors;
  79 *>          = 'V':  compute the right generalized eigenvectors.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] N
  83 *> \verbatim
  84 *>          N is INTEGER
  85 *>          The order of the matrices A, B, VL, and VR.  N >= 0.
  86 *> \endverbatim
  87 *>
  88 *> \param[in,out] A
  89 *> \verbatim
  90 *>          A is REAL array, dimension (LDA, N)
  91 *>          On entry, the matrix A in the pair (A,B).
  92 *>          On exit, A has been overwritten.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] LDA
  96 *> \verbatim
  97 *>          LDA is INTEGER
  98 *>          The leading dimension of A.  LDA >= max(1,N).
  99 *> \endverbatim
 100 *>
 101 *> \param[in,out] B
 102 *> \verbatim
 103 *>          B is REAL array, dimension (LDB, N)
 104 *>          On entry, the matrix B in the pair (A,B).
 105 *>          On exit, B has been overwritten.
 106 *> \endverbatim
 107 *>
 108 *> \param[in] LDB
 109 *> \verbatim
 110 *>          LDB is INTEGER
 111 *>          The leading dimension of B.  LDB >= max(1,N).
 112 *> \endverbatim
 113 *>
 114 *> \param[out] ALPHAR
 115 *> \verbatim
 116 *>          ALPHAR is REAL array, dimension (N)
 117 *> \endverbatim
 118 *>
 119 *> \param[out] ALPHAI
 120 *> \verbatim
 121 *>          ALPHAI is REAL array, dimension (N)
 122 *> \endverbatim
 123 *>
 124 *> \param[out] BETA
 125 *> \verbatim
 126 *>          BETA is REAL array, dimension (N)
 127 *>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
 128 *>          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
 129 *>          the j-th eigenvalue is real; if positive, then the j-th and
 130 *>          (j+1)-st eigenvalues are a complex conjugate pair, with
 131 *>          ALPHAI(j+1) negative.
 132 *>
 133 *>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
 134 *>          may easily over- or underflow, and BETA(j) may even be zero.
 135 *>          Thus, the user should avoid naively computing the ratio
 136 *>          alpha/beta.  However, ALPHAR and ALPHAI will be always less
 137 *>          than and usually comparable with norm(A) in magnitude, and
 138 *>          BETA always less than and usually comparable with norm(B).
 139 *> \endverbatim
 140 *>
 141 *> \param[out] VL
 142 *> \verbatim
 143 *>          VL is REAL array, dimension (LDVL,N)
 144 *>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
 145 *>          after another in the columns of VL, in the same order as
 146 *>          their eigenvalues. If the j-th eigenvalue is real, then
 147 *>          u(j) = VL(:,j), the j-th column of VL. If the j-th and
 148 *>          (j+1)-th eigenvalues form a complex conjugate pair, then
 149 *>          u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
 150 *>          Each eigenvector is scaled so the largest component has
 151 *>          abs(real part)+abs(imag. part)=1.
 152 *>          Not referenced if JOBVL = 'N'.
 153 *> \endverbatim
 154 *>
 155 *> \param[in] LDVL
 156 *> \verbatim
 157 *>          LDVL is INTEGER
 158 *>          The leading dimension of the matrix VL. LDVL >= 1, and
 159 *>          if JOBVL = 'V', LDVL >= N.
 160 *> \endverbatim
 161 *>
 162 *> \param[out] VR
 163 *> \verbatim
 164 *>          VR is REAL array, dimension (LDVR,N)
 165 *>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
 166 *>          after another in the columns of VR, in the same order as
 167 *>          their eigenvalues. If the j-th eigenvalue is real, then
 168 *>          v(j) = VR(:,j), the j-th column of VR. If the j-th and
 169 *>          (j+1)-th eigenvalues form a complex conjugate pair, then
 170 *>          v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
 171 *>          Each eigenvector is scaled so the largest component has
 172 *>          abs(real part)+abs(imag. part)=1.
 173 *>          Not referenced if JOBVR = 'N'.
 174 *> \endverbatim
 175 *>
 176 *> \param[in] LDVR
 177 *> \verbatim
 178 *>          LDVR is INTEGER
 179 *>          The leading dimension of the matrix VR. LDVR >= 1, and
 180 *>          if JOBVR = 'V', LDVR >= N.
 181 *> \endverbatim
 182 *>
 183 *> \param[out] WORK
 184 *> \verbatim
 185 *>          WORK is REAL array, dimension (MAX(1,LWORK))
 186 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 187 *> \endverbatim
 188 *>
 189 *> \param[in] LWORK
 190 *> \verbatim
 191 *>          LWORK is INTEGER
 192 *>
 193 *>          If LWORK = -1, then a workspace query is assumed; the routine
 194 *>          only calculates the optimal size of the WORK array, returns
 195 *>          this value as the first entry of the WORK array, and no error
 196 *>          message related to LWORK is issued by XERBLA.
 197 *> \endverbatim
 198 *>
 199 *> \param[out] INFO
 200 *> \verbatim
 201 *>          INFO is INTEGER
 202 *>          = 0:  successful exit
 203 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 204 *>          = 1,...,N:
 205 *>                The QZ iteration failed.  No eigenvectors have been
 206 *>                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
 207 *>                should be correct for j=INFO+1,...,N.
 208 *>          > N:  =N+1: other than QZ iteration failed in SHGEQZ.
 209 *>                =N+2: error return from STGEVC.
 210 *> \endverbatim
 211 *
 212 *  Authors:
 213 *  ========
 214 *
 215 *> \author Univ. of Tennessee
 216 *> \author Univ. of California Berkeley
 217 *> \author Univ. of Colorado Denver
 218 *> \author NAG Ltd.
 219 *
 220 *> \date January 2015
 221 *
 222 *> \ingroup realGEeigen
 223 *
 224 *  =====================================================================
 225       SUBROUTINE SGGEV3( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR,
 226      $                   ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK,
 227      $                   INFO )
 228 *
 229 *  -- LAPACK driver routine (version 3.6.0) --
 230 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 231 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 232 *     January 2015
 233 *
 234 *     .. Scalar Arguments ..
 235       CHARACTER          JOBVL, JOBVR
 236       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
 237 *     ..
 238 *     .. Array Arguments ..
 239       REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
 240      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
 241      $                   VR( LDVR, * ), WORK( * )
 242 *     ..
 243 *
 244 *  =====================================================================
 245 *
 246 *     .. Parameters ..
 247       REAL               ZERO, ONE
 248       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 249 *     ..
 250 *     .. Local Scalars ..
 251       LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
 252       CHARACTER          CHTEMP
 253       INTEGER            ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
 254      $                   IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, LWKOPT
 255       REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
 256      $                   SMLNUM, TEMP
 257 *     ..
 258 *     .. Local Arrays ..
 259       LOGICAL            LDUMMA( 1 )
 260 *     ..
 261 *     .. External Subroutines ..
 262       EXTERNAL           SGEQRF, SGGBAK, SGGBAL, SGGHD3, SHGEQZ, SLABAD,
 263      $                   SLACPY, SLASCL, SLASET, SORGQR, SORMQR, STGEVC,
 264      $                   XERBLA
 265 *     ..
 266 *     .. External Functions ..
 267       LOGICAL            LSAME
 268       REAL               SLAMCH, SLANGE
 269       EXTERNAL           LSAME, SLAMCH, SLANGE
 270 *     ..
 271 *     .. Intrinsic Functions ..
 272       INTRINSIC          ABS, MAX, SQRT
 273 *     ..
 274 *     .. Executable Statements ..
 275 *
 276 *     Decode the input arguments
 277 *
 278       IF( LSAME( JOBVL, 'N' ) ) THEN
 279          IJOBVL = 1
 280          ILVL = .FALSE.
 281       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
 282          IJOBVL = 2
 283          ILVL = .TRUE.
 284       ELSE
 285          IJOBVL = -1
 286          ILVL = .FALSE.
 287       END IF
 288 *
 289       IF( LSAME( JOBVR, 'N' ) ) THEN
 290          IJOBVR = 1
 291          ILVR = .FALSE.
 292       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
 293          IJOBVR = 2
 294          ILVR = .TRUE.
 295       ELSE
 296          IJOBVR = -1
 297          ILVR = .FALSE.
 298       END IF
 299       ILV = ILVL .OR. ILVR
 300 *
 301 *     Test the input arguments
 302 *
 303       INFO = 0
 304       LQUERY = ( LWORK.EQ.-1 )
 305       IF( IJOBVL.LE.0 ) THEN
 306          INFO = -1
 307       ELSE IF( IJOBVR.LE.0 ) THEN
 308          INFO = -2
 309       ELSE IF( N.LT.0 ) THEN
 310          INFO = -3
 311       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 312          INFO = -5
 313       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 314          INFO = -7
 315       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
 316          INFO = -12
 317       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
 318          INFO = -14
 319       ELSE IF( LWORK.LT.MAX( 1, 8*N ) .AND. .NOT.LQUERY ) THEN
 320          INFO = -16
 321       END IF
 322 *
 323 *     Compute workspace
 324 *
 325       IF( INFO.EQ.0 ) THEN
 326          CALL SGEQRF( N, N, B, LDB, WORK, WORK, -1, IERR )
 327          LWKOPT = MAX( 1, 8*N, 3*N+INT ( WORK( 1 ) ) )
 328          CALL SORMQR( 'L', 'T', N, N, N, B, LDB, WORK, A, LDA, WORK,
 329      $                -1, IERR )
 330          LWKOPT = MAX( LWKOPT, 3*N+INT ( WORK( 1 ) ) )
 331          CALL SGGHD3( JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, VL, LDVL,
 332      $                VR, LDVR, WORK, -1, IERR )
 333          LWKOPT = MAX( LWKOPT, 3*N+INT ( WORK( 1 ) ) )
 334          IF( ILVL ) THEN
 335             CALL SORGQR( N, N, N, VL, LDVL, WORK, WORK, -1, IERR )
 336             LWKOPT = MAX( LWKOPT, 3*N+INT ( WORK( 1 ) ) )
 337             CALL SHGEQZ( 'S', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
 338      $                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
 339      $                   WORK, -1, IERR )
 340             LWKOPT = MAX( LWKOPT, 2*N+INT ( WORK( 1 ) ) )
 341          ELSE
 342             CALL SHGEQZ( 'E', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
 343      $                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
 344      $                   WORK, -1, IERR )
 345             LWKOPT = MAX( LWKOPT, 2*N+INT ( WORK( 1 ) ) )
 346          END IF
 347          WORK( 1 ) = REAL( LWKOPT )
 348 *
 349       END IF
 350 *
 351       IF( INFO.NE.0 ) THEN
 352          CALL XERBLA( 'SGGEV3 ', -INFO )
 353          RETURN
 354       ELSE IF( LQUERY ) THEN
 355          RETURN
 356       END IF
 357 *
 358 *     Quick return if possible
 359 *
 360       IF( N.EQ.0 )
 361      $   RETURN
 362 *
 363 *     Get machine constants
 364 *
 365       EPS = SLAMCH( 'P' )
 366       SMLNUM = SLAMCH( 'S' )
 367       BIGNUM = ONE / SMLNUM
 368       CALL SLABAD( SMLNUM, BIGNUM )
 369       SMLNUM = SQRT( SMLNUM ) / EPS
 370       BIGNUM = ONE / SMLNUM
 371 *
 372 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 373 *
 374       ANRM = SLANGE( 'M', N, N, A, LDA, WORK )
 375       ILASCL = .FALSE.
 376       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 377          ANRMTO = SMLNUM
 378          ILASCL = .TRUE.
 379       ELSE IF( ANRM.GT.BIGNUM ) THEN
 380          ANRMTO = BIGNUM
 381          ILASCL = .TRUE.
 382       END IF
 383       IF( ILASCL )
 384      $   CALL SLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
 385 *
 386 *     Scale B if max element outside range [SMLNUM,BIGNUM]
 387 *
 388       BNRM = SLANGE( 'M', N, N, B, LDB, WORK )
 389       ILBSCL = .FALSE.
 390       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
 391          BNRMTO = SMLNUM
 392          ILBSCL = .TRUE.
 393       ELSE IF( BNRM.GT.BIGNUM ) THEN
 394          BNRMTO = BIGNUM
 395          ILBSCL = .TRUE.
 396       END IF
 397       IF( ILBSCL )
 398      $   CALL SLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
 399 *
 400 *     Permute the matrices A, B to isolate eigenvalues if possible
 401 *
 402       ILEFT = 1
 403       IRIGHT = N + 1
 404       IWRK = IRIGHT + N
 405       CALL SGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
 406      $             WORK( IRIGHT ), WORK( IWRK ), IERR )
 407 *
 408 *     Reduce B to triangular form (QR decomposition of B)
 409 *
 410       IROWS = IHI + 1 - ILO
 411       IF( ILV ) THEN
 412          ICOLS = N + 1 - ILO
 413       ELSE
 414          ICOLS = IROWS
 415       END IF
 416       ITAU = IWRK
 417       IWRK = ITAU + IROWS
 418       CALL SGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
 419      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
 420 *
 421 *     Apply the orthogonal transformation to matrix A
 422 *
 423       CALL SORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
 424      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
 425      $             LWORK+1-IWRK, IERR )
 426 *
 427 *     Initialize VL
 428 *
 429       IF( ILVL ) THEN
 430          CALL SLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
 431          IF( IROWS.GT.1 ) THEN
 432             CALL SLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
 433      $                   VL( ILO+1, ILO ), LDVL )
 434          END IF
 435          CALL SORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
 436      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
 437       END IF
 438 *
 439 *     Initialize VR
 440 *
 441       IF( ILVR )
 442      $   CALL SLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
 443 *
 444 *     Reduce to generalized Hessenberg form
 445 *
 446       IF( ILV ) THEN
 447 *
 448 *        Eigenvectors requested -- work on whole matrix.
 449 *
 450          CALL SGGHD3( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
 451      $                LDVL, VR, LDVR, WORK( IWRK ), LWORK+1-IWRK, IERR )
 452       ELSE
 453          CALL SGGHD3( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
 454      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR,
 455      $                WORK( IWRK ), LWORK+1-IWRK, IERR )
 456       END IF
 457 *
 458 *     Perform QZ algorithm (Compute eigenvalues, and optionally, the
 459 *     Schur forms and Schur vectors)
 460 *
 461       IWRK = ITAU
 462       IF( ILV ) THEN
 463          CHTEMP = 'S'
 464       ELSE
 465          CHTEMP = 'E'
 466       END IF
 467       CALL SHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
 468      $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
 469      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
 470       IF( IERR.NE.0 ) THEN
 471          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
 472             INFO = IERR
 473          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
 474             INFO = IERR - N
 475          ELSE
 476             INFO = N + 1
 477          END IF
 478          GO TO 110
 479       END IF
 480 *
 481 *     Compute Eigenvectors
 482 *
 483       IF( ILV ) THEN
 484          IF( ILVL ) THEN
 485             IF( ILVR ) THEN
 486                CHTEMP = 'B'
 487             ELSE
 488                CHTEMP = 'L'
 489             END IF
 490          ELSE
 491             CHTEMP = 'R'
 492          END IF
 493          CALL STGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
 494      $                VR, LDVR, N, IN, WORK( IWRK ), IERR )
 495          IF( IERR.NE.0 ) THEN
 496             INFO = N + 2
 497             GO TO 110
 498          END IF
 499 *
 500 *        Undo balancing on VL and VR and normalization
 501 *
 502          IF( ILVL ) THEN
 503             CALL SGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
 504      $                   WORK( IRIGHT ), N, VL, LDVL, IERR )
 505             DO 50 JC = 1, N
 506                IF( ALPHAI( JC ).LT.ZERO )
 507      $            GO TO 50
 508                TEMP = ZERO
 509                IF( ALPHAI( JC ).EQ.ZERO ) THEN
 510                   DO 10 JR = 1, N
 511                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
 512    10             CONTINUE
 513                ELSE
 514                   DO 20 JR = 1, N
 515                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
 516      $                      ABS( VL( JR, JC+1 ) ) )
 517    20             CONTINUE
 518                END IF
 519                IF( TEMP.LT.SMLNUM )
 520      $            GO TO 50
 521                TEMP = ONE / TEMP
 522                IF( ALPHAI( JC ).EQ.ZERO ) THEN
 523                   DO 30 JR = 1, N
 524                      VL( JR, JC ) = VL( JR, JC )*TEMP
 525    30             CONTINUE
 526                ELSE
 527                   DO 40 JR = 1, N
 528                      VL( JR, JC ) = VL( JR, JC )*TEMP
 529                      VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
 530    40             CONTINUE
 531                END IF
 532    50       CONTINUE
 533          END IF
 534          IF( ILVR ) THEN
 535             CALL SGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
 536      $                   WORK( IRIGHT ), N, VR, LDVR, IERR )
 537             DO 100 JC = 1, N
 538                IF( ALPHAI( JC ).LT.ZERO )
 539      $            GO TO 100
 540                TEMP = ZERO
 541                IF( ALPHAI( JC ).EQ.ZERO ) THEN
 542                   DO 60 JR = 1, N
 543                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
 544    60             CONTINUE
 545                ELSE
 546                   DO 70 JR = 1, N
 547                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
 548      $                      ABS( VR( JR, JC+1 ) ) )
 549    70             CONTINUE
 550                END IF
 551                IF( TEMP.LT.SMLNUM )
 552      $            GO TO 100
 553                TEMP = ONE / TEMP
 554                IF( ALPHAI( JC ).EQ.ZERO ) THEN
 555                   DO 80 JR = 1, N
 556                      VR( JR, JC ) = VR( JR, JC )*TEMP
 557    80             CONTINUE
 558                ELSE
 559                   DO 90 JR = 1, N
 560                      VR( JR, JC ) = VR( JR, JC )*TEMP
 561                      VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
 562    90             CONTINUE
 563                END IF
 564   100       CONTINUE
 565          END IF
 566 *
 567 *        End of eigenvector calculation
 568 *
 569       END IF
 570 *
 571 *     Undo scaling if necessary
 572 *
 573   110 CONTINUE
 574 *
 575       IF( ILASCL ) THEN
 576          CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
 577          CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
 578       END IF
 579 *
 580       IF( ILBSCL ) THEN
 581          CALL SLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
 582       END IF
 583 *
 584       WORK( 1 ) = REAL( LWKOPT )
 585       RETURN
 586 *
 587 *     End of SGGEV3
 588 *
 589       END