LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
stgevc.f
Go to the documentation of this file.
00001 *> \brief \b STGEVC
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download STGEVC + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgevc.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgevc.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgevc.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE STGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
00022 *                          LDVL, VR, LDVR, MM, M, WORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          HOWMNY, SIDE
00026 *       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       LOGICAL            SELECT( * )
00030 *       REAL               P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
00031 *      $                   VR( LDVR, * ), WORK( * )
00032 *       ..
00033 *  
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> STGEVC computes some or all of the right and/or left eigenvectors of
00042 *> a pair of real matrices (S,P), where S is a quasi-triangular matrix
00043 *> and P is upper triangular.  Matrix pairs of this type are produced by
00044 *> the generalized Schur factorization of a matrix pair (A,B):
00045 *>
00046 *>    A = Q*S*Z**T,  B = Q*P*Z**T
00047 *>
00048 *> as computed by SGGHRD + SHGEQZ.
00049 *>
00050 *> The right eigenvector x and the left eigenvector y of (S,P)
00051 *> corresponding to an eigenvalue w are defined by:
00052 *> 
00053 *>    S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
00054 *> 
00055 *> where y**H denotes the conjugate tranpose of y.
00056 *> The eigenvalues are not input to this routine, but are computed
00057 *> directly from the diagonal blocks of S and P.
00058 *> 
00059 *> This routine returns the matrices X and/or Y of right and left
00060 *> eigenvectors of (S,P), or the products Z*X and/or Q*Y,
00061 *> where Z and Q are input matrices.
00062 *> If Q and Z are the orthogonal factors from the generalized Schur
00063 *> factorization of a matrix pair (A,B), then Z*X and Q*Y
00064 *> are the matrices of right and left eigenvectors of (A,B).
00065 *> 
00066 *> \endverbatim
00067 *
00068 *  Arguments:
00069 *  ==========
00070 *
00071 *> \param[in] SIDE
00072 *> \verbatim
00073 *>          SIDE is CHARACTER*1
00074 *>          = 'R': compute right eigenvectors only;
00075 *>          = 'L': compute left eigenvectors only;
00076 *>          = 'B': compute both right and left eigenvectors.
00077 *> \endverbatim
00078 *>
00079 *> \param[in] HOWMNY
00080 *> \verbatim
00081 *>          HOWMNY is CHARACTER*1
00082 *>          = 'A': compute all right and/or left eigenvectors;
00083 *>          = 'B': compute all right and/or left eigenvectors,
00084 *>                 backtransformed by the matrices in VR and/or VL;
00085 *>          = 'S': compute selected right and/or left eigenvectors,
00086 *>                 specified by the logical array SELECT.
00087 *> \endverbatim
00088 *>
00089 *> \param[in] SELECT
00090 *> \verbatim
00091 *>          SELECT is LOGICAL array, dimension (N)
00092 *>          If HOWMNY='S', SELECT specifies the eigenvectors to be
00093 *>          computed.  If w(j) is a real eigenvalue, the corresponding
00094 *>          real eigenvector is computed if SELECT(j) is .TRUE..
00095 *>          If w(j) and w(j+1) are the real and imaginary parts of a
00096 *>          complex eigenvalue, the corresponding complex eigenvector
00097 *>          is computed if either SELECT(j) or SELECT(j+1) is .TRUE.,
00098 *>          and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is
00099 *>          set to .FALSE..
00100 *>          Not referenced if HOWMNY = 'A' or 'B'.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] N
00104 *> \verbatim
00105 *>          N is INTEGER
00106 *>          The order of the matrices S and P.  N >= 0.
00107 *> \endverbatim
00108 *>
00109 *> \param[in] S
00110 *> \verbatim
00111 *>          S is REAL array, dimension (LDS,N)
00112 *>          The upper quasi-triangular matrix S from a generalized Schur
00113 *>          factorization, as computed by SHGEQZ.
00114 *> \endverbatim
00115 *>
00116 *> \param[in] LDS
00117 *> \verbatim
00118 *>          LDS is INTEGER
00119 *>          The leading dimension of array S.  LDS >= max(1,N).
00120 *> \endverbatim
00121 *>
00122 *> \param[in] P
00123 *> \verbatim
00124 *>          P is REAL array, dimension (LDP,N)
00125 *>          The upper triangular matrix P from a generalized Schur
00126 *>          factorization, as computed by SHGEQZ.
00127 *>          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
00128 *>          of S must be in positive diagonal form.
00129 *> \endverbatim
00130 *>
00131 *> \param[in] LDP
00132 *> \verbatim
00133 *>          LDP is INTEGER
00134 *>          The leading dimension of array P.  LDP >= max(1,N).
00135 *> \endverbatim
00136 *>
00137 *> \param[in,out] VL
00138 *> \verbatim
00139 *>          VL is REAL array, dimension (LDVL,MM)
00140 *>          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
00141 *>          contain an N-by-N matrix Q (usually the orthogonal matrix Q
00142 *>          of left Schur vectors returned by SHGEQZ).
00143 *>          On exit, if SIDE = 'L' or 'B', VL contains:
00144 *>          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
00145 *>          if HOWMNY = 'B', the matrix Q*Y;
00146 *>          if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
00147 *>                      SELECT, stored consecutively in the columns of
00148 *>                      VL, in the same order as their eigenvalues.
00149 *>
00150 *>          A complex eigenvector corresponding to a complex eigenvalue
00151 *>          is stored in two consecutive columns, the first holding the
00152 *>          real part, and the second the imaginary part.
00153 *>
00154 *>          Not referenced if SIDE = 'R'.
00155 *> \endverbatim
00156 *>
00157 *> \param[in] LDVL
00158 *> \verbatim
00159 *>          LDVL is INTEGER
00160 *>          The leading dimension of array VL.  LDVL >= 1, and if
00161 *>          SIDE = 'L' or 'B', LDVL >= N.
00162 *> \endverbatim
00163 *>
00164 *> \param[in,out] VR
00165 *> \verbatim
00166 *>          VR is REAL array, dimension (LDVR,MM)
00167 *>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
00168 *>          contain an N-by-N matrix Z (usually the orthogonal matrix Z
00169 *>          of right Schur vectors returned by SHGEQZ).
00170 *>
00171 *>          On exit, if SIDE = 'R' or 'B', VR contains:
00172 *>          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
00173 *>          if HOWMNY = 'B' or 'b', the matrix Z*X;
00174 *>          if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
00175 *>                      specified by SELECT, stored consecutively in the
00176 *>                      columns of VR, in the same order as their
00177 *>                      eigenvalues.
00178 *>
00179 *>          A complex eigenvector corresponding to a complex eigenvalue
00180 *>          is stored in two consecutive columns, the first holding the
00181 *>          real part and the second the imaginary part.
00182 *>          
00183 *>          Not referenced if SIDE = 'L'.
00184 *> \endverbatim
00185 *>
00186 *> \param[in] LDVR
00187 *> \verbatim
00188 *>          LDVR is INTEGER
00189 *>          The leading dimension of the array VR.  LDVR >= 1, and if
00190 *>          SIDE = 'R' or 'B', LDVR >= N.
00191 *> \endverbatim
00192 *>
00193 *> \param[in] MM
00194 *> \verbatim
00195 *>          MM is INTEGER
00196 *>          The number of columns in the arrays VL and/or VR. MM >= M.
00197 *> \endverbatim
00198 *>
00199 *> \param[out] M
00200 *> \verbatim
00201 *>          M is INTEGER
00202 *>          The number of columns in the arrays VL and/or VR actually
00203 *>          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
00204 *>          is set to N.  Each selected real eigenvector occupies one
00205 *>          column and each selected complex eigenvector occupies two
00206 *>          columns.
00207 *> \endverbatim
00208 *>
00209 *> \param[out] WORK
00210 *> \verbatim
00211 *>          WORK is REAL array, dimension (6*N)
00212 *> \endverbatim
00213 *>
00214 *> \param[out] INFO
00215 *> \verbatim
00216 *>          INFO is INTEGER
00217 *>          = 0:  successful exit.
00218 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00219 *>          > 0:  the 2-by-2 block (INFO:INFO+1) does not have a complex
00220 *>                eigenvalue.
00221 *> \endverbatim
00222 *
00223 *  Authors:
00224 *  ========
00225 *
00226 *> \author Univ. of Tennessee 
00227 *> \author Univ. of California Berkeley 
00228 *> \author Univ. of Colorado Denver 
00229 *> \author NAG Ltd. 
00230 *
00231 *> \date November 2011
00232 *
00233 *> \ingroup realGEcomputational
00234 *
00235 *> \par Further Details:
00236 *  =====================
00237 *>
00238 *> \verbatim
00239 *>
00240 *>  Allocation of workspace:
00241 *>  ---------- -- ---------
00242 *>
00243 *>     WORK( j ) = 1-norm of j-th column of A, above the diagonal
00244 *>     WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
00245 *>     WORK( 2*N+1:3*N ) = real part of eigenvector
00246 *>     WORK( 3*N+1:4*N ) = imaginary part of eigenvector
00247 *>     WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
00248 *>     WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
00249 *>
00250 *>  Rowwise vs. columnwise solution methods:
00251 *>  ------- --  ---------- -------- -------
00252 *>
00253 *>  Finding a generalized eigenvector consists basically of solving the
00254 *>  singular triangular system
00255 *>
00256 *>   (A - w B) x = 0     (for right) or:   (A - w B)**H y = 0  (for left)
00257 *>
00258 *>  Consider finding the i-th right eigenvector (assume all eigenvalues
00259 *>  are real). The equation to be solved is:
00260 *>       n                   i
00261 *>  0 = sum  C(j,k) v(k)  = sum  C(j,k) v(k)     for j = i,. . .,1
00262 *>      k=j                 k=j
00263 *>
00264 *>  where  C = (A - w B)  (The components v(i+1:n) are 0.)
00265 *>
00266 *>  The "rowwise" method is:
00267 *>
00268 *>  (1)  v(i) := 1
00269 *>  for j = i-1,. . .,1:
00270 *>                          i
00271 *>      (2) compute  s = - sum C(j,k) v(k)   and
00272 *>                        k=j+1
00273 *>
00274 *>      (3) v(j) := s / C(j,j)
00275 *>
00276 *>  Step 2 is sometimes called the "dot product" step, since it is an
00277 *>  inner product between the j-th row and the portion of the eigenvector
00278 *>  that has been computed so far.
00279 *>
00280 *>  The "columnwise" method consists basically in doing the sums
00281 *>  for all the rows in parallel.  As each v(j) is computed, the
00282 *>  contribution of v(j) times the j-th column of C is added to the
00283 *>  partial sums.  Since FORTRAN arrays are stored columnwise, this has
00284 *>  the advantage that at each step, the elements of C that are accessed
00285 *>  are adjacent to one another, whereas with the rowwise method, the
00286 *>  elements accessed at a step are spaced LDS (and LDP) words apart.
00287 *>
00288 *>  When finding left eigenvectors, the matrix in question is the
00289 *>  transpose of the one in storage, so the rowwise method then
00290 *>  actually accesses columns of A and B at each step, and so is the
00291 *>  preferred method.
00292 *> \endverbatim
00293 *>
00294 *  =====================================================================
00295       SUBROUTINE STGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
00296      $                   LDVL, VR, LDVR, MM, M, WORK, INFO )
00297 *
00298 *  -- LAPACK computational routine (version 3.4.0) --
00299 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00300 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00301 *     November 2011
00302 *
00303 *     .. Scalar Arguments ..
00304       CHARACTER          HOWMNY, SIDE
00305       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
00306 *     ..
00307 *     .. Array Arguments ..
00308       LOGICAL            SELECT( * )
00309       REAL               P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
00310      $                   VR( LDVR, * ), WORK( * )
00311 *     ..
00312 *
00313 *
00314 *  =====================================================================
00315 *
00316 *     .. Parameters ..
00317       REAL               ZERO, ONE, SAFETY
00318       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0,
00319      $                   SAFETY = 1.0E+2 )
00320 *     ..
00321 *     .. Local Scalars ..
00322       LOGICAL            COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
00323      $                   ILBBAD, ILCOMP, ILCPLX, LSA, LSB
00324       INTEGER            I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
00325      $                   J, JA, JC, JE, JR, JW, NA, NW
00326       REAL               ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
00327      $                   BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A,
00328      $                   CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA,
00329      $                   CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE,
00330      $                   SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX,
00331      $                   XSCALE
00332 *     ..
00333 *     .. Local Arrays ..
00334       REAL               BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
00335      $                   SUMP( 2, 2 )
00336 *     ..
00337 *     .. External Functions ..
00338       LOGICAL            LSAME
00339       REAL               SLAMCH
00340       EXTERNAL           LSAME, SLAMCH
00341 *     ..
00342 *     .. External Subroutines ..
00343       EXTERNAL           SGEMV, SLABAD, SLACPY, SLAG2, SLALN2, XERBLA
00344 *     ..
00345 *     .. Intrinsic Functions ..
00346       INTRINSIC          ABS, MAX, MIN
00347 *     ..
00348 *     .. Executable Statements ..
00349 *
00350 *     Decode and Test the input parameters
00351 *
00352       IF( LSAME( HOWMNY, 'A' ) ) THEN
00353          IHWMNY = 1
00354          ILALL = .TRUE.
00355          ILBACK = .FALSE.
00356       ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
00357          IHWMNY = 2
00358          ILALL = .FALSE.
00359          ILBACK = .FALSE.
00360       ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
00361          IHWMNY = 3
00362          ILALL = .TRUE.
00363          ILBACK = .TRUE.
00364       ELSE
00365          IHWMNY = -1
00366          ILALL = .TRUE.
00367       END IF
00368 *
00369       IF( LSAME( SIDE, 'R' ) ) THEN
00370          ISIDE = 1
00371          COMPL = .FALSE.
00372          COMPR = .TRUE.
00373       ELSE IF( LSAME( SIDE, 'L' ) ) THEN
00374          ISIDE = 2
00375          COMPL = .TRUE.
00376          COMPR = .FALSE.
00377       ELSE IF( LSAME( SIDE, 'B' ) ) THEN
00378          ISIDE = 3
00379          COMPL = .TRUE.
00380          COMPR = .TRUE.
00381       ELSE
00382          ISIDE = -1
00383       END IF
00384 *
00385       INFO = 0
00386       IF( ISIDE.LT.0 ) THEN
00387          INFO = -1
00388       ELSE IF( IHWMNY.LT.0 ) THEN
00389          INFO = -2
00390       ELSE IF( N.LT.0 ) THEN
00391          INFO = -4
00392       ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
00393          INFO = -6
00394       ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
00395          INFO = -8
00396       END IF
00397       IF( INFO.NE.0 ) THEN
00398          CALL XERBLA( 'STGEVC', -INFO )
00399          RETURN
00400       END IF
00401 *
00402 *     Count the number of eigenvectors to be computed
00403 *
00404       IF( .NOT.ILALL ) THEN
00405          IM = 0
00406          ILCPLX = .FALSE.
00407          DO 10 J = 1, N
00408             IF( ILCPLX ) THEN
00409                ILCPLX = .FALSE.
00410                GO TO 10
00411             END IF
00412             IF( J.LT.N ) THEN
00413                IF( S( J+1, J ).NE.ZERO )
00414      $            ILCPLX = .TRUE.
00415             END IF
00416             IF( ILCPLX ) THEN
00417                IF( SELECT( J ) .OR. SELECT( J+1 ) )
00418      $            IM = IM + 2
00419             ELSE
00420                IF( SELECT( J ) )
00421      $            IM = IM + 1
00422             END IF
00423    10    CONTINUE
00424       ELSE
00425          IM = N
00426       END IF
00427 *
00428 *     Check 2-by-2 diagonal blocks of A, B
00429 *
00430       ILABAD = .FALSE.
00431       ILBBAD = .FALSE.
00432       DO 20 J = 1, N - 1
00433          IF( S( J+1, J ).NE.ZERO ) THEN
00434             IF( P( J, J ).EQ.ZERO .OR. P( J+1, J+1 ).EQ.ZERO .OR.
00435      $          P( J, J+1 ).NE.ZERO )ILBBAD = .TRUE.
00436             IF( J.LT.N-1 ) THEN
00437                IF( S( J+2, J+1 ).NE.ZERO )
00438      $            ILABAD = .TRUE.
00439             END IF
00440          END IF
00441    20 CONTINUE
00442 *
00443       IF( ILABAD ) THEN
00444          INFO = -5
00445       ELSE IF( ILBBAD ) THEN
00446          INFO = -7
00447       ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
00448          INFO = -10
00449       ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
00450          INFO = -12
00451       ELSE IF( MM.LT.IM ) THEN
00452          INFO = -13
00453       END IF
00454       IF( INFO.NE.0 ) THEN
00455          CALL XERBLA( 'STGEVC', -INFO )
00456          RETURN
00457       END IF
00458 *
00459 *     Quick return if possible
00460 *
00461       M = IM
00462       IF( N.EQ.0 )
00463      $   RETURN
00464 *
00465 *     Machine Constants
00466 *
00467       SAFMIN = SLAMCH( 'Safe minimum' )
00468       BIG = ONE / SAFMIN
00469       CALL SLABAD( SAFMIN, BIG )
00470       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
00471       SMALL = SAFMIN*N / ULP
00472       BIG = ONE / SMALL
00473       BIGNUM = ONE / ( SAFMIN*N )
00474 *
00475 *     Compute the 1-norm of each column of the strictly upper triangular
00476 *     part (i.e., excluding all elements belonging to the diagonal
00477 *     blocks) of A and B to check for possible overflow in the
00478 *     triangular solver.
00479 *
00480       ANORM = ABS( S( 1, 1 ) )
00481       IF( N.GT.1 )
00482      $   ANORM = ANORM + ABS( S( 2, 1 ) )
00483       BNORM = ABS( P( 1, 1 ) )
00484       WORK( 1 ) = ZERO
00485       WORK( N+1 ) = ZERO
00486 *
00487       DO 50 J = 2, N
00488          TEMP = ZERO
00489          TEMP2 = ZERO
00490          IF( S( J, J-1 ).EQ.ZERO ) THEN
00491             IEND = J - 1
00492          ELSE
00493             IEND = J - 2
00494          END IF
00495          DO 30 I = 1, IEND
00496             TEMP = TEMP + ABS( S( I, J ) )
00497             TEMP2 = TEMP2 + ABS( P( I, J ) )
00498    30    CONTINUE
00499          WORK( J ) = TEMP
00500          WORK( N+J ) = TEMP2
00501          DO 40 I = IEND + 1, MIN( J+1, N )
00502             TEMP = TEMP + ABS( S( I, J ) )
00503             TEMP2 = TEMP2 + ABS( P( I, J ) )
00504    40    CONTINUE
00505          ANORM = MAX( ANORM, TEMP )
00506          BNORM = MAX( BNORM, TEMP2 )
00507    50 CONTINUE
00508 *
00509       ASCALE = ONE / MAX( ANORM, SAFMIN )
00510       BSCALE = ONE / MAX( BNORM, SAFMIN )
00511 *
00512 *     Left eigenvectors
00513 *
00514       IF( COMPL ) THEN
00515          IEIG = 0
00516 *
00517 *        Main loop over eigenvalues
00518 *
00519          ILCPLX = .FALSE.
00520          DO 220 JE = 1, N
00521 *
00522 *           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
00523 *           (b) this would be the second of a complex pair.
00524 *           Check for complex eigenvalue, so as to be sure of which
00525 *           entry(-ies) of SELECT to look at.
00526 *
00527             IF( ILCPLX ) THEN
00528                ILCPLX = .FALSE.
00529                GO TO 220
00530             END IF
00531             NW = 1
00532             IF( JE.LT.N ) THEN
00533                IF( S( JE+1, JE ).NE.ZERO ) THEN
00534                   ILCPLX = .TRUE.
00535                   NW = 2
00536                END IF
00537             END IF
00538             IF( ILALL ) THEN
00539                ILCOMP = .TRUE.
00540             ELSE IF( ILCPLX ) THEN
00541                ILCOMP = SELECT( JE ) .OR. SELECT( JE+1 )
00542             ELSE
00543                ILCOMP = SELECT( JE )
00544             END IF
00545             IF( .NOT.ILCOMP )
00546      $         GO TO 220
00547 *
00548 *           Decide if (a) singular pencil, (b) real eigenvalue, or
00549 *           (c) complex eigenvalue.
00550 *
00551             IF( .NOT.ILCPLX ) THEN
00552                IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
00553      $             ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
00554 *
00555 *                 Singular matrix pencil -- return unit eigenvector
00556 *
00557                   IEIG = IEIG + 1
00558                   DO 60 JR = 1, N
00559                      VL( JR, IEIG ) = ZERO
00560    60             CONTINUE
00561                   VL( IEIG, IEIG ) = ONE
00562                   GO TO 220
00563                END IF
00564             END IF
00565 *
00566 *           Clear vector
00567 *
00568             DO 70 JR = 1, NW*N
00569                WORK( 2*N+JR ) = ZERO
00570    70       CONTINUE
00571 *                                                 T
00572 *           Compute coefficients in  ( a A - b B )  y = 0
00573 *              a  is  ACOEF
00574 *              b  is  BCOEFR + i*BCOEFI
00575 *
00576             IF( .NOT.ILCPLX ) THEN
00577 *
00578 *              Real eigenvalue
00579 *
00580                TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
00581      $                ABS( P( JE, JE ) )*BSCALE, SAFMIN )
00582                SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
00583                SBETA = ( TEMP*P( JE, JE ) )*BSCALE
00584                ACOEF = SBETA*ASCALE
00585                BCOEFR = SALFAR*BSCALE
00586                BCOEFI = ZERO
00587 *
00588 *              Scale to avoid underflow
00589 *
00590                SCALE = ONE
00591                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
00592                LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
00593      $               SMALL
00594                IF( LSA )
00595      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
00596                IF( LSB )
00597      $            SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
00598      $                    MIN( BNORM, BIG ) )
00599                IF( LSA .OR. LSB ) THEN
00600                   SCALE = MIN( SCALE, ONE /
00601      $                    ( SAFMIN*MAX( ONE, ABS( ACOEF ),
00602      $                    ABS( BCOEFR ) ) ) )
00603                   IF( LSA ) THEN
00604                      ACOEF = ASCALE*( SCALE*SBETA )
00605                   ELSE
00606                      ACOEF = SCALE*ACOEF
00607                   END IF
00608                   IF( LSB ) THEN
00609                      BCOEFR = BSCALE*( SCALE*SALFAR )
00610                   ELSE
00611                      BCOEFR = SCALE*BCOEFR
00612                   END IF
00613                END IF
00614                ACOEFA = ABS( ACOEF )
00615                BCOEFA = ABS( BCOEFR )
00616 *
00617 *              First component is 1
00618 *
00619                WORK( 2*N+JE ) = ONE
00620                XMAX = ONE
00621             ELSE
00622 *
00623 *              Complex eigenvalue
00624 *
00625                CALL SLAG2( S( JE, JE ), LDS, P( JE, JE ), LDP,
00626      $                     SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
00627      $                     BCOEFI )
00628                BCOEFI = -BCOEFI
00629                IF( BCOEFI.EQ.ZERO ) THEN
00630                   INFO = JE
00631                   RETURN
00632                END IF
00633 *
00634 *              Scale to avoid over/underflow
00635 *
00636                ACOEFA = ABS( ACOEF )
00637                BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
00638                SCALE = ONE
00639                IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
00640      $            SCALE = ( SAFMIN / ULP ) / ACOEFA
00641                IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
00642      $            SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
00643                IF( SAFMIN*ACOEFA.GT.ASCALE )
00644      $            SCALE = ASCALE / ( SAFMIN*ACOEFA )
00645                IF( SAFMIN*BCOEFA.GT.BSCALE )
00646      $            SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
00647                IF( SCALE.NE.ONE ) THEN
00648                   ACOEF = SCALE*ACOEF
00649                   ACOEFA = ABS( ACOEF )
00650                   BCOEFR = SCALE*BCOEFR
00651                   BCOEFI = SCALE*BCOEFI
00652                   BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
00653                END IF
00654 *
00655 *              Compute first two components of eigenvector
00656 *
00657                TEMP = ACOEF*S( JE+1, JE )
00658                TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
00659                TEMP2I = -BCOEFI*P( JE, JE )
00660                IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
00661                   WORK( 2*N+JE ) = ONE
00662                   WORK( 3*N+JE ) = ZERO
00663                   WORK( 2*N+JE+1 ) = -TEMP2R / TEMP
00664                   WORK( 3*N+JE+1 ) = -TEMP2I / TEMP
00665                ELSE
00666                   WORK( 2*N+JE+1 ) = ONE
00667                   WORK( 3*N+JE+1 ) = ZERO
00668                   TEMP = ACOEF*S( JE, JE+1 )
00669                   WORK( 2*N+JE ) = ( BCOEFR*P( JE+1, JE+1 )-ACOEF*
00670      $                             S( JE+1, JE+1 ) ) / TEMP
00671                   WORK( 3*N+JE ) = BCOEFI*P( JE+1, JE+1 ) / TEMP
00672                END IF
00673                XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
00674      $                ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) )
00675             END IF
00676 *
00677             DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
00678 *
00679 *                                           T
00680 *           Triangular solve of  (a A - b B)  y = 0
00681 *
00682 *                                   T
00683 *           (rowwise in  (a A - b B) , or columnwise in (a A - b B) )
00684 *
00685             IL2BY2 = .FALSE.
00686 *
00687             DO 160 J = JE + NW, N
00688                IF( IL2BY2 ) THEN
00689                   IL2BY2 = .FALSE.
00690                   GO TO 160
00691                END IF
00692 *
00693                NA = 1
00694                BDIAG( 1 ) = P( J, J )
00695                IF( J.LT.N ) THEN
00696                   IF( S( J+1, J ).NE.ZERO ) THEN
00697                      IL2BY2 = .TRUE.
00698                      BDIAG( 2 ) = P( J+1, J+1 )
00699                      NA = 2
00700                   END IF
00701                END IF
00702 *
00703 *              Check whether scaling is necessary for dot products
00704 *
00705                XSCALE = ONE / MAX( ONE, XMAX )
00706                TEMP = MAX( WORK( J ), WORK( N+J ),
00707      $                ACOEFA*WORK( J )+BCOEFA*WORK( N+J ) )
00708                IF( IL2BY2 )
00709      $            TEMP = MAX( TEMP, WORK( J+1 ), WORK( N+J+1 ),
00710      $                   ACOEFA*WORK( J+1 )+BCOEFA*WORK( N+J+1 ) )
00711                IF( TEMP.GT.BIGNUM*XSCALE ) THEN
00712                   DO 90 JW = 0, NW - 1
00713                      DO 80 JR = JE, J - 1
00714                         WORK( ( JW+2 )*N+JR ) = XSCALE*
00715      $                     WORK( ( JW+2 )*N+JR )
00716    80                CONTINUE
00717    90             CONTINUE
00718                   XMAX = XMAX*XSCALE
00719                END IF
00720 *
00721 *              Compute dot products
00722 *
00723 *                    j-1
00724 *              SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
00725 *                    k=je
00726 *
00727 *              To reduce the op count, this is done as
00728 *
00729 *              _        j-1                  _        j-1
00730 *              a*conjg( sum  S(k,j)*x(k) ) - b*conjg( sum  P(k,j)*x(k) )
00731 *                       k=je                          k=je
00732 *
00733 *              which may cause underflow problems if A or B are close
00734 *              to underflow.  (E.g., less than SMALL.)
00735 *
00736 *
00737                DO 120 JW = 1, NW
00738                   DO 110 JA = 1, NA
00739                      SUMS( JA, JW ) = ZERO
00740                      SUMP( JA, JW ) = ZERO
00741 *
00742                      DO 100 JR = JE, J - 1
00743                         SUMS( JA, JW ) = SUMS( JA, JW ) +
00744      $                                   S( JR, J+JA-1 )*
00745      $                                   WORK( ( JW+1 )*N+JR )
00746                         SUMP( JA, JW ) = SUMP( JA, JW ) +
00747      $                                   P( JR, J+JA-1 )*
00748      $                                   WORK( ( JW+1 )*N+JR )
00749   100                CONTINUE
00750   110             CONTINUE
00751   120          CONTINUE
00752 *
00753                DO 130 JA = 1, NA
00754                   IF( ILCPLX ) THEN
00755                      SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
00756      $                              BCOEFR*SUMP( JA, 1 ) -
00757      $                              BCOEFI*SUMP( JA, 2 )
00758                      SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) +
00759      $                              BCOEFR*SUMP( JA, 2 ) +
00760      $                              BCOEFI*SUMP( JA, 1 )
00761                   ELSE
00762                      SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
00763      $                              BCOEFR*SUMP( JA, 1 )
00764                   END IF
00765   130          CONTINUE
00766 *
00767 *                                  T
00768 *              Solve  ( a A - b B )  y = SUM(,)
00769 *              with scaling and perturbation of the denominator
00770 *
00771                CALL SLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS,
00772      $                      BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR,
00773      $                      BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP,
00774      $                      IINFO )
00775                IF( SCALE.LT.ONE ) THEN
00776                   DO 150 JW = 0, NW - 1
00777                      DO 140 JR = JE, J - 1
00778                         WORK( ( JW+2 )*N+JR ) = SCALE*
00779      $                     WORK( ( JW+2 )*N+JR )
00780   140                CONTINUE
00781   150             CONTINUE
00782                   XMAX = SCALE*XMAX
00783                END IF
00784                XMAX = MAX( XMAX, TEMP )
00785   160       CONTINUE
00786 *
00787 *           Copy eigenvector to VL, back transforming if
00788 *           HOWMNY='B'.
00789 *
00790             IEIG = IEIG + 1
00791             IF( ILBACK ) THEN
00792                DO 170 JW = 0, NW - 1
00793                   CALL SGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL,
00794      $                        WORK( ( JW+2 )*N+JE ), 1, ZERO,
00795      $                        WORK( ( JW+4 )*N+1 ), 1 )
00796   170          CONTINUE
00797                CALL SLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ),
00798      $                      LDVL )
00799                IBEG = 1
00800             ELSE
00801                CALL SLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ),
00802      $                      LDVL )
00803                IBEG = JE
00804             END IF
00805 *
00806 *           Scale eigenvector
00807 *
00808             XMAX = ZERO
00809             IF( ILCPLX ) THEN
00810                DO 180 J = IBEG, N
00811                   XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+
00812      $                   ABS( VL( J, IEIG+1 ) ) )
00813   180          CONTINUE
00814             ELSE
00815                DO 190 J = IBEG, N
00816                   XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) )
00817   190          CONTINUE
00818             END IF
00819 *
00820             IF( XMAX.GT.SAFMIN ) THEN
00821                XSCALE = ONE / XMAX
00822 *
00823                DO 210 JW = 0, NW - 1
00824                   DO 200 JR = IBEG, N
00825                      VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW )
00826   200             CONTINUE
00827   210          CONTINUE
00828             END IF
00829             IEIG = IEIG + NW - 1
00830 *
00831   220    CONTINUE
00832       END IF
00833 *
00834 *     Right eigenvectors
00835 *
00836       IF( COMPR ) THEN
00837          IEIG = IM + 1
00838 *
00839 *        Main loop over eigenvalues
00840 *
00841          ILCPLX = .FALSE.
00842          DO 500 JE = N, 1, -1
00843 *
00844 *           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
00845 *           (b) this would be the second of a complex pair.
00846 *           Check for complex eigenvalue, so as to be sure of which
00847 *           entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
00848 *           or SELECT(JE-1).
00849 *           If this is a complex pair, the 2-by-2 diagonal block
00850 *           corresponding to the eigenvalue is in rows/columns JE-1:JE
00851 *
00852             IF( ILCPLX ) THEN
00853                ILCPLX = .FALSE.
00854                GO TO 500
00855             END IF
00856             NW = 1
00857             IF( JE.GT.1 ) THEN
00858                IF( S( JE, JE-1 ).NE.ZERO ) THEN
00859                   ILCPLX = .TRUE.
00860                   NW = 2
00861                END IF
00862             END IF
00863             IF( ILALL ) THEN
00864                ILCOMP = .TRUE.
00865             ELSE IF( ILCPLX ) THEN
00866                ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 )
00867             ELSE
00868                ILCOMP = SELECT( JE )
00869             END IF
00870             IF( .NOT.ILCOMP )
00871      $         GO TO 500
00872 *
00873 *           Decide if (a) singular pencil, (b) real eigenvalue, or
00874 *           (c) complex eigenvalue.
00875 *
00876             IF( .NOT.ILCPLX ) THEN
00877                IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
00878      $             ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
00879 *
00880 *                 Singular matrix pencil -- unit eigenvector
00881 *
00882                   IEIG = IEIG - 1
00883                   DO 230 JR = 1, N
00884                      VR( JR, IEIG ) = ZERO
00885   230             CONTINUE
00886                   VR( IEIG, IEIG ) = ONE
00887                   GO TO 500
00888                END IF
00889             END IF
00890 *
00891 *           Clear vector
00892 *
00893             DO 250 JW = 0, NW - 1
00894                DO 240 JR = 1, N
00895                   WORK( ( JW+2 )*N+JR ) = ZERO
00896   240          CONTINUE
00897   250       CONTINUE
00898 *
00899 *           Compute coefficients in  ( a A - b B ) x = 0
00900 *              a  is  ACOEF
00901 *              b  is  BCOEFR + i*BCOEFI
00902 *
00903             IF( .NOT.ILCPLX ) THEN
00904 *
00905 *              Real eigenvalue
00906 *
00907                TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
00908      $                ABS( P( JE, JE ) )*BSCALE, SAFMIN )
00909                SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
00910                SBETA = ( TEMP*P( JE, JE ) )*BSCALE
00911                ACOEF = SBETA*ASCALE
00912                BCOEFR = SALFAR*BSCALE
00913                BCOEFI = ZERO
00914 *
00915 *              Scale to avoid underflow
00916 *
00917                SCALE = ONE
00918                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
00919                LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
00920      $               SMALL
00921                IF( LSA )
00922      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
00923                IF( LSB )
00924      $            SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
00925      $                    MIN( BNORM, BIG ) )
00926                IF( LSA .OR. LSB ) THEN
00927                   SCALE = MIN( SCALE, ONE /
00928      $                    ( SAFMIN*MAX( ONE, ABS( ACOEF ),
00929      $                    ABS( BCOEFR ) ) ) )
00930                   IF( LSA ) THEN
00931                      ACOEF = ASCALE*( SCALE*SBETA )
00932                   ELSE
00933                      ACOEF = SCALE*ACOEF
00934                   END IF
00935                   IF( LSB ) THEN
00936                      BCOEFR = BSCALE*( SCALE*SALFAR )
00937                   ELSE
00938                      BCOEFR = SCALE*BCOEFR
00939                   END IF
00940                END IF
00941                ACOEFA = ABS( ACOEF )
00942                BCOEFA = ABS( BCOEFR )
00943 *
00944 *              First component is 1
00945 *
00946                WORK( 2*N+JE ) = ONE
00947                XMAX = ONE
00948 *
00949 *              Compute contribution from column JE of A and B to sum
00950 *              (See "Further Details", above.)
00951 *
00952                DO 260 JR = 1, JE - 1
00953                   WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) -
00954      $                             ACOEF*S( JR, JE )
00955   260          CONTINUE
00956             ELSE
00957 *
00958 *              Complex eigenvalue
00959 *
00960                CALL SLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP,
00961      $                     SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
00962      $                     BCOEFI )
00963                IF( BCOEFI.EQ.ZERO ) THEN
00964                   INFO = JE - 1
00965                   RETURN
00966                END IF
00967 *
00968 *              Scale to avoid over/underflow
00969 *
00970                ACOEFA = ABS( ACOEF )
00971                BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
00972                SCALE = ONE
00973                IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
00974      $            SCALE = ( SAFMIN / ULP ) / ACOEFA
00975                IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
00976      $            SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
00977                IF( SAFMIN*ACOEFA.GT.ASCALE )
00978      $            SCALE = ASCALE / ( SAFMIN*ACOEFA )
00979                IF( SAFMIN*BCOEFA.GT.BSCALE )
00980      $            SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
00981                IF( SCALE.NE.ONE ) THEN
00982                   ACOEF = SCALE*ACOEF
00983                   ACOEFA = ABS( ACOEF )
00984                   BCOEFR = SCALE*BCOEFR
00985                   BCOEFI = SCALE*BCOEFI
00986                   BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
00987                END IF
00988 *
00989 *              Compute first two components of eigenvector
00990 *              and contribution to sums
00991 *
00992                TEMP = ACOEF*S( JE, JE-1 )
00993                TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
00994                TEMP2I = -BCOEFI*P( JE, JE )
00995                IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
00996                   WORK( 2*N+JE ) = ONE
00997                   WORK( 3*N+JE ) = ZERO
00998                   WORK( 2*N+JE-1 ) = -TEMP2R / TEMP
00999                   WORK( 3*N+JE-1 ) = -TEMP2I / TEMP
01000                ELSE
01001                   WORK( 2*N+JE-1 ) = ONE
01002                   WORK( 3*N+JE-1 ) = ZERO
01003                   TEMP = ACOEF*S( JE-1, JE )
01004                   WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF*
01005      $                             S( JE-1, JE-1 ) ) / TEMP
01006                   WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP
01007                END IF
01008 *
01009                XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
01010      $                ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) )
01011 *
01012 *              Compute contribution from columns JE and JE-1
01013 *              of A and B to the sums.
01014 *
01015                CREALA = ACOEF*WORK( 2*N+JE-1 )
01016                CIMAGA = ACOEF*WORK( 3*N+JE-1 )
01017                CREALB = BCOEFR*WORK( 2*N+JE-1 ) -
01018      $                  BCOEFI*WORK( 3*N+JE-1 )
01019                CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) +
01020      $                  BCOEFR*WORK( 3*N+JE-1 )
01021                CRE2A = ACOEF*WORK( 2*N+JE )
01022                CIM2A = ACOEF*WORK( 3*N+JE )
01023                CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE )
01024                CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE )
01025                DO 270 JR = 1, JE - 2
01026                   WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) +
01027      $                             CREALB*P( JR, JE-1 ) -
01028      $                             CRE2A*S( JR, JE ) + CRE2B*P( JR, JE )
01029                   WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) +
01030      $                             CIMAGB*P( JR, JE-1 ) -
01031      $                             CIM2A*S( JR, JE ) + CIM2B*P( JR, JE )
01032   270          CONTINUE
01033             END IF
01034 *
01035             DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
01036 *
01037 *           Columnwise triangular solve of  (a A - b B)  x = 0
01038 *
01039             IL2BY2 = .FALSE.
01040             DO 370 J = JE - NW, 1, -1
01041 *
01042 *              If a 2-by-2 block, is in position j-1:j, wait until
01043 *              next iteration to process it (when it will be j:j+1)
01044 *
01045                IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN
01046                   IF( S( J, J-1 ).NE.ZERO ) THEN
01047                      IL2BY2 = .TRUE.
01048                      GO TO 370
01049                   END IF
01050                END IF
01051                BDIAG( 1 ) = P( J, J )
01052                IF( IL2BY2 ) THEN
01053                   NA = 2
01054                   BDIAG( 2 ) = P( J+1, J+1 )
01055                ELSE
01056                   NA = 1
01057                END IF
01058 *
01059 *              Compute x(j) (and x(j+1), if 2-by-2 block)
01060 *
01061                CALL SLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ),
01062      $                      LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ),
01063      $                      N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP,
01064      $                      IINFO )
01065                IF( SCALE.LT.ONE ) THEN
01066 *
01067                   DO 290 JW = 0, NW - 1
01068                      DO 280 JR = 1, JE
01069                         WORK( ( JW+2 )*N+JR ) = SCALE*
01070      $                     WORK( ( JW+2 )*N+JR )
01071   280                CONTINUE
01072   290             CONTINUE
01073                END IF
01074                XMAX = MAX( SCALE*XMAX, TEMP )
01075 *
01076                DO 310 JW = 1, NW
01077                   DO 300 JA = 1, NA
01078                      WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW )
01079   300             CONTINUE
01080   310          CONTINUE
01081 *
01082 *              w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
01083 *
01084                IF( J.GT.1 ) THEN
01085 *
01086 *                 Check whether scaling is necessary for sum.
01087 *
01088                   XSCALE = ONE / MAX( ONE, XMAX )
01089                   TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J )
01090                   IF( IL2BY2 )
01091      $               TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA*
01092      $                      WORK( N+J+1 ) )
01093                   TEMP = MAX( TEMP, ACOEFA, BCOEFA )
01094                   IF( TEMP.GT.BIGNUM*XSCALE ) THEN
01095 *
01096                      DO 330 JW = 0, NW - 1
01097                         DO 320 JR = 1, JE
01098                            WORK( ( JW+2 )*N+JR ) = XSCALE*
01099      $                        WORK( ( JW+2 )*N+JR )
01100   320                   CONTINUE
01101   330                CONTINUE
01102                      XMAX = XMAX*XSCALE
01103                   END IF
01104 *
01105 *                 Compute the contributions of the off-diagonals of
01106 *                 column j (and j+1, if 2-by-2 block) of A and B to the
01107 *                 sums.
01108 *
01109 *
01110                   DO 360 JA = 1, NA
01111                      IF( ILCPLX ) THEN
01112                         CREALA = ACOEF*WORK( 2*N+J+JA-1 )
01113                         CIMAGA = ACOEF*WORK( 3*N+J+JA-1 )
01114                         CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) -
01115      $                           BCOEFI*WORK( 3*N+J+JA-1 )
01116                         CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) +
01117      $                           BCOEFR*WORK( 3*N+J+JA-1 )
01118                         DO 340 JR = 1, J - 1
01119                            WORK( 2*N+JR ) = WORK( 2*N+JR ) -
01120      $                                      CREALA*S( JR, J+JA-1 ) +
01121      $                                      CREALB*P( JR, J+JA-1 )
01122                            WORK( 3*N+JR ) = WORK( 3*N+JR ) -
01123      $                                      CIMAGA*S( JR, J+JA-1 ) +
01124      $                                      CIMAGB*P( JR, J+JA-1 )
01125   340                   CONTINUE
01126                      ELSE
01127                         CREALA = ACOEF*WORK( 2*N+J+JA-1 )
01128                         CREALB = BCOEFR*WORK( 2*N+J+JA-1 )
01129                         DO 350 JR = 1, J - 1
01130                            WORK( 2*N+JR ) = WORK( 2*N+JR ) -
01131      $                                      CREALA*S( JR, J+JA-1 ) +
01132      $                                      CREALB*P( JR, J+JA-1 )
01133   350                   CONTINUE
01134                      END IF
01135   360             CONTINUE
01136                END IF
01137 *
01138                IL2BY2 = .FALSE.
01139   370       CONTINUE
01140 *
01141 *           Copy eigenvector to VR, back transforming if
01142 *           HOWMNY='B'.
01143 *
01144             IEIG = IEIG - NW
01145             IF( ILBACK ) THEN
01146 *
01147                DO 410 JW = 0, NW - 1
01148                   DO 380 JR = 1, N
01149                      WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )*
01150      $                                       VR( JR, 1 )
01151   380             CONTINUE
01152 *
01153 *                 A series of compiler directives to defeat
01154 *                 vectorization for the next loop
01155 *
01156 *
01157                   DO 400 JC = 2, JE
01158                      DO 390 JR = 1, N
01159                         WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) +
01160      $                     WORK( ( JW+2 )*N+JC )*VR( JR, JC )
01161   390                CONTINUE
01162   400             CONTINUE
01163   410          CONTINUE
01164 *
01165                DO 430 JW = 0, NW - 1
01166                   DO 420 JR = 1, N
01167                      VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR )
01168   420             CONTINUE
01169   430          CONTINUE
01170 *
01171                IEND = N
01172             ELSE
01173                DO 450 JW = 0, NW - 1
01174                   DO 440 JR = 1, N
01175                      VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR )
01176   440             CONTINUE
01177   450          CONTINUE
01178 *
01179                IEND = JE
01180             END IF
01181 *
01182 *           Scale eigenvector
01183 *
01184             XMAX = ZERO
01185             IF( ILCPLX ) THEN
01186                DO 460 J = 1, IEND
01187                   XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+
01188      $                   ABS( VR( J, IEIG+1 ) ) )
01189   460          CONTINUE
01190             ELSE
01191                DO 470 J = 1, IEND
01192                   XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) )
01193   470          CONTINUE
01194             END IF
01195 *
01196             IF( XMAX.GT.SAFMIN ) THEN
01197                XSCALE = ONE / XMAX
01198                DO 490 JW = 0, NW - 1
01199                   DO 480 JR = 1, IEND
01200                      VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW )
01201   480             CONTINUE
01202   490          CONTINUE
01203             END IF
01204   500    CONTINUE
01205       END IF
01206 *
01207       RETURN
01208 *
01209 *     End of STGEVC
01210 *
01211       END
 All Files Functions