Annotation of rpl/lapack/lapack/dhgeqz.f, revision 1.18

1.9       bertrand    1: *> \brief \b DHGEQZ
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.18    ! bertrand    5: * Online html documentation available at
        !             6: *            http://www.netlib.org/lapack/explore-html/
1.9       bertrand    7: *
                      8: *> \htmlonly
1.18    ! bertrand    9: *> Download DHGEQZ + dependencies
        !            10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dhgeqz.f">
        !            11: *> [TGZ]</a>
        !            12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dhgeqz.f">
        !            13: *> [ZIP]</a>
        !            14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dhgeqz.f">
1.9       bertrand   15: *> [TXT]</a>
1.18    ! bertrand   16: *> \endhtmlonly
1.9       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
                     22: *                          ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
                     23: *                          LWORK, INFO )
1.18    ! bertrand   24: *
1.9       bertrand   25: *       .. Scalar Arguments ..
                     26: *       CHARACTER          COMPQ, COMPZ, JOB
                     27: *       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
                     28: *       ..
                     29: *       .. Array Arguments ..
                     30: *       DOUBLE PRECISION   ALPHAI( * ), ALPHAR( * ), BETA( * ),
                     31: *      $                   H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
                     32: *      $                   WORK( * ), Z( LDZ, * )
                     33: *       ..
1.18    ! bertrand   34: *
1.9       bertrand   35: *
                     36: *> \par Purpose:
                     37: *  =============
                     38: *>
                     39: *> \verbatim
                     40: *>
                     41: *> DHGEQZ computes the eigenvalues of a real matrix pair (H,T),
                     42: *> where H is an upper Hessenberg matrix and T is upper triangular,
                     43: *> using the double-shift QZ method.
                     44: *> Matrix pairs of this type are produced by the reduction to
                     45: *> generalized upper Hessenberg form of a real matrix pair (A,B):
                     46: *>
                     47: *>    A = Q1*H*Z1**T,  B = Q1*T*Z1**T,
                     48: *>
                     49: *> as computed by DGGHRD.
                     50: *>
                     51: *> If JOB='S', then the Hessenberg-triangular pair (H,T) is
                     52: *> also reduced to generalized Schur form,
1.18    ! bertrand   53: *>
1.9       bertrand   54: *>    H = Q*S*Z**T,  T = Q*P*Z**T,
1.18    ! bertrand   55: *>
1.9       bertrand   56: *> where Q and Z are orthogonal matrices, P is an upper triangular
                     57: *> matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
                     58: *> diagonal blocks.
                     59: *>
                     60: *> The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
                     61: *> (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
                     62: *> eigenvalues.
                     63: *>
                     64: *> Additionally, the 2-by-2 upper triangular diagonal blocks of P
                     65: *> corresponding to 2-by-2 blocks of S are reduced to positive diagonal
                     66: *> form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
                     67: *> P(j,j) > 0, and P(j+1,j+1) > 0.
                     68: *>
                     69: *> Optionally, the orthogonal matrix Q from the generalized Schur
                     70: *> factorization may be postmultiplied into an input matrix Q1, and the
                     71: *> orthogonal matrix Z may be postmultiplied into an input matrix Z1.
                     72: *> If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced
                     73: *> the matrix pair (A,B) to generalized upper Hessenberg form, then the
                     74: *> output matrices Q1*Q and Z1*Z are the orthogonal factors from the
                     75: *> generalized Schur factorization of (A,B):
                     76: *>
                     77: *>    A = (Q1*Q)*S*(Z1*Z)**T,  B = (Q1*Q)*P*(Z1*Z)**T.
1.18    ! bertrand   78: *>
1.9       bertrand   79: *> To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
                     80: *> of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
                     81: *> complex and beta real.
                     82: *> If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
                     83: *> generalized nonsymmetric eigenvalue problem (GNEP)
                     84: *>    A*x = lambda*B*x
                     85: *> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
                     86: *> alternate form of the GNEP
                     87: *>    mu*A*y = B*y.
                     88: *> Real eigenvalues can be read directly from the generalized Schur
1.18    ! bertrand   89: *> form:
1.9       bertrand   90: *>   alpha = S(i,i), beta = P(i,i).
                     91: *>
                     92: *> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
                     93: *>      Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
                     94: *>      pp. 241--256.
                     95: *> \endverbatim
                     96: *
                     97: *  Arguments:
                     98: *  ==========
                     99: *
                    100: *> \param[in] JOB
                    101: *> \verbatim
                    102: *>          JOB is CHARACTER*1
                    103: *>          = 'E': Compute eigenvalues only;
1.18    ! bertrand  104: *>          = 'S': Compute eigenvalues and the Schur form.
1.9       bertrand  105: *> \endverbatim
                    106: *>
                    107: *> \param[in] COMPQ
                    108: *> \verbatim
                    109: *>          COMPQ is CHARACTER*1
                    110: *>          = 'N': Left Schur vectors (Q) are not computed;
                    111: *>          = 'I': Q is initialized to the unit matrix and the matrix Q
                    112: *>                 of left Schur vectors of (H,T) is returned;
                    113: *>          = 'V': Q must contain an orthogonal matrix Q1 on entry and
                    114: *>                 the product Q1*Q is returned.
                    115: *> \endverbatim
                    116: *>
                    117: *> \param[in] COMPZ
                    118: *> \verbatim
                    119: *>          COMPZ is CHARACTER*1
                    120: *>          = 'N': Right Schur vectors (Z) are not computed;
                    121: *>          = 'I': Z is initialized to the unit matrix and the matrix Z
                    122: *>                 of right Schur vectors of (H,T) is returned;
                    123: *>          = 'V': Z must contain an orthogonal matrix Z1 on entry and
                    124: *>                 the product Z1*Z is returned.
                    125: *> \endverbatim
                    126: *>
                    127: *> \param[in] N
                    128: *> \verbatim
                    129: *>          N is INTEGER
                    130: *>          The order of the matrices H, T, Q, and Z.  N >= 0.
                    131: *> \endverbatim
                    132: *>
                    133: *> \param[in] ILO
                    134: *> \verbatim
                    135: *>          ILO is INTEGER
                    136: *> \endverbatim
                    137: *>
                    138: *> \param[in] IHI
                    139: *> \verbatim
                    140: *>          IHI is INTEGER
                    141: *>          ILO and IHI mark the rows and columns of H which are in
                    142: *>          Hessenberg form.  It is assumed that A is already upper
                    143: *>          triangular in rows and columns 1:ILO-1 and IHI+1:N.
                    144: *>          If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
                    145: *> \endverbatim
                    146: *>
                    147: *> \param[in,out] H
                    148: *> \verbatim
                    149: *>          H is DOUBLE PRECISION array, dimension (LDH, N)
                    150: *>          On entry, the N-by-N upper Hessenberg matrix H.
                    151: *>          On exit, if JOB = 'S', H contains the upper quasi-triangular
                    152: *>          matrix S from the generalized Schur factorization.
                    153: *>          If JOB = 'E', the diagonal blocks of H match those of S, but
                    154: *>          the rest of H is unspecified.
                    155: *> \endverbatim
                    156: *>
                    157: *> \param[in] LDH
                    158: *> \verbatim
                    159: *>          LDH is INTEGER
                    160: *>          The leading dimension of the array H.  LDH >= max( 1, N ).
                    161: *> \endverbatim
                    162: *>
                    163: *> \param[in,out] T
                    164: *> \verbatim
                    165: *>          T is DOUBLE PRECISION array, dimension (LDT, N)
                    166: *>          On entry, the N-by-N upper triangular matrix T.
                    167: *>          On exit, if JOB = 'S', T contains the upper triangular
                    168: *>          matrix P from the generalized Schur factorization;
                    169: *>          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
                    170: *>          are reduced to positive diagonal form, i.e., if H(j+1,j) is
                    171: *>          non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and
                    172: *>          T(j+1,j+1) > 0.
                    173: *>          If JOB = 'E', the diagonal blocks of T match those of P, but
                    174: *>          the rest of T is unspecified.
                    175: *> \endverbatim
                    176: *>
                    177: *> \param[in] LDT
                    178: *> \verbatim
                    179: *>          LDT is INTEGER
                    180: *>          The leading dimension of the array T.  LDT >= max( 1, N ).
                    181: *> \endverbatim
                    182: *>
                    183: *> \param[out] ALPHAR
                    184: *> \verbatim
                    185: *>          ALPHAR is DOUBLE PRECISION array, dimension (N)
                    186: *>          The real parts of each scalar alpha defining an eigenvalue
                    187: *>          of GNEP.
                    188: *> \endverbatim
                    189: *>
                    190: *> \param[out] ALPHAI
                    191: *> \verbatim
                    192: *>          ALPHAI is DOUBLE PRECISION array, dimension (N)
                    193: *>          The imaginary parts of each scalar alpha defining an
                    194: *>          eigenvalue of GNEP.
                    195: *>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
                    196: *>          positive, then the j-th and (j+1)-st eigenvalues are a
                    197: *>          complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
                    198: *> \endverbatim
                    199: *>
                    200: *> \param[out] BETA
                    201: *> \verbatim
                    202: *>          BETA is DOUBLE PRECISION array, dimension (N)
                    203: *>          The scalars beta that define the eigenvalues of GNEP.
                    204: *>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
                    205: *>          beta = BETA(j) represent the j-th eigenvalue of the matrix
                    206: *>          pair (A,B), in one of the forms lambda = alpha/beta or
                    207: *>          mu = beta/alpha.  Since either lambda or mu may overflow,
                    208: *>          they should not, in general, be computed.
                    209: *> \endverbatim
                    210: *>
                    211: *> \param[in,out] Q
                    212: *> \verbatim
                    213: *>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
1.16      bertrand  214: *>          On entry, if COMPQ = 'V', the orthogonal matrix Q1 used in
1.9       bertrand  215: *>          the reduction of (A,B) to generalized Hessenberg form.
1.16      bertrand  216: *>          On exit, if COMPQ = 'I', the orthogonal matrix of left Schur
                    217: *>          vectors of (H,T), and if COMPQ = 'V', the orthogonal matrix
1.9       bertrand  218: *>          of left Schur vectors of (A,B).
1.16      bertrand  219: *>          Not referenced if COMPQ = 'N'.
1.9       bertrand  220: *> \endverbatim
                    221: *>
                    222: *> \param[in] LDQ
                    223: *> \verbatim
                    224: *>          LDQ is INTEGER
                    225: *>          The leading dimension of the array Q.  LDQ >= 1.
                    226: *>          If COMPQ='V' or 'I', then LDQ >= N.
                    227: *> \endverbatim
                    228: *>
                    229: *> \param[in,out] Z
                    230: *> \verbatim
                    231: *>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
                    232: *>          On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
                    233: *>          the reduction of (A,B) to generalized Hessenberg form.
                    234: *>          On exit, if COMPZ = 'I', the orthogonal matrix of
                    235: *>          right Schur vectors of (H,T), and if COMPZ = 'V', the
                    236: *>          orthogonal matrix of right Schur vectors of (A,B).
                    237: *>          Not referenced if COMPZ = 'N'.
                    238: *> \endverbatim
                    239: *>
                    240: *> \param[in] LDZ
                    241: *> \verbatim
                    242: *>          LDZ is INTEGER
                    243: *>          The leading dimension of the array Z.  LDZ >= 1.
                    244: *>          If COMPZ='V' or 'I', then LDZ >= N.
                    245: *> \endverbatim
                    246: *>
                    247: *> \param[out] WORK
                    248: *> \verbatim
                    249: *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
                    250: *>          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
                    251: *> \endverbatim
                    252: *>
                    253: *> \param[in] LWORK
                    254: *> \verbatim
                    255: *>          LWORK is INTEGER
                    256: *>          The dimension of the array WORK.  LWORK >= max(1,N).
                    257: *>
                    258: *>          If LWORK = -1, then a workspace query is assumed; the routine
                    259: *>          only calculates the optimal size of the WORK array, returns
                    260: *>          this value as the first entry of the WORK array, and no error
                    261: *>          message related to LWORK is issued by XERBLA.
                    262: *> \endverbatim
                    263: *>
                    264: *> \param[out] INFO
                    265: *> \verbatim
                    266: *>          INFO is INTEGER
                    267: *>          = 0: successful exit
                    268: *>          < 0: if INFO = -i, the i-th argument had an illegal value
                    269: *>          = 1,...,N: the QZ iteration did not converge.  (H,T) is not
                    270: *>                     in Schur form, but ALPHAR(i), ALPHAI(i), and
                    271: *>                     BETA(i), i=INFO+1,...,N should be correct.
                    272: *>          = N+1,...,2*N: the shift calculation failed.  (H,T) is not
                    273: *>                     in Schur form, but ALPHAR(i), ALPHAI(i), and
                    274: *>                     BETA(i), i=INFO-N+1,...,N should be correct.
                    275: *> \endverbatim
                    276: *
                    277: *  Authors:
                    278: *  ========
                    279: *
1.18    ! bertrand  280: *> \author Univ. of Tennessee
        !           281: *> \author Univ. of California Berkeley
        !           282: *> \author Univ. of Colorado Denver
        !           283: *> \author NAG Ltd.
1.9       bertrand  284: *
1.16      bertrand  285: *> \date June 2016
1.9       bertrand  286: *
                    287: *> \ingroup doubleGEcomputational
                    288: *
                    289: *> \par Further Details:
                    290: *  =====================
                    291: *>
                    292: *> \verbatim
                    293: *>
                    294: *>  Iteration counters:
                    295: *>
                    296: *>  JITER  -- counts iterations.
                    297: *>  IITER  -- counts iterations run since ILAST was last
                    298: *>            changed.  This is therefore reset only when a 1-by-1 or
                    299: *>            2-by-2 block deflates off the bottom.
                    300: *> \endverbatim
                    301: *>
                    302: *  =====================================================================
1.1       bertrand  303:       SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
                    304:      $                   ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
                    305:      $                   LWORK, INFO )
                    306: *
1.18    ! bertrand  307: *  -- LAPACK computational routine (version 3.7.0) --
1.1       bertrand  308: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    309: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.16      bertrand  310: *     June 2016
1.1       bertrand  311: *
                    312: *     .. Scalar Arguments ..
                    313:       CHARACTER          COMPQ, COMPZ, JOB
                    314:       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
                    315: *     ..
                    316: *     .. Array Arguments ..
                    317:       DOUBLE PRECISION   ALPHAI( * ), ALPHAR( * ), BETA( * ),
                    318:      $                   H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
                    319:      $                   WORK( * ), Z( LDZ, * )
                    320: *     ..
                    321: *
                    322: *  =====================================================================
                    323: *
                    324: *     .. Parameters ..
                    325: *    $                     SAFETY = 1.0E+0 )
                    326:       DOUBLE PRECISION   HALF, ZERO, ONE, SAFETY
                    327:       PARAMETER          ( HALF = 0.5D+0, ZERO = 0.0D+0, ONE = 1.0D+0,
                    328:      $                   SAFETY = 1.0D+2 )
                    329: *     ..
                    330: *     .. Local Scalars ..
                    331:       LOGICAL            ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ,
                    332:      $                   LQUERY
                    333:       INTEGER            ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
                    334:      $                   ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
                    335:      $                   JR, MAXIT
                    336:       DOUBLE PRECISION   A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
                    337:      $                   AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
                    338:      $                   AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I,
                    339:      $                   B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE,
                    340:      $                   BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL,
                    341:      $                   CQ, CR, CZ, ESHIFT, S, S1, S1INV, S2, SAFMAX,
                    342:      $                   SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, SZR, T1,
                    343:      $                   TAU, TEMP, TEMP2, TEMPI, TEMPR, U1, U12, U12L,
                    344:      $                   U2, ULP, VS, W11, W12, W21, W22, WABS, WI, WR,
                    345:      $                   WR2
                    346: *     ..
                    347: *     .. Local Arrays ..
                    348:       DOUBLE PRECISION   V( 3 )
                    349: *     ..
                    350: *     .. External Functions ..
                    351:       LOGICAL            LSAME
                    352:       DOUBLE PRECISION   DLAMCH, DLANHS, DLAPY2, DLAPY3
                    353:       EXTERNAL           LSAME, DLAMCH, DLANHS, DLAPY2, DLAPY3
                    354: *     ..
                    355: *     .. External Subroutines ..
                    356:       EXTERNAL           DLAG2, DLARFG, DLARTG, DLASET, DLASV2, DROT,
                    357:      $                   XERBLA
                    358: *     ..
                    359: *     .. Intrinsic Functions ..
                    360:       INTRINSIC          ABS, DBLE, MAX, MIN, SQRT
                    361: *     ..
                    362: *     .. Executable Statements ..
                    363: *
                    364: *     Decode JOB, COMPQ, COMPZ
                    365: *
                    366:       IF( LSAME( JOB, 'E' ) ) THEN
                    367:          ILSCHR = .FALSE.
                    368:          ISCHUR = 1
                    369:       ELSE IF( LSAME( JOB, 'S' ) ) THEN
                    370:          ILSCHR = .TRUE.
                    371:          ISCHUR = 2
                    372:       ELSE
                    373:          ISCHUR = 0
                    374:       END IF
                    375: *
                    376:       IF( LSAME( COMPQ, 'N' ) ) THEN
                    377:          ILQ = .FALSE.
                    378:          ICOMPQ = 1
                    379:       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
                    380:          ILQ = .TRUE.
                    381:          ICOMPQ = 2
                    382:       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
                    383:          ILQ = .TRUE.
                    384:          ICOMPQ = 3
                    385:       ELSE
                    386:          ICOMPQ = 0
                    387:       END IF
                    388: *
                    389:       IF( LSAME( COMPZ, 'N' ) ) THEN
                    390:          ILZ = .FALSE.
                    391:          ICOMPZ = 1
                    392:       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
                    393:          ILZ = .TRUE.
                    394:          ICOMPZ = 2
                    395:       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
                    396:          ILZ = .TRUE.
                    397:          ICOMPZ = 3
                    398:       ELSE
                    399:          ICOMPZ = 0
                    400:       END IF
                    401: *
                    402: *     Check Argument Values
                    403: *
                    404:       INFO = 0
                    405:       WORK( 1 ) = MAX( 1, N )
                    406:       LQUERY = ( LWORK.EQ.-1 )
                    407:       IF( ISCHUR.EQ.0 ) THEN
                    408:          INFO = -1
                    409:       ELSE IF( ICOMPQ.EQ.0 ) THEN
                    410:          INFO = -2
                    411:       ELSE IF( ICOMPZ.EQ.0 ) THEN
                    412:          INFO = -3
                    413:       ELSE IF( N.LT.0 ) THEN
                    414:          INFO = -4
                    415:       ELSE IF( ILO.LT.1 ) THEN
                    416:          INFO = -5
                    417:       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
                    418:          INFO = -6
                    419:       ELSE IF( LDH.LT.N ) THEN
                    420:          INFO = -8
                    421:       ELSE IF( LDT.LT.N ) THEN
                    422:          INFO = -10
                    423:       ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
                    424:          INFO = -15
                    425:       ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
                    426:          INFO = -17
                    427:       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
                    428:          INFO = -19
                    429:       END IF
                    430:       IF( INFO.NE.0 ) THEN
                    431:          CALL XERBLA( 'DHGEQZ', -INFO )
                    432:          RETURN
                    433:       ELSE IF( LQUERY ) THEN
                    434:          RETURN
                    435:       END IF
                    436: *
                    437: *     Quick return if possible
                    438: *
                    439:       IF( N.LE.0 ) THEN
                    440:          WORK( 1 ) = DBLE( 1 )
                    441:          RETURN
                    442:       END IF
                    443: *
                    444: *     Initialize Q and Z
                    445: *
                    446:       IF( ICOMPQ.EQ.3 )
                    447:      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
                    448:       IF( ICOMPZ.EQ.3 )
                    449:      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
                    450: *
                    451: *     Machine Constants
                    452: *
                    453:       IN = IHI + 1 - ILO
                    454:       SAFMIN = DLAMCH( 'S' )
                    455:       SAFMAX = ONE / SAFMIN
                    456:       ULP = DLAMCH( 'E' )*DLAMCH( 'B' )
                    457:       ANORM = DLANHS( 'F', IN, H( ILO, ILO ), LDH, WORK )
                    458:       BNORM = DLANHS( 'F', IN, T( ILO, ILO ), LDT, WORK )
                    459:       ATOL = MAX( SAFMIN, ULP*ANORM )
                    460:       BTOL = MAX( SAFMIN, ULP*BNORM )
                    461:       ASCALE = ONE / MAX( SAFMIN, ANORM )
                    462:       BSCALE = ONE / MAX( SAFMIN, BNORM )
                    463: *
                    464: *     Set Eigenvalues IHI+1:N
                    465: *
                    466:       DO 30 J = IHI + 1, N
                    467:          IF( T( J, J ).LT.ZERO ) THEN
                    468:             IF( ILSCHR ) THEN
                    469:                DO 10 JR = 1, J
                    470:                   H( JR, J ) = -H( JR, J )
                    471:                   T( JR, J ) = -T( JR, J )
                    472:    10          CONTINUE
                    473:             ELSE
                    474:                H( J, J ) = -H( J, J )
                    475:                T( J, J ) = -T( J, J )
                    476:             END IF
                    477:             IF( ILZ ) THEN
                    478:                DO 20 JR = 1, N
                    479:                   Z( JR, J ) = -Z( JR, J )
                    480:    20          CONTINUE
                    481:             END IF
                    482:          END IF
                    483:          ALPHAR( J ) = H( J, J )
                    484:          ALPHAI( J ) = ZERO
                    485:          BETA( J ) = T( J, J )
                    486:    30 CONTINUE
                    487: *
                    488: *     If IHI < ILO, skip QZ steps
                    489: *
                    490:       IF( IHI.LT.ILO )
                    491:      $   GO TO 380
                    492: *
                    493: *     MAIN QZ ITERATION LOOP
                    494: *
                    495: *     Initialize dynamic indices
                    496: *
                    497: *     Eigenvalues ILAST+1:N have been found.
                    498: *        Column operations modify rows IFRSTM:whatever.
                    499: *        Row operations modify columns whatever:ILASTM.
                    500: *
                    501: *     If only eigenvalues are being computed, then
                    502: *        IFRSTM is the row of the last splitting row above row ILAST;
                    503: *        this is always at least ILO.
                    504: *     IITER counts iterations since the last eigenvalue was found,
                    505: *        to tell when to use an extraordinary shift.
                    506: *     MAXIT is the maximum number of QZ sweeps allowed.
                    507: *
                    508:       ILAST = IHI
                    509:       IF( ILSCHR ) THEN
                    510:          IFRSTM = 1
                    511:          ILASTM = N
                    512:       ELSE
                    513:          IFRSTM = ILO
                    514:          ILASTM = IHI
                    515:       END IF
                    516:       IITER = 0
                    517:       ESHIFT = ZERO
                    518:       MAXIT = 30*( IHI-ILO+1 )
                    519: *
                    520:       DO 360 JITER = 1, MAXIT
                    521: *
                    522: *        Split the matrix if possible.
                    523: *
                    524: *        Two tests:
                    525: *           1: H(j,j-1)=0  or  j=ILO
                    526: *           2: T(j,j)=0
                    527: *
                    528:          IF( ILAST.EQ.ILO ) THEN
                    529: *
                    530: *           Special case: j=ILAST
                    531: *
                    532:             GO TO 80
                    533:          ELSE
                    534:             IF( ABS( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
                    535:                H( ILAST, ILAST-1 ) = ZERO
                    536:                GO TO 80
                    537:             END IF
                    538:          END IF
                    539: *
                    540:          IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
                    541:             T( ILAST, ILAST ) = ZERO
                    542:             GO TO 70
                    543:          END IF
                    544: *
                    545: *        General case: j<ILAST
                    546: *
                    547:          DO 60 J = ILAST - 1, ILO, -1
                    548: *
                    549: *           Test 1: for H(j,j-1)=0 or j=ILO
                    550: *
                    551:             IF( J.EQ.ILO ) THEN
                    552:                ILAZRO = .TRUE.
                    553:             ELSE
                    554:                IF( ABS( H( J, J-1 ) ).LE.ATOL ) THEN
                    555:                   H( J, J-1 ) = ZERO
                    556:                   ILAZRO = .TRUE.
                    557:                ELSE
                    558:                   ILAZRO = .FALSE.
                    559:                END IF
                    560:             END IF
                    561: *
                    562: *           Test 2: for T(j,j)=0
                    563: *
                    564:             IF( ABS( T( J, J ) ).LT.BTOL ) THEN
                    565:                T( J, J ) = ZERO
                    566: *
                    567: *              Test 1a: Check for 2 consecutive small subdiagonals in A
                    568: *
                    569:                ILAZR2 = .FALSE.
                    570:                IF( .NOT.ILAZRO ) THEN
                    571:                   TEMP = ABS( H( J, J-1 ) )
                    572:                   TEMP2 = ABS( H( J, J ) )
                    573:                   TEMPR = MAX( TEMP, TEMP2 )
                    574:                   IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
                    575:                      TEMP = TEMP / TEMPR
                    576:                      TEMP2 = TEMP2 / TEMPR
                    577:                   END IF
                    578:                   IF( TEMP*( ASCALE*ABS( H( J+1, J ) ) ).LE.TEMP2*
                    579:      $                ( ASCALE*ATOL ) )ILAZR2 = .TRUE.
                    580:                END IF
                    581: *
                    582: *              If both tests pass (1 & 2), i.e., the leading diagonal
                    583: *              element of B in the block is zero, split a 1x1 block off
                    584: *              at the top. (I.e., at the J-th row/column) The leading
                    585: *              diagonal element of the remainder can also be zero, so
                    586: *              this may have to be done repeatedly.
                    587: *
                    588:                IF( ILAZRO .OR. ILAZR2 ) THEN
                    589:                   DO 40 JCH = J, ILAST - 1
                    590:                      TEMP = H( JCH, JCH )
                    591:                      CALL DLARTG( TEMP, H( JCH+1, JCH ), C, S,
                    592:      $                            H( JCH, JCH ) )
                    593:                      H( JCH+1, JCH ) = ZERO
                    594:                      CALL DROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
                    595:      $                          H( JCH+1, JCH+1 ), LDH, C, S )
                    596:                      CALL DROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
                    597:      $                          T( JCH+1, JCH+1 ), LDT, C, S )
                    598:                      IF( ILQ )
                    599:      $                  CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
                    600:      $                             C, S )
                    601:                      IF( ILAZR2 )
                    602:      $                  H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
                    603:                      ILAZR2 = .FALSE.
                    604:                      IF( ABS( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
                    605:                         IF( JCH+1.GE.ILAST ) THEN
                    606:                            GO TO 80
                    607:                         ELSE
                    608:                            IFIRST = JCH + 1
                    609:                            GO TO 110
                    610:                         END IF
                    611:                      END IF
                    612:                      T( JCH+1, JCH+1 ) = ZERO
                    613:    40             CONTINUE
                    614:                   GO TO 70
                    615:                ELSE
                    616: *
                    617: *                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
                    618: *                 Then process as in the case T(ILAST,ILAST)=0
                    619: *
                    620:                   DO 50 JCH = J, ILAST - 1
                    621:                      TEMP = T( JCH, JCH+1 )
                    622:                      CALL DLARTG( TEMP, T( JCH+1, JCH+1 ), C, S,
                    623:      $                            T( JCH, JCH+1 ) )
                    624:                      T( JCH+1, JCH+1 ) = ZERO
                    625:                      IF( JCH.LT.ILASTM-1 )
                    626:      $                  CALL DROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
                    627:      $                             T( JCH+1, JCH+2 ), LDT, C, S )
                    628:                      CALL DROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
                    629:      $                          H( JCH+1, JCH-1 ), LDH, C, S )
                    630:                      IF( ILQ )
                    631:      $                  CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
                    632:      $                             C, S )
                    633:                      TEMP = H( JCH+1, JCH )
                    634:                      CALL DLARTG( TEMP, H( JCH+1, JCH-1 ), C, S,
                    635:      $                            H( JCH+1, JCH ) )
                    636:                      H( JCH+1, JCH-1 ) = ZERO
                    637:                      CALL DROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
                    638:      $                          H( IFRSTM, JCH-1 ), 1, C, S )
                    639:                      CALL DROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
                    640:      $                          T( IFRSTM, JCH-1 ), 1, C, S )
                    641:                      IF( ILZ )
                    642:      $                  CALL DROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
                    643:      $                             C, S )
                    644:    50             CONTINUE
                    645:                   GO TO 70
                    646:                END IF
                    647:             ELSE IF( ILAZRO ) THEN
                    648: *
                    649: *              Only test 1 passed -- work on J:ILAST
                    650: *
                    651:                IFIRST = J
                    652:                GO TO 110
                    653:             END IF
                    654: *
                    655: *           Neither test passed -- try next J
                    656: *
                    657:    60    CONTINUE
                    658: *
                    659: *        (Drop-through is "impossible")
                    660: *
                    661:          INFO = N + 1
                    662:          GO TO 420
                    663: *
                    664: *        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
                    665: *        1x1 block.
                    666: *
                    667:    70    CONTINUE
                    668:          TEMP = H( ILAST, ILAST )
                    669:          CALL DLARTG( TEMP, H( ILAST, ILAST-1 ), C, S,
                    670:      $                H( ILAST, ILAST ) )
                    671:          H( ILAST, ILAST-1 ) = ZERO
                    672:          CALL DROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
                    673:      $              H( IFRSTM, ILAST-1 ), 1, C, S )
                    674:          CALL DROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
                    675:      $              T( IFRSTM, ILAST-1 ), 1, C, S )
                    676:          IF( ILZ )
                    677:      $      CALL DROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
                    678: *
                    679: *        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
                    680: *                              and BETA
                    681: *
                    682:    80    CONTINUE
                    683:          IF( T( ILAST, ILAST ).LT.ZERO ) THEN
                    684:             IF( ILSCHR ) THEN
                    685:                DO 90 J = IFRSTM, ILAST
                    686:                   H( J, ILAST ) = -H( J, ILAST )
                    687:                   T( J, ILAST ) = -T( J, ILAST )
                    688:    90          CONTINUE
                    689:             ELSE
                    690:                H( ILAST, ILAST ) = -H( ILAST, ILAST )
                    691:                T( ILAST, ILAST ) = -T( ILAST, ILAST )
                    692:             END IF
                    693:             IF( ILZ ) THEN
                    694:                DO 100 J = 1, N
                    695:                   Z( J, ILAST ) = -Z( J, ILAST )
                    696:   100          CONTINUE
                    697:             END IF
                    698:          END IF
                    699:          ALPHAR( ILAST ) = H( ILAST, ILAST )
                    700:          ALPHAI( ILAST ) = ZERO
                    701:          BETA( ILAST ) = T( ILAST, ILAST )
                    702: *
                    703: *        Go to next block -- exit if finished.
                    704: *
                    705:          ILAST = ILAST - 1
                    706:          IF( ILAST.LT.ILO )
                    707:      $      GO TO 380
                    708: *
                    709: *        Reset counters
                    710: *
                    711:          IITER = 0
                    712:          ESHIFT = ZERO
                    713:          IF( .NOT.ILSCHR ) THEN
                    714:             ILASTM = ILAST
                    715:             IF( IFRSTM.GT.ILAST )
                    716:      $         IFRSTM = ILO
                    717:          END IF
                    718:          GO TO 350
                    719: *
                    720: *        QZ step
                    721: *
                    722: *        This iteration only involves rows/columns IFIRST:ILAST. We
                    723: *        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
                    724: *
                    725:   110    CONTINUE
                    726:          IITER = IITER + 1
                    727:          IF( .NOT.ILSCHR ) THEN
                    728:             IFRSTM = IFIRST
                    729:          END IF
                    730: *
                    731: *        Compute single shifts.
                    732: *
                    733: *        At this point, IFIRST < ILAST, and the diagonal elements of
                    734: *        T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
                    735: *        magnitude)
                    736: *
                    737:          IF( ( IITER / 10 )*10.EQ.IITER ) THEN
                    738: *
                    739: *           Exceptional shift.  Chosen for no particularly good reason.
                    740: *           (Single shift only.)
                    741: *
1.14      bertrand  742:             IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST, ILAST-1 ) ).LT.
1.1       bertrand  743:      $          ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
1.14      bertrand  744:                ESHIFT = H( ILAST, ILAST-1 ) /
1.1       bertrand  745:      $                  T( ILAST-1, ILAST-1 )
                    746:             ELSE
                    747:                ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) )
                    748:             END IF
                    749:             S1 = ONE
                    750:             WR = ESHIFT
                    751: *
                    752:          ELSE
                    753: *
                    754: *           Shifts based on the generalized eigenvalues of the
                    755: *           bottom-right 2x2 block of A and B. The first eigenvalue
                    756: *           returned by DLAG2 is the Wilkinson shift (AEP p.512),
                    757: *
                    758:             CALL DLAG2( H( ILAST-1, ILAST-1 ), LDH,
                    759:      $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
                    760:      $                  S2, WR, WR2, WI )
                    761: *
1.14      bertrand  762:             IF ( ABS( (WR/S1)*T( ILAST, ILAST ) - H( ILAST, ILAST ) )
1.18    ! bertrand  763:      $         .GT. ABS( (WR2/S2)*T( ILAST, ILAST )
1.14      bertrand  764:      $         - H( ILAST, ILAST ) ) ) THEN
                    765:                TEMP = WR
                    766:                WR = WR2
                    767:                WR2 = TEMP
                    768:                TEMP = S1
                    769:                S1 = S2
                    770:                S2 = TEMP
                    771:             END IF
1.1       bertrand  772:             TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) )
                    773:             IF( WI.NE.ZERO )
                    774:      $         GO TO 200
                    775:          END IF
                    776: *
                    777: *        Fiddle with shift to avoid overflow
                    778: *
                    779:          TEMP = MIN( ASCALE, ONE )*( HALF*SAFMAX )
                    780:          IF( S1.GT.TEMP ) THEN
                    781:             SCALE = TEMP / S1
                    782:          ELSE
                    783:             SCALE = ONE
                    784:          END IF
                    785: *
                    786:          TEMP = MIN( BSCALE, ONE )*( HALF*SAFMAX )
                    787:          IF( ABS( WR ).GT.TEMP )
                    788:      $      SCALE = MIN( SCALE, TEMP / ABS( WR ) )
                    789:          S1 = SCALE*S1
                    790:          WR = SCALE*WR
                    791: *
                    792: *        Now check for two consecutive small subdiagonals.
                    793: *
                    794:          DO 120 J = ILAST - 1, IFIRST + 1, -1
                    795:             ISTART = J
                    796:             TEMP = ABS( S1*H( J, J-1 ) )
                    797:             TEMP2 = ABS( S1*H( J, J )-WR*T( J, J ) )
                    798:             TEMPR = MAX( TEMP, TEMP2 )
                    799:             IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
                    800:                TEMP = TEMP / TEMPR
                    801:                TEMP2 = TEMP2 / TEMPR
                    802:             END IF
                    803:             IF( ABS( ( ASCALE*H( J+1, J ) )*TEMP ).LE.( ASCALE*ATOL )*
                    804:      $          TEMP2 )GO TO 130
                    805:   120    CONTINUE
                    806: *
                    807:          ISTART = IFIRST
                    808:   130    CONTINUE
                    809: *
                    810: *        Do an implicit single-shift QZ sweep.
                    811: *
                    812: *        Initial Q
                    813: *
                    814:          TEMP = S1*H( ISTART, ISTART ) - WR*T( ISTART, ISTART )
                    815:          TEMP2 = S1*H( ISTART+1, ISTART )
                    816:          CALL DLARTG( TEMP, TEMP2, C, S, TEMPR )
                    817: *
                    818: *        Sweep
                    819: *
                    820:          DO 190 J = ISTART, ILAST - 1
                    821:             IF( J.GT.ISTART ) THEN
                    822:                TEMP = H( J, J-1 )
                    823:                CALL DLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
                    824:                H( J+1, J-1 ) = ZERO
                    825:             END IF
                    826: *
                    827:             DO 140 JC = J, ILASTM
                    828:                TEMP = C*H( J, JC ) + S*H( J+1, JC )
                    829:                H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
                    830:                H( J, JC ) = TEMP
                    831:                TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
                    832:                T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
                    833:                T( J, JC ) = TEMP2
                    834:   140       CONTINUE
                    835:             IF( ILQ ) THEN
                    836:                DO 150 JR = 1, N
                    837:                   TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
                    838:                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
                    839:                   Q( JR, J ) = TEMP
                    840:   150          CONTINUE
                    841:             END IF
                    842: *
                    843:             TEMP = T( J+1, J+1 )
                    844:             CALL DLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
                    845:             T( J+1, J ) = ZERO
                    846: *
                    847:             DO 160 JR = IFRSTM, MIN( J+2, ILAST )
                    848:                TEMP = C*H( JR, J+1 ) + S*H( JR, J )
                    849:                H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
                    850:                H( JR, J+1 ) = TEMP
                    851:   160       CONTINUE
                    852:             DO 170 JR = IFRSTM, J
                    853:                TEMP = C*T( JR, J+1 ) + S*T( JR, J )
                    854:                T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
                    855:                T( JR, J+1 ) = TEMP
                    856:   170       CONTINUE
                    857:             IF( ILZ ) THEN
                    858:                DO 180 JR = 1, N
                    859:                   TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
                    860:                   Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
                    861:                   Z( JR, J+1 ) = TEMP
                    862:   180          CONTINUE
                    863:             END IF
                    864:   190    CONTINUE
                    865: *
                    866:          GO TO 350
                    867: *
                    868: *        Use Francis double-shift
                    869: *
                    870: *        Note: the Francis double-shift should work with real shifts,
                    871: *              but only if the block is at least 3x3.
                    872: *              This code may break if this point is reached with
                    873: *              a 2x2 block with real eigenvalues.
                    874: *
                    875:   200    CONTINUE
                    876:          IF( IFIRST+1.EQ.ILAST ) THEN
                    877: *
                    878: *           Special case -- 2x2 block with complex eigenvectors
                    879: *
                    880: *           Step 1: Standardize, that is, rotate so that
                    881: *
                    882: *                       ( B11  0  )
                    883: *                   B = (         )  with B11 non-negative.
                    884: *                       (  0  B22 )
                    885: *
                    886:             CALL DLASV2( T( ILAST-1, ILAST-1 ), T( ILAST-1, ILAST ),
                    887:      $                   T( ILAST, ILAST ), B22, B11, SR, CR, SL, CL )
                    888: *
                    889:             IF( B11.LT.ZERO ) THEN
                    890:                CR = -CR
                    891:                SR = -SR
                    892:                B11 = -B11
                    893:                B22 = -B22
                    894:             END IF
                    895: *
                    896:             CALL DROT( ILASTM+1-IFIRST, H( ILAST-1, ILAST-1 ), LDH,
                    897:      $                 H( ILAST, ILAST-1 ), LDH, CL, SL )
                    898:             CALL DROT( ILAST+1-IFRSTM, H( IFRSTM, ILAST-1 ), 1,
                    899:      $                 H( IFRSTM, ILAST ), 1, CR, SR )
                    900: *
                    901:             IF( ILAST.LT.ILASTM )
                    902:      $         CALL DROT( ILASTM-ILAST, T( ILAST-1, ILAST+1 ), LDT,
                    903:      $                    T( ILAST, ILAST+1 ), LDT, CL, SL )
                    904:             IF( IFRSTM.LT.ILAST-1 )
                    905:      $         CALL DROT( IFIRST-IFRSTM, T( IFRSTM, ILAST-1 ), 1,
                    906:      $                    T( IFRSTM, ILAST ), 1, CR, SR )
                    907: *
                    908:             IF( ILQ )
                    909:      $         CALL DROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL,
                    910:      $                    SL )
                    911:             IF( ILZ )
                    912:      $         CALL DROT( N, Z( 1, ILAST-1 ), 1, Z( 1, ILAST ), 1, CR,
                    913:      $                    SR )
                    914: *
                    915:             T( ILAST-1, ILAST-1 ) = B11
                    916:             T( ILAST-1, ILAST ) = ZERO
                    917:             T( ILAST, ILAST-1 ) = ZERO
                    918:             T( ILAST, ILAST ) = B22
                    919: *
                    920: *           If B22 is negative, negate column ILAST
                    921: *
                    922:             IF( B22.LT.ZERO ) THEN
                    923:                DO 210 J = IFRSTM, ILAST
                    924:                   H( J, ILAST ) = -H( J, ILAST )
                    925:                   T( J, ILAST ) = -T( J, ILAST )
                    926:   210          CONTINUE
                    927: *
                    928:                IF( ILZ ) THEN
                    929:                   DO 220 J = 1, N
                    930:                      Z( J, ILAST ) = -Z( J, ILAST )
                    931:   220             CONTINUE
                    932:                END IF
1.11      bertrand  933:                B22 = -B22
1.1       bertrand  934:             END IF
                    935: *
                    936: *           Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
                    937: *
                    938: *           Recompute shift
                    939: *
                    940:             CALL DLAG2( H( ILAST-1, ILAST-1 ), LDH,
                    941:      $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
                    942:      $                  TEMP, WR, TEMP2, WI )
                    943: *
                    944: *           If standardization has perturbed the shift onto real line,
                    945: *           do another (real single-shift) QR step.
                    946: *
                    947:             IF( WI.EQ.ZERO )
                    948:      $         GO TO 350
                    949:             S1INV = ONE / S1
                    950: *
                    951: *           Do EISPACK (QZVAL) computation of alpha and beta
                    952: *
                    953:             A11 = H( ILAST-1, ILAST-1 )
                    954:             A21 = H( ILAST, ILAST-1 )
                    955:             A12 = H( ILAST-1, ILAST )
                    956:             A22 = H( ILAST, ILAST )
                    957: *
                    958: *           Compute complex Givens rotation on right
                    959: *           (Assume some element of C = (sA - wB) > unfl )
                    960: *                            __
                    961: *           (sA - wB) ( CZ   -SZ )
                    962: *                     ( SZ    CZ )
                    963: *
                    964:             C11R = S1*A11 - WR*B11
                    965:             C11I = -WI*B11
                    966:             C12 = S1*A12
                    967:             C21 = S1*A21
                    968:             C22R = S1*A22 - WR*B22
                    969:             C22I = -WI*B22
                    970: *
                    971:             IF( ABS( C11R )+ABS( C11I )+ABS( C12 ).GT.ABS( C21 )+
                    972:      $          ABS( C22R )+ABS( C22I ) ) THEN
                    973:                T1 = DLAPY3( C12, C11R, C11I )
                    974:                CZ = C12 / T1
                    975:                SZR = -C11R / T1
                    976:                SZI = -C11I / T1
                    977:             ELSE
                    978:                CZ = DLAPY2( C22R, C22I )
                    979:                IF( CZ.LE.SAFMIN ) THEN
                    980:                   CZ = ZERO
                    981:                   SZR = ONE
                    982:                   SZI = ZERO
                    983:                ELSE
                    984:                   TEMPR = C22R / CZ
                    985:                   TEMPI = C22I / CZ
                    986:                   T1 = DLAPY2( CZ, C21 )
                    987:                   CZ = CZ / T1
                    988:                   SZR = -C21*TEMPR / T1
                    989:                   SZI = C21*TEMPI / T1
                    990:                END IF
                    991:             END IF
                    992: *
                    993: *           Compute Givens rotation on left
                    994: *
                    995: *           (  CQ   SQ )
                    996: *           (  __      )  A or B
                    997: *           ( -SQ   CQ )
                    998: *
                    999:             AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 )
                   1000:             BN = ABS( B11 ) + ABS( B22 )
                   1001:             WABS = ABS( WR ) + ABS( WI )
                   1002:             IF( S1*AN.GT.WABS*BN ) THEN
                   1003:                CQ = CZ*B11
                   1004:                SQR = SZR*B22
                   1005:                SQI = -SZI*B22
                   1006:             ELSE
                   1007:                A1R = CZ*A11 + SZR*A12
                   1008:                A1I = SZI*A12
                   1009:                A2R = CZ*A21 + SZR*A22
                   1010:                A2I = SZI*A22
                   1011:                CQ = DLAPY2( A1R, A1I )
                   1012:                IF( CQ.LE.SAFMIN ) THEN
                   1013:                   CQ = ZERO
                   1014:                   SQR = ONE
                   1015:                   SQI = ZERO
                   1016:                ELSE
                   1017:                   TEMPR = A1R / CQ
                   1018:                   TEMPI = A1I / CQ
                   1019:                   SQR = TEMPR*A2R + TEMPI*A2I
                   1020:                   SQI = TEMPI*A2R - TEMPR*A2I
                   1021:                END IF
                   1022:             END IF
                   1023:             T1 = DLAPY3( CQ, SQR, SQI )
                   1024:             CQ = CQ / T1
                   1025:             SQR = SQR / T1
                   1026:             SQI = SQI / T1
                   1027: *
                   1028: *           Compute diagonal elements of QBZ
                   1029: *
                   1030:             TEMPR = SQR*SZR - SQI*SZI
                   1031:             TEMPI = SQR*SZI + SQI*SZR
                   1032:             B1R = CQ*CZ*B11 + TEMPR*B22
                   1033:             B1I = TEMPI*B22
                   1034:             B1A = DLAPY2( B1R, B1I )
                   1035:             B2R = CQ*CZ*B22 + TEMPR*B11
                   1036:             B2I = -TEMPI*B11
                   1037:             B2A = DLAPY2( B2R, B2I )
                   1038: *
                   1039: *           Normalize so beta > 0, and Im( alpha1 ) > 0
                   1040: *
                   1041:             BETA( ILAST-1 ) = B1A
                   1042:             BETA( ILAST ) = B2A
                   1043:             ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV
                   1044:             ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV
                   1045:             ALPHAR( ILAST ) = ( WR*B2A )*S1INV
                   1046:             ALPHAI( ILAST ) = -( WI*B2A )*S1INV
                   1047: *
                   1048: *           Step 3: Go to next block -- exit if finished.
                   1049: *
                   1050:             ILAST = IFIRST - 1
                   1051:             IF( ILAST.LT.ILO )
                   1052:      $         GO TO 380
                   1053: *
                   1054: *           Reset counters
                   1055: *
                   1056:             IITER = 0
                   1057:             ESHIFT = ZERO
                   1058:             IF( .NOT.ILSCHR ) THEN
                   1059:                ILASTM = ILAST
                   1060:                IF( IFRSTM.GT.ILAST )
                   1061:      $            IFRSTM = ILO
                   1062:             END IF
                   1063:             GO TO 350
                   1064:          ELSE
                   1065: *
                   1066: *           Usual case: 3x3 or larger block, using Francis implicit
                   1067: *                       double-shift
                   1068: *
                   1069: *                                    2
                   1070: *           Eigenvalue equation is  w  - c w + d = 0,
                   1071: *
                   1072: *                                         -1 2        -1
                   1073: *           so compute 1st column of  (A B  )  - c A B   + d
                   1074: *           using the formula in QZIT (from EISPACK)
                   1075: *
                   1076: *           We assume that the block is at least 3x3
                   1077: *
                   1078:             AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
                   1079:      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
                   1080:             AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
                   1081:      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
                   1082:             AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
                   1083:      $             ( BSCALE*T( ILAST, ILAST ) )
                   1084:             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
                   1085:      $             ( BSCALE*T( ILAST, ILAST ) )
                   1086:             U12 = T( ILAST-1, ILAST ) / T( ILAST, ILAST )
                   1087:             AD11L = ( ASCALE*H( IFIRST, IFIRST ) ) /
                   1088:      $              ( BSCALE*T( IFIRST, IFIRST ) )
                   1089:             AD21L = ( ASCALE*H( IFIRST+1, IFIRST ) ) /
                   1090:      $              ( BSCALE*T( IFIRST, IFIRST ) )
                   1091:             AD12L = ( ASCALE*H( IFIRST, IFIRST+1 ) ) /
                   1092:      $              ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
                   1093:             AD22L = ( ASCALE*H( IFIRST+1, IFIRST+1 ) ) /
                   1094:      $              ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
                   1095:             AD32L = ( ASCALE*H( IFIRST+2, IFIRST+1 ) ) /
                   1096:      $              ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
                   1097:             U12L = T( IFIRST, IFIRST+1 ) / T( IFIRST+1, IFIRST+1 )
                   1098: *
                   1099:             V( 1 ) = ( AD11-AD11L )*( AD22-AD11L ) - AD12*AD21 +
                   1100:      $               AD21*U12*AD11L + ( AD12L-AD11L*U12L )*AD21L
                   1101:             V( 2 ) = ( ( AD22L-AD11L )-AD21L*U12L-( AD11-AD11L )-
                   1102:      $               ( AD22-AD11L )+AD21*U12 )*AD21L
                   1103:             V( 3 ) = AD32L*AD21L
                   1104: *
                   1105:             ISTART = IFIRST
                   1106: *
                   1107:             CALL DLARFG( 3, V( 1 ), V( 2 ), 1, TAU )
                   1108:             V( 1 ) = ONE
                   1109: *
                   1110: *           Sweep
                   1111: *
                   1112:             DO 290 J = ISTART, ILAST - 2
                   1113: *
                   1114: *              All but last elements: use 3x3 Householder transforms.
                   1115: *
                   1116: *              Zero (j-1)st column of A
                   1117: *
                   1118:                IF( J.GT.ISTART ) THEN
                   1119:                   V( 1 ) = H( J, J-1 )
                   1120:                   V( 2 ) = H( J+1, J-1 )
                   1121:                   V( 3 ) = H( J+2, J-1 )
                   1122: *
                   1123:                   CALL DLARFG( 3, H( J, J-1 ), V( 2 ), 1, TAU )
                   1124:                   V( 1 ) = ONE
                   1125:                   H( J+1, J-1 ) = ZERO
                   1126:                   H( J+2, J-1 ) = ZERO
                   1127:                END IF
                   1128: *
                   1129:                DO 230 JC = J, ILASTM
                   1130:                   TEMP = TAU*( H( J, JC )+V( 2 )*H( J+1, JC )+V( 3 )*
                   1131:      $                   H( J+2, JC ) )
                   1132:                   H( J, JC ) = H( J, JC ) - TEMP
                   1133:                   H( J+1, JC ) = H( J+1, JC ) - TEMP*V( 2 )
                   1134:                   H( J+2, JC ) = H( J+2, JC ) - TEMP*V( 3 )
                   1135:                   TEMP2 = TAU*( T( J, JC )+V( 2 )*T( J+1, JC )+V( 3 )*
                   1136:      $                    T( J+2, JC ) )
                   1137:                   T( J, JC ) = T( J, JC ) - TEMP2
                   1138:                   T( J+1, JC ) = T( J+1, JC ) - TEMP2*V( 2 )
                   1139:                   T( J+2, JC ) = T( J+2, JC ) - TEMP2*V( 3 )
                   1140:   230          CONTINUE
                   1141:                IF( ILQ ) THEN
                   1142:                   DO 240 JR = 1, N
                   1143:                      TEMP = TAU*( Q( JR, J )+V( 2 )*Q( JR, J+1 )+V( 3 )*
                   1144:      $                      Q( JR, J+2 ) )
                   1145:                      Q( JR, J ) = Q( JR, J ) - TEMP
                   1146:                      Q( JR, J+1 ) = Q( JR, J+1 ) - TEMP*V( 2 )
                   1147:                      Q( JR, J+2 ) = Q( JR, J+2 ) - TEMP*V( 3 )
                   1148:   240             CONTINUE
                   1149:                END IF
                   1150: *
                   1151: *              Zero j-th column of B (see DLAGBC for details)
                   1152: *
                   1153: *              Swap rows to pivot
                   1154: *
                   1155:                ILPIVT = .FALSE.
                   1156:                TEMP = MAX( ABS( T( J+1, J+1 ) ), ABS( T( J+1, J+2 ) ) )
                   1157:                TEMP2 = MAX( ABS( T( J+2, J+1 ) ), ABS( T( J+2, J+2 ) ) )
                   1158:                IF( MAX( TEMP, TEMP2 ).LT.SAFMIN ) THEN
                   1159:                   SCALE = ZERO
                   1160:                   U1 = ONE
                   1161:                   U2 = ZERO
                   1162:                   GO TO 250
                   1163:                ELSE IF( TEMP.GE.TEMP2 ) THEN
                   1164:                   W11 = T( J+1, J+1 )
                   1165:                   W21 = T( J+2, J+1 )
                   1166:                   W12 = T( J+1, J+2 )
                   1167:                   W22 = T( J+2, J+2 )
                   1168:                   U1 = T( J+1, J )
                   1169:                   U2 = T( J+2, J )
                   1170:                ELSE
                   1171:                   W21 = T( J+1, J+1 )
                   1172:                   W11 = T( J+2, J+1 )
                   1173:                   W22 = T( J+1, J+2 )
                   1174:                   W12 = T( J+2, J+2 )
                   1175:                   U2 = T( J+1, J )
                   1176:                   U1 = T( J+2, J )
                   1177:                END IF
                   1178: *
                   1179: *              Swap columns if nec.
                   1180: *
                   1181:                IF( ABS( W12 ).GT.ABS( W11 ) ) THEN
                   1182:                   ILPIVT = .TRUE.
                   1183:                   TEMP = W12
                   1184:                   TEMP2 = W22
                   1185:                   W12 = W11
                   1186:                   W22 = W21
                   1187:                   W11 = TEMP
                   1188:                   W21 = TEMP2
                   1189:                END IF
                   1190: *
                   1191: *              LU-factor
                   1192: *
                   1193:                TEMP = W21 / W11
                   1194:                U2 = U2 - TEMP*U1
                   1195:                W22 = W22 - TEMP*W12
                   1196:                W21 = ZERO
                   1197: *
                   1198: *              Compute SCALE
                   1199: *
                   1200:                SCALE = ONE
                   1201:                IF( ABS( W22 ).LT.SAFMIN ) THEN
                   1202:                   SCALE = ZERO
                   1203:                   U2 = ONE
                   1204:                   U1 = -W12 / W11
                   1205:                   GO TO 250
                   1206:                END IF
                   1207:                IF( ABS( W22 ).LT.ABS( U2 ) )
                   1208:      $            SCALE = ABS( W22 / U2 )
                   1209:                IF( ABS( W11 ).LT.ABS( U1 ) )
                   1210:      $            SCALE = MIN( SCALE, ABS( W11 / U1 ) )
                   1211: *
                   1212: *              Solve
                   1213: *
                   1214:                U2 = ( SCALE*U2 ) / W22
                   1215:                U1 = ( SCALE*U1-W12*U2 ) / W11
                   1216: *
                   1217:   250          CONTINUE
                   1218:                IF( ILPIVT ) THEN
                   1219:                   TEMP = U2
                   1220:                   U2 = U1
                   1221:                   U1 = TEMP
                   1222:                END IF
                   1223: *
                   1224: *              Compute Householder Vector
                   1225: *
                   1226:                T1 = SQRT( SCALE**2+U1**2+U2**2 )
                   1227:                TAU = ONE + SCALE / T1
                   1228:                VS = -ONE / ( SCALE+T1 )
                   1229:                V( 1 ) = ONE
                   1230:                V( 2 ) = VS*U1
                   1231:                V( 3 ) = VS*U2
                   1232: *
                   1233: *              Apply transformations from the right.
                   1234: *
                   1235:                DO 260 JR = IFRSTM, MIN( J+3, ILAST )
                   1236:                   TEMP = TAU*( H( JR, J )+V( 2 )*H( JR, J+1 )+V( 3 )*
                   1237:      $                   H( JR, J+2 ) )
                   1238:                   H( JR, J ) = H( JR, J ) - TEMP
                   1239:                   H( JR, J+1 ) = H( JR, J+1 ) - TEMP*V( 2 )
                   1240:                   H( JR, J+2 ) = H( JR, J+2 ) - TEMP*V( 3 )
                   1241:   260          CONTINUE
                   1242:                DO 270 JR = IFRSTM, J + 2
                   1243:                   TEMP = TAU*( T( JR, J )+V( 2 )*T( JR, J+1 )+V( 3 )*
                   1244:      $                   T( JR, J+2 ) )
                   1245:                   T( JR, J ) = T( JR, J ) - TEMP
                   1246:                   T( JR, J+1 ) = T( JR, J+1 ) - TEMP*V( 2 )
                   1247:                   T( JR, J+2 ) = T( JR, J+2 ) - TEMP*V( 3 )
                   1248:   270          CONTINUE
                   1249:                IF( ILZ ) THEN
                   1250:                   DO 280 JR = 1, N
                   1251:                      TEMP = TAU*( Z( JR, J )+V( 2 )*Z( JR, J+1 )+V( 3 )*
                   1252:      $                      Z( JR, J+2 ) )
                   1253:                      Z( JR, J ) = Z( JR, J ) - TEMP
                   1254:                      Z( JR, J+1 ) = Z( JR, J+1 ) - TEMP*V( 2 )
                   1255:                      Z( JR, J+2 ) = Z( JR, J+2 ) - TEMP*V( 3 )
                   1256:   280             CONTINUE
                   1257:                END IF
                   1258:                T( J+1, J ) = ZERO
                   1259:                T( J+2, J ) = ZERO
                   1260:   290       CONTINUE
                   1261: *
                   1262: *           Last elements: Use Givens rotations
                   1263: *
                   1264: *           Rotations from the left
                   1265: *
                   1266:             J = ILAST - 1
                   1267:             TEMP = H( J, J-1 )
                   1268:             CALL DLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
                   1269:             H( J+1, J-1 ) = ZERO
                   1270: *
                   1271:             DO 300 JC = J, ILASTM
                   1272:                TEMP = C*H( J, JC ) + S*H( J+1, JC )
                   1273:                H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
                   1274:                H( J, JC ) = TEMP
                   1275:                TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
                   1276:                T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
                   1277:                T( J, JC ) = TEMP2
                   1278:   300       CONTINUE
                   1279:             IF( ILQ ) THEN
                   1280:                DO 310 JR = 1, N
                   1281:                   TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
                   1282:                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
                   1283:                   Q( JR, J ) = TEMP
                   1284:   310          CONTINUE
                   1285:             END IF
                   1286: *
                   1287: *           Rotations from the right.
                   1288: *
                   1289:             TEMP = T( J+1, J+1 )
                   1290:             CALL DLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
                   1291:             T( J+1, J ) = ZERO
                   1292: *
                   1293:             DO 320 JR = IFRSTM, ILAST
                   1294:                TEMP = C*H( JR, J+1 ) + S*H( JR, J )
                   1295:                H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
                   1296:                H( JR, J+1 ) = TEMP
                   1297:   320       CONTINUE
                   1298:             DO 330 JR = IFRSTM, ILAST - 1
                   1299:                TEMP = C*T( JR, J+1 ) + S*T( JR, J )
                   1300:                T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
                   1301:                T( JR, J+1 ) = TEMP
                   1302:   330       CONTINUE
                   1303:             IF( ILZ ) THEN
                   1304:                DO 340 JR = 1, N
                   1305:                   TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
                   1306:                   Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
                   1307:                   Z( JR, J+1 ) = TEMP
                   1308:   340          CONTINUE
                   1309:             END IF
                   1310: *
                   1311: *           End of Double-Shift code
                   1312: *
                   1313:          END IF
                   1314: *
                   1315:          GO TO 350
                   1316: *
                   1317: *        End of iteration loop
                   1318: *
                   1319:   350    CONTINUE
                   1320:   360 CONTINUE
                   1321: *
                   1322: *     Drop-through = non-convergence
                   1323: *
                   1324:       INFO = ILAST
                   1325:       GO TO 420
                   1326: *
                   1327: *     Successful completion of all QZ steps
                   1328: *
                   1329:   380 CONTINUE
                   1330: *
                   1331: *     Set Eigenvalues 1:ILO-1
                   1332: *
                   1333:       DO 410 J = 1, ILO - 1
                   1334:          IF( T( J, J ).LT.ZERO ) THEN
                   1335:             IF( ILSCHR ) THEN
                   1336:                DO 390 JR = 1, J
                   1337:                   H( JR, J ) = -H( JR, J )
                   1338:                   T( JR, J ) = -T( JR, J )
                   1339:   390          CONTINUE
                   1340:             ELSE
                   1341:                H( J, J ) = -H( J, J )
                   1342:                T( J, J ) = -T( J, J )
                   1343:             END IF
                   1344:             IF( ILZ ) THEN
                   1345:                DO 400 JR = 1, N
                   1346:                   Z( JR, J ) = -Z( JR, J )
                   1347:   400          CONTINUE
                   1348:             END IF
                   1349:          END IF
                   1350:          ALPHAR( J ) = H( J, J )
                   1351:          ALPHAI( J ) = ZERO
                   1352:          BETA( J ) = T( J, J )
                   1353:   410 CONTINUE
                   1354: *
                   1355: *     Normal Termination
                   1356: *
                   1357:       INFO = 0
                   1358: *
                   1359: *     Exit (other than argument error) -- return optimal workspace size
                   1360: *
                   1361:   420 CONTINUE
                   1362:       WORK( 1 ) = DBLE( N )
                   1363:       RETURN
                   1364: *
                   1365: *     End of DHGEQZ
                   1366: *
                   1367:       END

CVSweb interface <joel.bertrand@systella.fr>