SRC/shgeqz.f

   1 *> \brief \b SHGEQZ
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SHGEQZ + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/shgeqz.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/shgeqz.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/shgeqz.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
  22 *                          ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
  23 *                          LWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          COMPQ, COMPZ, JOB
  27 *       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       REAL               ALPHAI( * ), ALPHAR( * ), BETA( * ),
  31 *      $                   H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
  32 *      $                   WORK( * ), Z( LDZ, * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> SHGEQZ 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 SGGHRD.
  50 *>
  51 *> If JOB='S', then the Hessenberg-triangular pair (H,T) is
  52 *> also reduced to generalized Schur form,
  53 *>
  54 *>    H = Q*S*Z**T,  T = Q*P*Z**T,
  55 *>
  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 SGGHRD 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.
  78 *>
  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
  89 *> form:
  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;
 104 *>          = 'S': Compute eigenvalues and the Schur form.
 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL array, dimension (LDQ, N)
 214 *>          On entry, if COMPQ = 'V', the orthogonal matrix Q1 used in
 215 *>          the reduction of (A,B) to generalized Hessenberg form.
 216 *>          On exit, if COMPQ = 'I', the orthogonal matrix of left Schur
 217 *>          vectors of (H,T), and if COMPQ = 'V', the orthogonal matrix
 218 *>          of left Schur vectors of (A,B).
 219 *>          Not referenced if COMPQ = 'N'.
 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 REAL 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 REAL 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 *
 280 *> \author Univ. of Tennessee
 281 *> \author Univ. of California Berkeley
 282 *> \author Univ. of Colorado Denver
 283 *> \author NAG Ltd.
 284 *
 285 *> \date June 2016
 286 *
 287 *> \ingroup realGEcomputational
 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 *  =====================================================================
 303       SUBROUTINE SHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
 304      $                   ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
 305      $                   LWORK, INFO )
 306 *
 307 *  -- LAPACK computational routine (version 3.6.1) --
 308 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 309 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 310 *     June 2016
 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       REAL               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       REAL               HALF, ZERO, ONE, SAFETY
 327       PARAMETER          ( HALF = 0.5E+0, ZERO = 0.0E+0, ONE = 1.0E+0,
 328      $                   SAFETY = 1.0E+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       REAL               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       REAL               V( 3 )
 349 *     ..
 350 *     .. External Functions ..
 351       LOGICAL            LSAME
 352       REAL               SLAMCH, SLANHS, SLAPY2, SLAPY3
 353       EXTERNAL           LSAME, SLAMCH, SLANHS, SLAPY2, SLAPY3
 354 *     ..
 355 *     .. External Subroutines ..
 356       EXTERNAL           SLAG2, SLARFG, SLARTG, SLASET, SLASV2, SROT,
 357      $                   XERBLA
 358 *     ..
 359 *     .. Intrinsic Functions ..
 360       INTRINSIC          ABS, MAX, MIN, REAL, 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( 'SHGEQZ', -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 ) = REAL( 1 )
 441          RETURN
 442       END IF
 443 *
 444 *     Initialize Q and Z
 445 *
 446       IF( ICOMPQ.EQ.3 )
 447      $   CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
 448       IF( ICOMPZ.EQ.3 )
 449      $   CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
 450 *
 451 *     Machine Constants
 452 *
 453       IN = IHI + 1 - ILO
 454       SAFMIN = SLAMCH( 'S' )
 455       SAFMAX = ONE / SAFMIN
 456       ULP = SLAMCH( 'E' )*SLAMCH( 'B' )
 457       ANORM = SLANHS( 'F', IN, H( ILO, ILO ), LDH, WORK )
 458       BNORM = SLANHS( '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 SLARTG( TEMP, H( JCH+1, JCH ), C, S,
 592      $                            H( JCH, JCH ) )
 593                      H( JCH+1, JCH ) = ZERO
 594                      CALL SROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
 595      $                          H( JCH+1, JCH+1 ), LDH, C, S )
 596                      CALL SROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
 597      $                          T( JCH+1, JCH+1 ), LDT, C, S )
 598                      IF( ILQ )
 599      $                  CALL SROT( 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 SLARTG( 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 SROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
 627      $                             T( JCH+1, JCH+2 ), LDT, C, S )
 628                      CALL SROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
 629      $                          H( JCH+1, JCH-1 ), LDH, C, S )
 630                      IF( ILQ )
 631      $                  CALL SROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
 632      $                             C, S )
 633                      TEMP = H( JCH+1, JCH )
 634                      CALL SLARTG( TEMP, H( JCH+1, JCH-1 ), C, S,
 635      $                            H( JCH+1, JCH ) )
 636                      H( JCH+1, JCH-1 ) = ZERO
 637                      CALL SROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
 638      $                          H( IFRSTM, JCH-1 ), 1, C, S )
 639                      CALL SROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
 640      $                          T( IFRSTM, JCH-1 ), 1, C, S )
 641                      IF( ILZ )
 642      $                  CALL SROT( 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 SLARTG( TEMP, H( ILAST, ILAST-1 ), C, S,
 670      $                H( ILAST, ILAST ) )
 671          H( ILAST, ILAST-1 ) = ZERO
 672          CALL SROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
 673      $              H( IFRSTM, ILAST-1 ), 1, C, S )
 674          CALL SROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
 675      $              T( IFRSTM, ILAST-1 ), 1, C, S )
 676          IF( ILZ )
 677      $      CALL SROT( 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 *
 742             IF( ( REAL( MAXIT )*SAFMIN )*ABS( H( ILAST, ILAST-1 ) ).LT.
 743      $          ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
 744                ESHIFT = H( ILAST, ILAST-1 ) /
 745      $                  T( ILAST-1, ILAST-1 )
 746             ELSE
 747                ESHIFT = ESHIFT + ONE / ( SAFMIN*REAL( 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 SLAG2 is the Wilkinson shift (AEP p.512),
 757 *
 758             CALL SLAG2( H( ILAST-1, ILAST-1 ), LDH,
 759      $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
 760      $                  S2, WR, WR2, WI )
 761 *
 762             IF ( ABS( (WR/S1)*T( ILAST, ILAST ) - H( ILAST, ILAST ) )
 763      $         .GT. ABS( (WR2/S2)*T( ILAST, ILAST )
 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
 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 SLARTG( 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 SLARTG( 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 SLARTG( 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 SLASV2( 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 SROT( ILASTM+1-IFIRST, H( ILAST-1, ILAST-1 ), LDH,
 897      $                 H( ILAST, ILAST-1 ), LDH, CL, SL )
 898             CALL SROT( ILAST+1-IFRSTM, H( IFRSTM, ILAST-1 ), 1,
 899      $                 H( IFRSTM, ILAST ), 1, CR, SR )
 900 *
 901             IF( ILAST.LT.ILASTM )
 902      $         CALL SROT( ILASTM-ILAST, T( ILAST-1, ILAST+1 ), LDT,
 903      $                    T( ILAST, ILAST+1 ), LDT, CL, SL )
 904             IF( IFRSTM.LT.ILAST-1 )
 905      $         CALL SROT( IFIRST-IFRSTM, T( IFRSTM, ILAST-1 ), 1,
 906      $                    T( IFRSTM, ILAST ), 1, CR, SR )
 907 *
 908             IF( ILQ )
 909      $         CALL SROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL,
 910      $                    SL )
 911             IF( ILZ )
 912      $         CALL SROT( 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
 933                B22 = -B22
 934             END IF
 935 *
 936 *           Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
 937 *
 938 *           Recompute shift
 939 *
 940             CALL SLAG2( 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 = SLAPY3( C12, C11R, C11I )
 974                CZ = C12 / T1
 975                SZR = -C11R / T1
 976                SZI = -C11I / T1
 977             ELSE
 978                CZ = SLAPY2( 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 = SLAPY2( 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 = SLAPY2( 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 = SLAPY3( 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 = SLAPY2( B1R, B1I )
1035             B2R = CQ*CZ*B22 + TEMPR*B11
1036             B2I = -TEMPI*B11
1037             B2A = SLAPY2( 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 SLARFG( 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 SLARFG( 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 SLAGBC 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 SLARTG( 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 SLARTG( 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 ) = REAL( N )
1363       RETURN
1364 *
1365 *     End of SHGEQZ
1366 *
1367       END