LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ztgevc.f
Go to the documentation of this file.
00001 *> \brief \b ZTGEVC
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZTGEVC + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgevc.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgevc.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgevc.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
00022 *                          LDVL, VR, LDVR, MM, M, WORK, RWORK, 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 *       DOUBLE PRECISION   RWORK( * )
00031 *       COMPLEX*16         P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
00032 *      $                   VR( LDVR, * ), WORK( * )
00033 *       ..
00034 *  
00035 *  
00036 *
00037 *> \par Purpose:
00038 *  =============
00039 *>
00040 *> \verbatim
00041 *>
00042 *> ZTGEVC computes some or all of the right and/or left eigenvectors of
00043 *> a pair of complex matrices (S,P), where S and P are upper triangular.
00044 *> Matrix pairs of this type are produced by the generalized Schur
00045 *> factorization of a complex matrix pair (A,B):
00046 *> 
00047 *>    A = Q*S*Z**H,  B = Q*P*Z**H
00048 *> 
00049 *> as computed by ZGGHRD + ZHGEQZ.
00050 *> 
00051 *> The right eigenvector x and the left eigenvector y of (S,P)
00052 *> corresponding to an eigenvalue w are defined by:
00053 *> 
00054 *>    S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
00055 *> 
00056 *> where y**H denotes the conjugate tranpose of y.
00057 *> The eigenvalues are not input to this routine, but are computed
00058 *> directly from the diagonal elements of S and P.
00059 *> 
00060 *> This routine returns the matrices X and/or Y of right and left
00061 *> eigenvectors of (S,P), or the products Z*X and/or Q*Y,
00062 *> where Z and Q are input matrices.
00063 *> If Q and Z are the unitary factors from the generalized Schur
00064 *> factorization of a matrix pair (A,B), then Z*X and Q*Y
00065 *> are the matrices of right and left eigenvectors of (A,B).
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.  The eigenvector corresponding to the j-th
00094 *>          eigenvalue is computed if SELECT(j) = .TRUE..
00095 *>          Not referenced if HOWMNY = 'A' or 'B'.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] N
00099 *> \verbatim
00100 *>          N is INTEGER
00101 *>          The order of the matrices S and P.  N >= 0.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] S
00105 *> \verbatim
00106 *>          S is COMPLEX*16 array, dimension (LDS,N)
00107 *>          The upper triangular matrix S from a generalized Schur
00108 *>          factorization, as computed by ZHGEQZ.
00109 *> \endverbatim
00110 *>
00111 *> \param[in] LDS
00112 *> \verbatim
00113 *>          LDS is INTEGER
00114 *>          The leading dimension of array S.  LDS >= max(1,N).
00115 *> \endverbatim
00116 *>
00117 *> \param[in] P
00118 *> \verbatim
00119 *>          P is COMPLEX*16 array, dimension (LDP,N)
00120 *>          The upper triangular matrix P from a generalized Schur
00121 *>          factorization, as computed by ZHGEQZ.  P must have real
00122 *>          diagonal elements.
00123 *> \endverbatim
00124 *>
00125 *> \param[in] LDP
00126 *> \verbatim
00127 *>          LDP is INTEGER
00128 *>          The leading dimension of array P.  LDP >= max(1,N).
00129 *> \endverbatim
00130 *>
00131 *> \param[in,out] VL
00132 *> \verbatim
00133 *>          VL is COMPLEX*16 array, dimension (LDVL,MM)
00134 *>          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
00135 *>          contain an N-by-N matrix Q (usually the unitary matrix Q
00136 *>          of left Schur vectors returned by ZHGEQZ).
00137 *>          On exit, if SIDE = 'L' or 'B', VL contains:
00138 *>          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
00139 *>          if HOWMNY = 'B', the matrix Q*Y;
00140 *>          if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
00141 *>                      SELECT, stored consecutively in the columns of
00142 *>                      VL, in the same order as their eigenvalues.
00143 *>          Not referenced if SIDE = 'R'.
00144 *> \endverbatim
00145 *>
00146 *> \param[in] LDVL
00147 *> \verbatim
00148 *>          LDVL is INTEGER
00149 *>          The leading dimension of array VL.  LDVL >= 1, and if
00150 *>          SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
00151 *> \endverbatim
00152 *>
00153 *> \param[in,out] VR
00154 *> \verbatim
00155 *>          VR is COMPLEX*16 array, dimension (LDVR,MM)
00156 *>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
00157 *>          contain an N-by-N matrix Q (usually the unitary matrix Z
00158 *>          of right Schur vectors returned by ZHGEQZ).
00159 *>          On exit, if SIDE = 'R' or 'B', VR contains:
00160 *>          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
00161 *>          if HOWMNY = 'B', the matrix Z*X;
00162 *>          if HOWMNY = 'S', the right eigenvectors of (S,P) specified by
00163 *>                      SELECT, stored consecutively in the columns of
00164 *>                      VR, in the same order as their eigenvalues.
00165 *>          Not referenced if SIDE = 'L'.
00166 *> \endverbatim
00167 *>
00168 *> \param[in] LDVR
00169 *> \verbatim
00170 *>          LDVR is INTEGER
00171 *>          The leading dimension of the array VR.  LDVR >= 1, and if
00172 *>          SIDE = 'R' or 'B', LDVR >= N.
00173 *> \endverbatim
00174 *>
00175 *> \param[in] MM
00176 *> \verbatim
00177 *>          MM is INTEGER
00178 *>          The number of columns in the arrays VL and/or VR. MM >= M.
00179 *> \endverbatim
00180 *>
00181 *> \param[out] M
00182 *> \verbatim
00183 *>          M is INTEGER
00184 *>          The number of columns in the arrays VL and/or VR actually
00185 *>          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
00186 *>          is set to N.  Each selected eigenvector occupies one column.
00187 *> \endverbatim
00188 *>
00189 *> \param[out] WORK
00190 *> \verbatim
00191 *>          WORK is COMPLEX*16 array, dimension (2*N)
00192 *> \endverbatim
00193 *>
00194 *> \param[out] RWORK
00195 *> \verbatim
00196 *>          RWORK is DOUBLE PRECISION array, dimension (2*N)
00197 *> \endverbatim
00198 *>
00199 *> \param[out] INFO
00200 *> \verbatim
00201 *>          INFO is INTEGER
00202 *>          = 0:  successful exit.
00203 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00204 *> \endverbatim
00205 *
00206 *  Authors:
00207 *  ========
00208 *
00209 *> \author Univ. of Tennessee 
00210 *> \author Univ. of California Berkeley 
00211 *> \author Univ. of Colorado Denver 
00212 *> \author NAG Ltd. 
00213 *
00214 *> \date November 2011
00215 *
00216 *> \ingroup complex16GEcomputational
00217 *
00218 *  =====================================================================
00219       SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
00220      $                   LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
00221 *
00222 *  -- LAPACK computational routine (version 3.4.0) --
00223 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00224 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00225 *     November 2011
00226 *
00227 *     .. Scalar Arguments ..
00228       CHARACTER          HOWMNY, SIDE
00229       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
00230 *     ..
00231 *     .. Array Arguments ..
00232       LOGICAL            SELECT( * )
00233       DOUBLE PRECISION   RWORK( * )
00234       COMPLEX*16         P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
00235      $                   VR( LDVR, * ), WORK( * )
00236 *     ..
00237 *
00238 *
00239 *  =====================================================================
00240 *
00241 *     .. Parameters ..
00242       DOUBLE PRECISION   ZERO, ONE
00243       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00244       COMPLEX*16         CZERO, CONE
00245       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00246      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00247 *     ..
00248 *     .. Local Scalars ..
00249       LOGICAL            COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
00250      $                   LSA, LSB
00251       INTEGER            I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
00252      $                   J, JE, JR
00253       DOUBLE PRECISION   ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
00254      $                   BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
00255      $                   SCALE, SMALL, TEMP, ULP, XMAX
00256       COMPLEX*16         BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
00257 *     ..
00258 *     .. External Functions ..
00259       LOGICAL            LSAME
00260       DOUBLE PRECISION   DLAMCH
00261       COMPLEX*16         ZLADIV
00262       EXTERNAL           LSAME, DLAMCH, ZLADIV
00263 *     ..
00264 *     .. External Subroutines ..
00265       EXTERNAL           DLABAD, XERBLA, ZGEMV
00266 *     ..
00267 *     .. Intrinsic Functions ..
00268       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
00269 *     ..
00270 *     .. Statement Functions ..
00271       DOUBLE PRECISION   ABS1
00272 *     ..
00273 *     .. Statement Function definitions ..
00274       ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
00275 *     ..
00276 *     .. Executable Statements ..
00277 *
00278 *     Decode and Test the input parameters
00279 *
00280       IF( LSAME( HOWMNY, 'A' ) ) THEN
00281          IHWMNY = 1
00282          ILALL = .TRUE.
00283          ILBACK = .FALSE.
00284       ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
00285          IHWMNY = 2
00286          ILALL = .FALSE.
00287          ILBACK = .FALSE.
00288       ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
00289          IHWMNY = 3
00290          ILALL = .TRUE.
00291          ILBACK = .TRUE.
00292       ELSE
00293          IHWMNY = -1
00294       END IF
00295 *
00296       IF( LSAME( SIDE, 'R' ) ) THEN
00297          ISIDE = 1
00298          COMPL = .FALSE.
00299          COMPR = .TRUE.
00300       ELSE IF( LSAME( SIDE, 'L' ) ) THEN
00301          ISIDE = 2
00302          COMPL = .TRUE.
00303          COMPR = .FALSE.
00304       ELSE IF( LSAME( SIDE, 'B' ) ) THEN
00305          ISIDE = 3
00306          COMPL = .TRUE.
00307          COMPR = .TRUE.
00308       ELSE
00309          ISIDE = -1
00310       END IF
00311 *
00312       INFO = 0
00313       IF( ISIDE.LT.0 ) THEN
00314          INFO = -1
00315       ELSE IF( IHWMNY.LT.0 ) THEN
00316          INFO = -2
00317       ELSE IF( N.LT.0 ) THEN
00318          INFO = -4
00319       ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
00320          INFO = -6
00321       ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
00322          INFO = -8
00323       END IF
00324       IF( INFO.NE.0 ) THEN
00325          CALL XERBLA( 'ZTGEVC', -INFO )
00326          RETURN
00327       END IF
00328 *
00329 *     Count the number of eigenvectors
00330 *
00331       IF( .NOT.ILALL ) THEN
00332          IM = 0
00333          DO 10 J = 1, N
00334             IF( SELECT( J ) )
00335      $         IM = IM + 1
00336    10    CONTINUE
00337       ELSE
00338          IM = N
00339       END IF
00340 *
00341 *     Check diagonal of B
00342 *
00343       ILBBAD = .FALSE.
00344       DO 20 J = 1, N
00345          IF( DIMAG( P( J, J ) ).NE.ZERO )
00346      $      ILBBAD = .TRUE.
00347    20 CONTINUE
00348 *
00349       IF( ILBBAD ) THEN
00350          INFO = -7
00351       ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
00352          INFO = -10
00353       ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
00354          INFO = -12
00355       ELSE IF( MM.LT.IM ) THEN
00356          INFO = -13
00357       END IF
00358       IF( INFO.NE.0 ) THEN
00359          CALL XERBLA( 'ZTGEVC', -INFO )
00360          RETURN
00361       END IF
00362 *
00363 *     Quick return if possible
00364 *
00365       M = IM
00366       IF( N.EQ.0 )
00367      $   RETURN
00368 *
00369 *     Machine Constants
00370 *
00371       SAFMIN = DLAMCH( 'Safe minimum' )
00372       BIG = ONE / SAFMIN
00373       CALL DLABAD( SAFMIN, BIG )
00374       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00375       SMALL = SAFMIN*N / ULP
00376       BIG = ONE / SMALL
00377       BIGNUM = ONE / ( SAFMIN*N )
00378 *
00379 *     Compute the 1-norm of each column of the strictly upper triangular
00380 *     part of A and B to check for possible overflow in the triangular
00381 *     solver.
00382 *
00383       ANORM = ABS1( S( 1, 1 ) )
00384       BNORM = ABS1( P( 1, 1 ) )
00385       RWORK( 1 ) = ZERO
00386       RWORK( N+1 ) = ZERO
00387       DO 40 J = 2, N
00388          RWORK( J ) = ZERO
00389          RWORK( N+J ) = ZERO
00390          DO 30 I = 1, J - 1
00391             RWORK( J ) = RWORK( J ) + ABS1( S( I, J ) )
00392             RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) )
00393    30    CONTINUE
00394          ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) )
00395          BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) )
00396    40 CONTINUE
00397 *
00398       ASCALE = ONE / MAX( ANORM, SAFMIN )
00399       BSCALE = ONE / MAX( BNORM, SAFMIN )
00400 *
00401 *     Left eigenvectors
00402 *
00403       IF( COMPL ) THEN
00404          IEIG = 0
00405 *
00406 *        Main loop over eigenvalues
00407 *
00408          DO 140 JE = 1, N
00409             IF( ILALL ) THEN
00410                ILCOMP = .TRUE.
00411             ELSE
00412                ILCOMP = SELECT( JE )
00413             END IF
00414             IF( ILCOMP ) THEN
00415                IEIG = IEIG + 1
00416 *
00417                IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
00418      $             ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
00419 *
00420 *                 Singular matrix pencil -- return unit eigenvector
00421 *
00422                   DO 50 JR = 1, N
00423                      VL( JR, IEIG ) = CZERO
00424    50             CONTINUE
00425                   VL( IEIG, IEIG ) = CONE
00426                   GO TO 140
00427                END IF
00428 *
00429 *              Non-singular eigenvalue:
00430 *              Compute coefficients  a  and  b  in
00431 *                   H
00432 *                 y  ( a A - b B ) = 0
00433 *
00434                TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
00435      $                ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
00436                SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
00437                SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE
00438                ACOEFF = SBETA*ASCALE
00439                BCOEFF = SALPHA*BSCALE
00440 *
00441 *              Scale to avoid underflow
00442 *
00443                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
00444                LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
00445      $               SMALL
00446 *
00447                SCALE = ONE
00448                IF( LSA )
00449      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
00450                IF( LSB )
00451      $            SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
00452      $                    MIN( BNORM, BIG ) )
00453                IF( LSA .OR. LSB ) THEN
00454                   SCALE = MIN( SCALE, ONE /
00455      $                    ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
00456      $                    ABS1( BCOEFF ) ) ) )
00457                   IF( LSA ) THEN
00458                      ACOEFF = ASCALE*( SCALE*SBETA )
00459                   ELSE
00460                      ACOEFF = SCALE*ACOEFF
00461                   END IF
00462                   IF( LSB ) THEN
00463                      BCOEFF = BSCALE*( SCALE*SALPHA )
00464                   ELSE
00465                      BCOEFF = SCALE*BCOEFF
00466                   END IF
00467                END IF
00468 *
00469                ACOEFA = ABS( ACOEFF )
00470                BCOEFA = ABS1( BCOEFF )
00471                XMAX = ONE
00472                DO 60 JR = 1, N
00473                   WORK( JR ) = CZERO
00474    60          CONTINUE
00475                WORK( JE ) = CONE
00476                DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
00477 *
00478 *                                              H
00479 *              Triangular solve of  (a A - b B)  y = 0
00480 *
00481 *                                      H
00482 *              (rowwise in  (a A - b B) , or columnwise in a A - b B)
00483 *
00484                DO 100 J = JE + 1, N
00485 *
00486 *                 Compute
00487 *                       j-1
00488 *                 SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
00489 *                       k=je
00490 *                 (Scale if necessary)
00491 *
00492                   TEMP = ONE / XMAX
00493                   IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM*
00494      $                TEMP ) THEN
00495                      DO 70 JR = JE, J - 1
00496                         WORK( JR ) = TEMP*WORK( JR )
00497    70                CONTINUE
00498                      XMAX = ONE
00499                   END IF
00500                   SUMA = CZERO
00501                   SUMB = CZERO
00502 *
00503                   DO 80 JR = JE, J - 1
00504                      SUMA = SUMA + DCONJG( S( JR, J ) )*WORK( JR )
00505                      SUMB = SUMB + DCONJG( P( JR, J ) )*WORK( JR )
00506    80             CONTINUE
00507                   SUM = ACOEFF*SUMA - DCONJG( BCOEFF )*SUMB
00508 *
00509 *                 Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
00510 *
00511 *                 with scaling and perturbation of the denominator
00512 *
00513                   D = DCONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
00514                   IF( ABS1( D ).LE.DMIN )
00515      $               D = DCMPLX( DMIN )
00516 *
00517                   IF( ABS1( D ).LT.ONE ) THEN
00518                      IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN
00519                         TEMP = ONE / ABS1( SUM )
00520                         DO 90 JR = JE, J - 1
00521                            WORK( JR ) = TEMP*WORK( JR )
00522    90                   CONTINUE
00523                         XMAX = TEMP*XMAX
00524                         SUM = TEMP*SUM
00525                      END IF
00526                   END IF
00527                   WORK( J ) = ZLADIV( -SUM, D )
00528                   XMAX = MAX( XMAX, ABS1( WORK( J ) ) )
00529   100          CONTINUE
00530 *
00531 *              Back transform eigenvector if HOWMNY='B'.
00532 *
00533                IF( ILBACK ) THEN
00534                   CALL ZGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL,
00535      $                        WORK( JE ), 1, CZERO, WORK( N+1 ), 1 )
00536                   ISRC = 2
00537                   IBEG = 1
00538                ELSE
00539                   ISRC = 1
00540                   IBEG = JE
00541                END IF
00542 *
00543 *              Copy and scale eigenvector into column of VL
00544 *
00545                XMAX = ZERO
00546                DO 110 JR = IBEG, N
00547                   XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
00548   110          CONTINUE
00549 *
00550                IF( XMAX.GT.SAFMIN ) THEN
00551                   TEMP = ONE / XMAX
00552                   DO 120 JR = IBEG, N
00553                      VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
00554   120             CONTINUE
00555                ELSE
00556                   IBEG = N + 1
00557                END IF
00558 *
00559                DO 130 JR = 1, IBEG - 1
00560                   VL( JR, IEIG ) = CZERO
00561   130          CONTINUE
00562 *
00563             END IF
00564   140    CONTINUE
00565       END IF
00566 *
00567 *     Right eigenvectors
00568 *
00569       IF( COMPR ) THEN
00570          IEIG = IM + 1
00571 *
00572 *        Main loop over eigenvalues
00573 *
00574          DO 250 JE = N, 1, -1
00575             IF( ILALL ) THEN
00576                ILCOMP = .TRUE.
00577             ELSE
00578                ILCOMP = SELECT( JE )
00579             END IF
00580             IF( ILCOMP ) THEN
00581                IEIG = IEIG - 1
00582 *
00583                IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
00584      $             ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
00585 *
00586 *                 Singular matrix pencil -- return unit eigenvector
00587 *
00588                   DO 150 JR = 1, N
00589                      VR( JR, IEIG ) = CZERO
00590   150             CONTINUE
00591                   VR( IEIG, IEIG ) = CONE
00592                   GO TO 250
00593                END IF
00594 *
00595 *              Non-singular eigenvalue:
00596 *              Compute coefficients  a  and  b  in
00597 *
00598 *              ( a A - b B ) x  = 0
00599 *
00600                TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
00601      $                ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
00602                SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
00603                SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE
00604                ACOEFF = SBETA*ASCALE
00605                BCOEFF = SALPHA*BSCALE
00606 *
00607 *              Scale to avoid underflow
00608 *
00609                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
00610                LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
00611      $               SMALL
00612 *
00613                SCALE = ONE
00614                IF( LSA )
00615      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
00616                IF( LSB )
00617      $            SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
00618      $                    MIN( BNORM, BIG ) )
00619                IF( LSA .OR. LSB ) THEN
00620                   SCALE = MIN( SCALE, ONE /
00621      $                    ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
00622      $                    ABS1( BCOEFF ) ) ) )
00623                   IF( LSA ) THEN
00624                      ACOEFF = ASCALE*( SCALE*SBETA )
00625                   ELSE
00626                      ACOEFF = SCALE*ACOEFF
00627                   END IF
00628                   IF( LSB ) THEN
00629                      BCOEFF = BSCALE*( SCALE*SALPHA )
00630                   ELSE
00631                      BCOEFF = SCALE*BCOEFF
00632                   END IF
00633                END IF
00634 *
00635                ACOEFA = ABS( ACOEFF )
00636                BCOEFA = ABS1( BCOEFF )
00637                XMAX = ONE
00638                DO 160 JR = 1, N
00639                   WORK( JR ) = CZERO
00640   160          CONTINUE
00641                WORK( JE ) = CONE
00642                DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
00643 *
00644 *              Triangular solve of  (a A - b B) x = 0  (columnwise)
00645 *
00646 *              WORK(1:j-1) contains sums w,
00647 *              WORK(j+1:JE) contains x
00648 *
00649                DO 170 JR = 1, JE - 1
00650                   WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE )
00651   170          CONTINUE
00652                WORK( JE ) = CONE
00653 *
00654                DO 210 J = JE - 1, 1, -1
00655 *
00656 *                 Form x(j) := - w(j) / d
00657 *                 with scaling and perturbation of the denominator
00658 *
00659                   D = ACOEFF*S( J, J ) - BCOEFF*P( J, J )
00660                   IF( ABS1( D ).LE.DMIN )
00661      $               D = DCMPLX( DMIN )
00662 *
00663                   IF( ABS1( D ).LT.ONE ) THEN
00664                      IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN
00665                         TEMP = ONE / ABS1( WORK( J ) )
00666                         DO 180 JR = 1, JE
00667                            WORK( JR ) = TEMP*WORK( JR )
00668   180                   CONTINUE
00669                      END IF
00670                   END IF
00671 *
00672                   WORK( J ) = ZLADIV( -WORK( J ), D )
00673 *
00674                   IF( J.GT.1 ) THEN
00675 *
00676 *                    w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
00677 *
00678                      IF( ABS1( WORK( J ) ).GT.ONE ) THEN
00679                         TEMP = ONE / ABS1( WORK( J ) )
00680                         IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE.
00681      $                      BIGNUM*TEMP ) THEN
00682                            DO 190 JR = 1, JE
00683                               WORK( JR ) = TEMP*WORK( JR )
00684   190                      CONTINUE
00685                         END IF
00686                      END IF
00687 *
00688                      CA = ACOEFF*WORK( J )
00689                      CB = BCOEFF*WORK( J )
00690                      DO 200 JR = 1, J - 1
00691                         WORK( JR ) = WORK( JR ) + CA*S( JR, J ) -
00692      $                               CB*P( JR, J )
00693   200                CONTINUE
00694                   END IF
00695   210          CONTINUE
00696 *
00697 *              Back transform eigenvector if HOWMNY='B'.
00698 *
00699                IF( ILBACK ) THEN
00700                   CALL ZGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1,
00701      $                        CZERO, WORK( N+1 ), 1 )
00702                   ISRC = 2
00703                   IEND = N
00704                ELSE
00705                   ISRC = 1
00706                   IEND = JE
00707                END IF
00708 *
00709 *              Copy and scale eigenvector into column of VR
00710 *
00711                XMAX = ZERO
00712                DO 220 JR = 1, IEND
00713                   XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
00714   220          CONTINUE
00715 *
00716                IF( XMAX.GT.SAFMIN ) THEN
00717                   TEMP = ONE / XMAX
00718                   DO 230 JR = 1, IEND
00719                      VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
00720   230             CONTINUE
00721                ELSE
00722                   IEND = 0
00723                END IF
00724 *
00725                DO 240 JR = IEND + 1, N
00726                   VR( JR, IEIG ) = CZERO
00727   240          CONTINUE
00728 *
00729             END IF
00730   250    CONTINUE
00731       END IF
00732 *
00733       RETURN
00734 *
00735 *     End of ZTGEVC
00736 *
00737       END
 All Files Functions