LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dhgeqz.f
Go to the documentation of this file.
00001 *> \brief \b DHGEQZ
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DHGEQZ + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dhgeqz.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dhgeqz.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dhgeqz.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
00022 *                          ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
00023 *                          LWORK, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          COMPQ, COMPZ, JOB
00027 *       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       DOUBLE PRECISION   ALPHAI( * ), ALPHAR( * ), BETA( * ),
00031 *      $                   H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
00032 *      $                   WORK( * ), Z( LDZ, * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> DHGEQZ computes the eigenvalues of a real matrix pair (H,T),
00042 *> where H is an upper Hessenberg matrix and T is upper triangular,
00043 *> using the double-shift QZ method.
00044 *> Matrix pairs of this type are produced by the reduction to
00045 *> generalized upper Hessenberg form of a real matrix pair (A,B):
00046 *>
00047 *>    A = Q1*H*Z1**T,  B = Q1*T*Z1**T,
00048 *>
00049 *> as computed by DGGHRD.
00050 *>
00051 *> If JOB='S', then the Hessenberg-triangular pair (H,T) is
00052 *> also reduced to generalized Schur form,
00053 *> 
00054 *>    H = Q*S*Z**T,  T = Q*P*Z**T,
00055 *> 
00056 *> where Q and Z are orthogonal matrices, P is an upper triangular
00057 *> matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
00058 *> diagonal blocks.
00059 *>
00060 *> The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
00061 *> (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
00062 *> eigenvalues.
00063 *>
00064 *> Additionally, the 2-by-2 upper triangular diagonal blocks of P
00065 *> corresponding to 2-by-2 blocks of S are reduced to positive diagonal
00066 *> form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
00067 *> P(j,j) > 0, and P(j+1,j+1) > 0.
00068 *>
00069 *> Optionally, the orthogonal matrix Q from the generalized Schur
00070 *> factorization may be postmultiplied into an input matrix Q1, and the
00071 *> orthogonal matrix Z may be postmultiplied into an input matrix Z1.
00072 *> If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced
00073 *> the matrix pair (A,B) to generalized upper Hessenberg form, then the
00074 *> output matrices Q1*Q and Z1*Z are the orthogonal factors from the
00075 *> generalized Schur factorization of (A,B):
00076 *>
00077 *>    A = (Q1*Q)*S*(Z1*Z)**T,  B = (Q1*Q)*P*(Z1*Z)**T.
00078 *> 
00079 *> To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
00080 *> of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
00081 *> complex and beta real.
00082 *> If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
00083 *> generalized nonsymmetric eigenvalue problem (GNEP)
00084 *>    A*x = lambda*B*x
00085 *> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
00086 *> alternate form of the GNEP
00087 *>    mu*A*y = B*y.
00088 *> Real eigenvalues can be read directly from the generalized Schur
00089 *> form: 
00090 *>   alpha = S(i,i), beta = P(i,i).
00091 *>
00092 *> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
00093 *>      Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
00094 *>      pp. 241--256.
00095 *> \endverbatim
00096 *
00097 *  Arguments:
00098 *  ==========
00099 *
00100 *> \param[in] JOB
00101 *> \verbatim
00102 *>          JOB is CHARACTER*1
00103 *>          = 'E': Compute eigenvalues only;
00104 *>          = 'S': Compute eigenvalues and the Schur form. 
00105 *> \endverbatim
00106 *>
00107 *> \param[in] COMPQ
00108 *> \verbatim
00109 *>          COMPQ is CHARACTER*1
00110 *>          = 'N': Left Schur vectors (Q) are not computed;
00111 *>          = 'I': Q is initialized to the unit matrix and the matrix Q
00112 *>                 of left Schur vectors of (H,T) is returned;
00113 *>          = 'V': Q must contain an orthogonal matrix Q1 on entry and
00114 *>                 the product Q1*Q is returned.
00115 *> \endverbatim
00116 *>
00117 *> \param[in] COMPZ
00118 *> \verbatim
00119 *>          COMPZ is CHARACTER*1
00120 *>          = 'N': Right Schur vectors (Z) are not computed;
00121 *>          = 'I': Z is initialized to the unit matrix and the matrix Z
00122 *>                 of right Schur vectors of (H,T) is returned;
00123 *>          = 'V': Z must contain an orthogonal matrix Z1 on entry and
00124 *>                 the product Z1*Z is returned.
00125 *> \endverbatim
00126 *>
00127 *> \param[in] N
00128 *> \verbatim
00129 *>          N is INTEGER
00130 *>          The order of the matrices H, T, Q, and Z.  N >= 0.
00131 *> \endverbatim
00132 *>
00133 *> \param[in] ILO
00134 *> \verbatim
00135 *>          ILO is INTEGER
00136 *> \endverbatim
00137 *>
00138 *> \param[in] IHI
00139 *> \verbatim
00140 *>          IHI is INTEGER
00141 *>          ILO and IHI mark the rows and columns of H which are in
00142 *>          Hessenberg form.  It is assumed that A is already upper
00143 *>          triangular in rows and columns 1:ILO-1 and IHI+1:N.
00144 *>          If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
00145 *> \endverbatim
00146 *>
00147 *> \param[in,out] H
00148 *> \verbatim
00149 *>          H is DOUBLE PRECISION array, dimension (LDH, N)
00150 *>          On entry, the N-by-N upper Hessenberg matrix H.
00151 *>          On exit, if JOB = 'S', H contains the upper quasi-triangular
00152 *>          matrix S from the generalized Schur factorization.
00153 *>          If JOB = 'E', the diagonal blocks of H match those of S, but
00154 *>          the rest of H is unspecified.
00155 *> \endverbatim
00156 *>
00157 *> \param[in] LDH
00158 *> \verbatim
00159 *>          LDH is INTEGER
00160 *>          The leading dimension of the array H.  LDH >= max( 1, N ).
00161 *> \endverbatim
00162 *>
00163 *> \param[in,out] T
00164 *> \verbatim
00165 *>          T is DOUBLE PRECISION array, dimension (LDT, N)
00166 *>          On entry, the N-by-N upper triangular matrix T.
00167 *>          On exit, if JOB = 'S', T contains the upper triangular
00168 *>          matrix P from the generalized Schur factorization;
00169 *>          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
00170 *>          are reduced to positive diagonal form, i.e., if H(j+1,j) is
00171 *>          non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and
00172 *>          T(j+1,j+1) > 0.
00173 *>          If JOB = 'E', the diagonal blocks of T match those of P, but
00174 *>          the rest of T is unspecified.
00175 *> \endverbatim
00176 *>
00177 *> \param[in] LDT
00178 *> \verbatim
00179 *>          LDT is INTEGER
00180 *>          The leading dimension of the array T.  LDT >= max( 1, N ).
00181 *> \endverbatim
00182 *>
00183 *> \param[out] ALPHAR
00184 *> \verbatim
00185 *>          ALPHAR is DOUBLE PRECISION array, dimension (N)
00186 *>          The real parts of each scalar alpha defining an eigenvalue
00187 *>          of GNEP.
00188 *> \endverbatim
00189 *>
00190 *> \param[out] ALPHAI
00191 *> \verbatim
00192 *>          ALPHAI is DOUBLE PRECISION array, dimension (N)
00193 *>          The imaginary parts of each scalar alpha defining an
00194 *>          eigenvalue of GNEP.
00195 *>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
00196 *>          positive, then the j-th and (j+1)-st eigenvalues are a
00197 *>          complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
00198 *> \endverbatim
00199 *>
00200 *> \param[out] BETA
00201 *> \verbatim
00202 *>          BETA is DOUBLE PRECISION array, dimension (N)
00203 *>          The scalars beta that define the eigenvalues of GNEP.
00204 *>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
00205 *>          beta = BETA(j) represent the j-th eigenvalue of the matrix
00206 *>          pair (A,B), in one of the forms lambda = alpha/beta or
00207 *>          mu = beta/alpha.  Since either lambda or mu may overflow,
00208 *>          they should not, in general, be computed.
00209 *> \endverbatim
00210 *>
00211 *> \param[in,out] Q
00212 *> \verbatim
00213 *>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
00214 *>          On entry, if COMPZ = 'V', the orthogonal matrix Q1 used in
00215 *>          the reduction of (A,B) to generalized Hessenberg form.
00216 *>          On exit, if COMPZ = 'I', the orthogonal matrix of left Schur
00217 *>          vectors of (H,T), and if COMPZ = 'V', the orthogonal matrix
00218 *>          of left Schur vectors of (A,B).
00219 *>          Not referenced if COMPZ = 'N'.
00220 *> \endverbatim
00221 *>
00222 *> \param[in] LDQ
00223 *> \verbatim
00224 *>          LDQ is INTEGER
00225 *>          The leading dimension of the array Q.  LDQ >= 1.
00226 *>          If COMPQ='V' or 'I', then LDQ >= N.
00227 *> \endverbatim
00228 *>
00229 *> \param[in,out] Z
00230 *> \verbatim
00231 *>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
00232 *>          On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
00233 *>          the reduction of (A,B) to generalized Hessenberg form.
00234 *>          On exit, if COMPZ = 'I', the orthogonal matrix of
00235 *>          right Schur vectors of (H,T), and if COMPZ = 'V', the
00236 *>          orthogonal matrix of right Schur vectors of (A,B).
00237 *>          Not referenced if COMPZ = 'N'.
00238 *> \endverbatim
00239 *>
00240 *> \param[in] LDZ
00241 *> \verbatim
00242 *>          LDZ is INTEGER
00243 *>          The leading dimension of the array Z.  LDZ >= 1.
00244 *>          If COMPZ='V' or 'I', then LDZ >= N.
00245 *> \endverbatim
00246 *>
00247 *> \param[out] WORK
00248 *> \verbatim
00249 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00250 *>          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
00251 *> \endverbatim
00252 *>
00253 *> \param[in] LWORK
00254 *> \verbatim
00255 *>          LWORK is INTEGER
00256 *>          The dimension of the array WORK.  LWORK >= max(1,N).
00257 *>
00258 *>          If LWORK = -1, then a workspace query is assumed; the routine
00259 *>          only calculates the optimal size of the WORK array, returns
00260 *>          this value as the first entry of the WORK array, and no error
00261 *>          message related to LWORK is issued by XERBLA.
00262 *> \endverbatim
00263 *>
00264 *> \param[out] INFO
00265 *> \verbatim
00266 *>          INFO is INTEGER
00267 *>          = 0: successful exit
00268 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00269 *>          = 1,...,N: the QZ iteration did not converge.  (H,T) is not
00270 *>                     in Schur form, but ALPHAR(i), ALPHAI(i), and
00271 *>                     BETA(i), i=INFO+1,...,N should be correct.
00272 *>          = N+1,...,2*N: the shift calculation failed.  (H,T) is not
00273 *>                     in Schur form, but ALPHAR(i), ALPHAI(i), and
00274 *>                     BETA(i), i=INFO-N+1,...,N should be correct.
00275 *> \endverbatim
00276 *
00277 *  Authors:
00278 *  ========
00279 *
00280 *> \author Univ. of Tennessee 
00281 *> \author Univ. of California Berkeley 
00282 *> \author Univ. of Colorado Denver 
00283 *> \author NAG Ltd. 
00284 *
00285 *> \date April 2012
00286 *
00287 *> \ingroup doubleGEcomputational
00288 *
00289 *> \par Further Details:
00290 *  =====================
00291 *>
00292 *> \verbatim
00293 *>
00294 *>  Iteration counters:
00295 *>
00296 *>  JITER  -- counts iterations.
00297 *>  IITER  -- counts iterations run since ILAST was last
00298 *>            changed.  This is therefore reset only when a 1-by-1 or
00299 *>            2-by-2 block deflates off the bottom.
00300 *> \endverbatim
00301 *>
00302 *  =====================================================================
00303       SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
00304      $                   ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
00305      $                   LWORK, INFO )
00306 *
00307 *  -- LAPACK computational routine (version 3.4.1) --
00308 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00309 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00310 *     April 2012
00311 *
00312 *     .. Scalar Arguments ..
00313       CHARACTER          COMPQ, COMPZ, JOB
00314       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
00315 *     ..
00316 *     .. Array Arguments ..
00317       DOUBLE PRECISION   ALPHAI( * ), ALPHAR( * ), BETA( * ),
00318      $                   H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
00319      $                   WORK( * ), Z( LDZ, * )
00320 *     ..
00321 *
00322 *  =====================================================================
00323 *
00324 *     .. Parameters ..
00325 *    $                     SAFETY = 1.0E+0 )
00326       DOUBLE PRECISION   HALF, ZERO, ONE, SAFETY
00327       PARAMETER          ( HALF = 0.5D+0, ZERO = 0.0D+0, ONE = 1.0D+0,
00328      $                   SAFETY = 1.0D+2 )
00329 *     ..
00330 *     .. Local Scalars ..
00331       LOGICAL            ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ,
00332      $                   LQUERY
00333       INTEGER            ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
00334      $                   ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
00335      $                   JR, MAXIT
00336       DOUBLE PRECISION   A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
00337      $                   AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
00338      $                   AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I,
00339      $                   B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE,
00340      $                   BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL,
00341      $                   CQ, CR, CZ, ESHIFT, S, S1, S1INV, S2, SAFMAX,
00342      $                   SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, SZR, T1,
00343      $                   TAU, TEMP, TEMP2, TEMPI, TEMPR, U1, U12, U12L,
00344      $                   U2, ULP, VS, W11, W12, W21, W22, WABS, WI, WR,
00345      $                   WR2
00346 *     ..
00347 *     .. Local Arrays ..
00348       DOUBLE PRECISION   V( 3 )
00349 *     ..
00350 *     .. External Functions ..
00351       LOGICAL            LSAME
00352       DOUBLE PRECISION   DLAMCH, DLANHS, DLAPY2, DLAPY3
00353       EXTERNAL           LSAME, DLAMCH, DLANHS, DLAPY2, DLAPY3
00354 *     ..
00355 *     .. External Subroutines ..
00356       EXTERNAL           DLAG2, DLARFG, DLARTG, DLASET, DLASV2, DROT,
00357      $                   XERBLA
00358 *     ..
00359 *     .. Intrinsic Functions ..
00360       INTRINSIC          ABS, DBLE, MAX, MIN, SQRT
00361 *     ..
00362 *     .. Executable Statements ..
00363 *
00364 *     Decode JOB, COMPQ, COMPZ
00365 *
00366       IF( LSAME( JOB, 'E' ) ) THEN
00367          ILSCHR = .FALSE.
00368          ISCHUR = 1
00369       ELSE IF( LSAME( JOB, 'S' ) ) THEN
00370          ILSCHR = .TRUE.
00371          ISCHUR = 2
00372       ELSE
00373          ISCHUR = 0
00374       END IF
00375 *
00376       IF( LSAME( COMPQ, 'N' ) ) THEN
00377          ILQ = .FALSE.
00378          ICOMPQ = 1
00379       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
00380          ILQ = .TRUE.
00381          ICOMPQ = 2
00382       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
00383          ILQ = .TRUE.
00384          ICOMPQ = 3
00385       ELSE
00386          ICOMPQ = 0
00387       END IF
00388 *
00389       IF( LSAME( COMPZ, 'N' ) ) THEN
00390          ILZ = .FALSE.
00391          ICOMPZ = 1
00392       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
00393          ILZ = .TRUE.
00394          ICOMPZ = 2
00395       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
00396          ILZ = .TRUE.
00397          ICOMPZ = 3
00398       ELSE
00399          ICOMPZ = 0
00400       END IF
00401 *
00402 *     Check Argument Values
00403 *
00404       INFO = 0
00405       WORK( 1 ) = MAX( 1, N )
00406       LQUERY = ( LWORK.EQ.-1 )
00407       IF( ISCHUR.EQ.0 ) THEN
00408          INFO = -1
00409       ELSE IF( ICOMPQ.EQ.0 ) THEN
00410          INFO = -2
00411       ELSE IF( ICOMPZ.EQ.0 ) THEN
00412          INFO = -3
00413       ELSE IF( N.LT.0 ) THEN
00414          INFO = -4
00415       ELSE IF( ILO.LT.1 ) THEN
00416          INFO = -5
00417       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
00418          INFO = -6
00419       ELSE IF( LDH.LT.N ) THEN
00420          INFO = -8
00421       ELSE IF( LDT.LT.N ) THEN
00422          INFO = -10
00423       ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
00424          INFO = -15
00425       ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
00426          INFO = -17
00427       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
00428          INFO = -19
00429       END IF
00430       IF( INFO.NE.0 ) THEN
00431          CALL XERBLA( 'DHGEQZ', -INFO )
00432          RETURN
00433       ELSE IF( LQUERY ) THEN
00434          RETURN
00435       END IF
00436 *
00437 *     Quick return if possible
00438 *
00439       IF( N.LE.0 ) THEN
00440          WORK( 1 ) = DBLE( 1 )
00441          RETURN
00442       END IF
00443 *
00444 *     Initialize Q and Z
00445 *
00446       IF( ICOMPQ.EQ.3 )
00447      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
00448       IF( ICOMPZ.EQ.3 )
00449      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
00450 *
00451 *     Machine Constants
00452 *
00453       IN = IHI + 1 - ILO
00454       SAFMIN = DLAMCH( 'S' )
00455       SAFMAX = ONE / SAFMIN
00456       ULP = DLAMCH( 'E' )*DLAMCH( 'B' )
00457       ANORM = DLANHS( 'F', IN, H( ILO, ILO ), LDH, WORK )
00458       BNORM = DLANHS( 'F', IN, T( ILO, ILO ), LDT, WORK )
00459       ATOL = MAX( SAFMIN, ULP*ANORM )
00460       BTOL = MAX( SAFMIN, ULP*BNORM )
00461       ASCALE = ONE / MAX( SAFMIN, ANORM )
00462       BSCALE = ONE / MAX( SAFMIN, BNORM )
00463 *
00464 *     Set Eigenvalues IHI+1:N
00465 *
00466       DO 30 J = IHI + 1, N
00467          IF( T( J, J ).LT.ZERO ) THEN
00468             IF( ILSCHR ) THEN
00469                DO 10 JR = 1, J
00470                   H( JR, J ) = -H( JR, J )
00471                   T( JR, J ) = -T( JR, J )
00472    10          CONTINUE
00473             ELSE
00474                H( J, J ) = -H( J, J )
00475                T( J, J ) = -T( J, J )
00476             END IF
00477             IF( ILZ ) THEN
00478                DO 20 JR = 1, N
00479                   Z( JR, J ) = -Z( JR, J )
00480    20          CONTINUE
00481             END IF
00482          END IF
00483          ALPHAR( J ) = H( J, J )
00484          ALPHAI( J ) = ZERO
00485          BETA( J ) = T( J, J )
00486    30 CONTINUE
00487 *
00488 *     If IHI < ILO, skip QZ steps
00489 *
00490       IF( IHI.LT.ILO )
00491      $   GO TO 380
00492 *
00493 *     MAIN QZ ITERATION LOOP
00494 *
00495 *     Initialize dynamic indices
00496 *
00497 *     Eigenvalues ILAST+1:N have been found.
00498 *        Column operations modify rows IFRSTM:whatever.
00499 *        Row operations modify columns whatever:ILASTM.
00500 *
00501 *     If only eigenvalues are being computed, then
00502 *        IFRSTM is the row of the last splitting row above row ILAST;
00503 *        this is always at least ILO.
00504 *     IITER counts iterations since the last eigenvalue was found,
00505 *        to tell when to use an extraordinary shift.
00506 *     MAXIT is the maximum number of QZ sweeps allowed.
00507 *
00508       ILAST = IHI
00509       IF( ILSCHR ) THEN
00510          IFRSTM = 1
00511          ILASTM = N
00512       ELSE
00513          IFRSTM = ILO
00514          ILASTM = IHI
00515       END IF
00516       IITER = 0
00517       ESHIFT = ZERO
00518       MAXIT = 30*( IHI-ILO+1 )
00519 *
00520       DO 360 JITER = 1, MAXIT
00521 *
00522 *        Split the matrix if possible.
00523 *
00524 *        Two tests:
00525 *           1: H(j,j-1)=0  or  j=ILO
00526 *           2: T(j,j)=0
00527 *
00528          IF( ILAST.EQ.ILO ) THEN
00529 *
00530 *           Special case: j=ILAST
00531 *
00532             GO TO 80
00533          ELSE
00534             IF( ABS( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
00535                H( ILAST, ILAST-1 ) = ZERO
00536                GO TO 80
00537             END IF
00538          END IF
00539 *
00540          IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
00541             T( ILAST, ILAST ) = ZERO
00542             GO TO 70
00543          END IF
00544 *
00545 *        General case: j<ILAST
00546 *
00547          DO 60 J = ILAST - 1, ILO, -1
00548 *
00549 *           Test 1: for H(j,j-1)=0 or j=ILO
00550 *
00551             IF( J.EQ.ILO ) THEN
00552                ILAZRO = .TRUE.
00553             ELSE
00554                IF( ABS( H( J, J-1 ) ).LE.ATOL ) THEN
00555                   H( J, J-1 ) = ZERO
00556                   ILAZRO = .TRUE.
00557                ELSE
00558                   ILAZRO = .FALSE.
00559                END IF
00560             END IF
00561 *
00562 *           Test 2: for T(j,j)=0
00563 *
00564             IF( ABS( T( J, J ) ).LT.BTOL ) THEN
00565                T( J, J ) = ZERO
00566 *
00567 *              Test 1a: Check for 2 consecutive small subdiagonals in A
00568 *
00569                ILAZR2 = .FALSE.
00570                IF( .NOT.ILAZRO ) THEN
00571                   TEMP = ABS( H( J, J-1 ) )
00572                   TEMP2 = ABS( H( J, J ) )
00573                   TEMPR = MAX( TEMP, TEMP2 )
00574                   IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
00575                      TEMP = TEMP / TEMPR
00576                      TEMP2 = TEMP2 / TEMPR
00577                   END IF
00578                   IF( TEMP*( ASCALE*ABS( H( J+1, J ) ) ).LE.TEMP2*
00579      $                ( ASCALE*ATOL ) )ILAZR2 = .TRUE.
00580                END IF
00581 *
00582 *              If both tests pass (1 & 2), i.e., the leading diagonal
00583 *              element of B in the block is zero, split a 1x1 block off
00584 *              at the top. (I.e., at the J-th row/column) The leading
00585 *              diagonal element of the remainder can also be zero, so
00586 *              this may have to be done repeatedly.
00587 *
00588                IF( ILAZRO .OR. ILAZR2 ) THEN
00589                   DO 40 JCH = J, ILAST - 1
00590                      TEMP = H( JCH, JCH )
00591                      CALL DLARTG( TEMP, H( JCH+1, JCH ), C, S,
00592      $                            H( JCH, JCH ) )
00593                      H( JCH+1, JCH ) = ZERO
00594                      CALL DROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
00595      $                          H( JCH+1, JCH+1 ), LDH, C, S )
00596                      CALL DROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
00597      $                          T( JCH+1, JCH+1 ), LDT, C, S )
00598                      IF( ILQ )
00599      $                  CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
00600      $                             C, S )
00601                      IF( ILAZR2 )
00602      $                  H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
00603                      ILAZR2 = .FALSE.
00604                      IF( ABS( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
00605                         IF( JCH+1.GE.ILAST ) THEN
00606                            GO TO 80
00607                         ELSE
00608                            IFIRST = JCH + 1
00609                            GO TO 110
00610                         END IF
00611                      END IF
00612                      T( JCH+1, JCH+1 ) = ZERO
00613    40             CONTINUE
00614                   GO TO 70
00615                ELSE
00616 *
00617 *                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
00618 *                 Then process as in the case T(ILAST,ILAST)=0
00619 *
00620                   DO 50 JCH = J, ILAST - 1
00621                      TEMP = T( JCH, JCH+1 )
00622                      CALL DLARTG( TEMP, T( JCH+1, JCH+1 ), C, S,
00623      $                            T( JCH, JCH+1 ) )
00624                      T( JCH+1, JCH+1 ) = ZERO
00625                      IF( JCH.LT.ILASTM-1 )
00626      $                  CALL DROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
00627      $                             T( JCH+1, JCH+2 ), LDT, C, S )
00628                      CALL DROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
00629      $                          H( JCH+1, JCH-1 ), LDH, C, S )
00630                      IF( ILQ )
00631      $                  CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
00632      $                             C, S )
00633                      TEMP = H( JCH+1, JCH )
00634                      CALL DLARTG( TEMP, H( JCH+1, JCH-1 ), C, S,
00635      $                            H( JCH+1, JCH ) )
00636                      H( JCH+1, JCH-1 ) = ZERO
00637                      CALL DROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
00638      $                          H( IFRSTM, JCH-1 ), 1, C, S )
00639                      CALL DROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
00640      $                          T( IFRSTM, JCH-1 ), 1, C, S )
00641                      IF( ILZ )
00642      $                  CALL DROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
00643      $                             C, S )
00644    50             CONTINUE
00645                   GO TO 70
00646                END IF
00647             ELSE IF( ILAZRO ) THEN
00648 *
00649 *              Only test 1 passed -- work on J:ILAST
00650 *
00651                IFIRST = J
00652                GO TO 110
00653             END IF
00654 *
00655 *           Neither test passed -- try next J
00656 *
00657    60    CONTINUE
00658 *
00659 *        (Drop-through is "impossible")
00660 *
00661          INFO = N + 1
00662          GO TO 420
00663 *
00664 *        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
00665 *        1x1 block.
00666 *
00667    70    CONTINUE
00668          TEMP = H( ILAST, ILAST )
00669          CALL DLARTG( TEMP, H( ILAST, ILAST-1 ), C, S,
00670      $                H( ILAST, ILAST ) )
00671          H( ILAST, ILAST-1 ) = ZERO
00672          CALL DROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
00673      $              H( IFRSTM, ILAST-1 ), 1, C, S )
00674          CALL DROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
00675      $              T( IFRSTM, ILAST-1 ), 1, C, S )
00676          IF( ILZ )
00677      $      CALL DROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
00678 *
00679 *        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
00680 *                              and BETA
00681 *
00682    80    CONTINUE
00683          IF( T( ILAST, ILAST ).LT.ZERO ) THEN
00684             IF( ILSCHR ) THEN
00685                DO 90 J = IFRSTM, ILAST
00686                   H( J, ILAST ) = -H( J, ILAST )
00687                   T( J, ILAST ) = -T( J, ILAST )
00688    90          CONTINUE
00689             ELSE
00690                H( ILAST, ILAST ) = -H( ILAST, ILAST )
00691                T( ILAST, ILAST ) = -T( ILAST, ILAST )
00692             END IF
00693             IF( ILZ ) THEN
00694                DO 100 J = 1, N
00695                   Z( J, ILAST ) = -Z( J, ILAST )
00696   100          CONTINUE
00697             END IF
00698          END IF
00699          ALPHAR( ILAST ) = H( ILAST, ILAST )
00700          ALPHAI( ILAST ) = ZERO
00701          BETA( ILAST ) = T( ILAST, ILAST )
00702 *
00703 *        Go to next block -- exit if finished.
00704 *
00705          ILAST = ILAST - 1
00706          IF( ILAST.LT.ILO )
00707      $      GO TO 380
00708 *
00709 *        Reset counters
00710 *
00711          IITER = 0
00712          ESHIFT = ZERO
00713          IF( .NOT.ILSCHR ) THEN
00714             ILASTM = ILAST
00715             IF( IFRSTM.GT.ILAST )
00716      $         IFRSTM = ILO
00717          END IF
00718          GO TO 350
00719 *
00720 *        QZ step
00721 *
00722 *        This iteration only involves rows/columns IFIRST:ILAST. We
00723 *        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
00724 *
00725   110    CONTINUE
00726          IITER = IITER + 1
00727          IF( .NOT.ILSCHR ) THEN
00728             IFRSTM = IFIRST
00729          END IF
00730 *
00731 *        Compute single shifts.
00732 *
00733 *        At this point, IFIRST < ILAST, and the diagonal elements of
00734 *        T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
00735 *        magnitude)
00736 *
00737          IF( ( IITER / 10 )*10.EQ.IITER ) THEN
00738 *
00739 *           Exceptional shift.  Chosen for no particularly good reason.
00740 *           (Single shift only.)
00741 *
00742             IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST ) ).LT.
00743      $          ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
00744                ESHIFT = ESHIFT + H( ILAST, ILAST-1 ) /
00745      $                  T( ILAST-1, ILAST-1 )
00746             ELSE
00747                ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) )
00748             END IF
00749             S1 = ONE
00750             WR = ESHIFT
00751 *
00752          ELSE
00753 *
00754 *           Shifts based on the generalized eigenvalues of the
00755 *           bottom-right 2x2 block of A and B. The first eigenvalue
00756 *           returned by DLAG2 is the Wilkinson shift (AEP p.512),
00757 *
00758             CALL DLAG2( H( ILAST-1, ILAST-1 ), LDH,
00759      $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
00760      $                  S2, WR, WR2, WI )
00761 *
00762             TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) )
00763             IF( WI.NE.ZERO )
00764      $         GO TO 200
00765          END IF
00766 *
00767 *        Fiddle with shift to avoid overflow
00768 *
00769          TEMP = MIN( ASCALE, ONE )*( HALF*SAFMAX )
00770          IF( S1.GT.TEMP ) THEN
00771             SCALE = TEMP / S1
00772          ELSE
00773             SCALE = ONE
00774          END IF
00775 *
00776          TEMP = MIN( BSCALE, ONE )*( HALF*SAFMAX )
00777          IF( ABS( WR ).GT.TEMP )
00778      $      SCALE = MIN( SCALE, TEMP / ABS( WR ) )
00779          S1 = SCALE*S1
00780          WR = SCALE*WR
00781 *
00782 *        Now check for two consecutive small subdiagonals.
00783 *
00784          DO 120 J = ILAST - 1, IFIRST + 1, -1
00785             ISTART = J
00786             TEMP = ABS( S1*H( J, J-1 ) )
00787             TEMP2 = ABS( S1*H( J, J )-WR*T( J, J ) )
00788             TEMPR = MAX( TEMP, TEMP2 )
00789             IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
00790                TEMP = TEMP / TEMPR
00791                TEMP2 = TEMP2 / TEMPR
00792             END IF
00793             IF( ABS( ( ASCALE*H( J+1, J ) )*TEMP ).LE.( ASCALE*ATOL )*
00794      $          TEMP2 )GO TO 130
00795   120    CONTINUE
00796 *
00797          ISTART = IFIRST
00798   130    CONTINUE
00799 *
00800 *        Do an implicit single-shift QZ sweep.
00801 *
00802 *        Initial Q
00803 *
00804          TEMP = S1*H( ISTART, ISTART ) - WR*T( ISTART, ISTART )
00805          TEMP2 = S1*H( ISTART+1, ISTART )
00806          CALL DLARTG( TEMP, TEMP2, C, S, TEMPR )
00807 *
00808 *        Sweep
00809 *
00810          DO 190 J = ISTART, ILAST - 1
00811             IF( J.GT.ISTART ) THEN
00812                TEMP = H( J, J-1 )
00813                CALL DLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
00814                H( J+1, J-1 ) = ZERO
00815             END IF
00816 *
00817             DO 140 JC = J, ILASTM
00818                TEMP = C*H( J, JC ) + S*H( J+1, JC )
00819                H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
00820                H( J, JC ) = TEMP
00821                TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
00822                T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
00823                T( J, JC ) = TEMP2
00824   140       CONTINUE
00825             IF( ILQ ) THEN
00826                DO 150 JR = 1, N
00827                   TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
00828                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
00829                   Q( JR, J ) = TEMP
00830   150          CONTINUE
00831             END IF
00832 *
00833             TEMP = T( J+1, J+1 )
00834             CALL DLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
00835             T( J+1, J ) = ZERO
00836 *
00837             DO 160 JR = IFRSTM, MIN( J+2, ILAST )
00838                TEMP = C*H( JR, J+1 ) + S*H( JR, J )
00839                H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
00840                H( JR, J+1 ) = TEMP
00841   160       CONTINUE
00842             DO 170 JR = IFRSTM, J
00843                TEMP = C*T( JR, J+1 ) + S*T( JR, J )
00844                T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
00845                T( JR, J+1 ) = TEMP
00846   170       CONTINUE
00847             IF( ILZ ) THEN
00848                DO 180 JR = 1, N
00849                   TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
00850                   Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
00851                   Z( JR, J+1 ) = TEMP
00852   180          CONTINUE
00853             END IF
00854   190    CONTINUE
00855 *
00856          GO TO 350
00857 *
00858 *        Use Francis double-shift
00859 *
00860 *        Note: the Francis double-shift should work with real shifts,
00861 *              but only if the block is at least 3x3.
00862 *              This code may break if this point is reached with
00863 *              a 2x2 block with real eigenvalues.
00864 *
00865   200    CONTINUE
00866          IF( IFIRST+1.EQ.ILAST ) THEN
00867 *
00868 *           Special case -- 2x2 block with complex eigenvectors
00869 *
00870 *           Step 1: Standardize, that is, rotate so that
00871 *
00872 *                       ( B11  0  )
00873 *                   B = (         )  with B11 non-negative.
00874 *                       (  0  B22 )
00875 *
00876             CALL DLASV2( T( ILAST-1, ILAST-1 ), T( ILAST-1, ILAST ),
00877      $                   T( ILAST, ILAST ), B22, B11, SR, CR, SL, CL )
00878 *
00879             IF( B11.LT.ZERO ) THEN
00880                CR = -CR
00881                SR = -SR
00882                B11 = -B11
00883                B22 = -B22
00884             END IF
00885 *
00886             CALL DROT( ILASTM+1-IFIRST, H( ILAST-1, ILAST-1 ), LDH,
00887      $                 H( ILAST, ILAST-1 ), LDH, CL, SL )
00888             CALL DROT( ILAST+1-IFRSTM, H( IFRSTM, ILAST-1 ), 1,
00889      $                 H( IFRSTM, ILAST ), 1, CR, SR )
00890 *
00891             IF( ILAST.LT.ILASTM )
00892      $         CALL DROT( ILASTM-ILAST, T( ILAST-1, ILAST+1 ), LDT,
00893      $                    T( ILAST, ILAST+1 ), LDT, CL, SL )
00894             IF( IFRSTM.LT.ILAST-1 )
00895      $         CALL DROT( IFIRST-IFRSTM, T( IFRSTM, ILAST-1 ), 1,
00896      $                    T( IFRSTM, ILAST ), 1, CR, SR )
00897 *
00898             IF( ILQ )
00899      $         CALL DROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL,
00900      $                    SL )
00901             IF( ILZ )
00902      $         CALL DROT( N, Z( 1, ILAST-1 ), 1, Z( 1, ILAST ), 1, CR,
00903      $                    SR )
00904 *
00905             T( ILAST-1, ILAST-1 ) = B11
00906             T( ILAST-1, ILAST ) = ZERO
00907             T( ILAST, ILAST-1 ) = ZERO
00908             T( ILAST, ILAST ) = B22
00909 *
00910 *           If B22 is negative, negate column ILAST
00911 *
00912             IF( B22.LT.ZERO ) THEN
00913                DO 210 J = IFRSTM, ILAST
00914                   H( J, ILAST ) = -H( J, ILAST )
00915                   T( J, ILAST ) = -T( J, ILAST )
00916   210          CONTINUE
00917 *
00918                IF( ILZ ) THEN
00919                   DO 220 J = 1, N
00920                      Z( J, ILAST ) = -Z( J, ILAST )
00921   220             CONTINUE
00922                END IF
00923                B22 = -B22
00924             END IF
00925 *
00926 *           Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
00927 *
00928 *           Recompute shift
00929 *
00930             CALL DLAG2( H( ILAST-1, ILAST-1 ), LDH,
00931      $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
00932      $                  TEMP, WR, TEMP2, WI )
00933 *
00934 *           If standardization has perturbed the shift onto real line,
00935 *           do another (real single-shift) QR step.
00936 *
00937             IF( WI.EQ.ZERO )
00938      $         GO TO 350
00939             S1INV = ONE / S1
00940 *
00941 *           Do EISPACK (QZVAL) computation of alpha and beta
00942 *
00943             A11 = H( ILAST-1, ILAST-1 )
00944             A21 = H( ILAST, ILAST-1 )
00945             A12 = H( ILAST-1, ILAST )
00946             A22 = H( ILAST, ILAST )
00947 *
00948 *           Compute complex Givens rotation on right
00949 *           (Assume some element of C = (sA - wB) > unfl )
00950 *                            __
00951 *           (sA - wB) ( CZ   -SZ )
00952 *                     ( SZ    CZ )
00953 *
00954             C11R = S1*A11 - WR*B11
00955             C11I = -WI*B11
00956             C12 = S1*A12
00957             C21 = S1*A21
00958             C22R = S1*A22 - WR*B22
00959             C22I = -WI*B22
00960 *
00961             IF( ABS( C11R )+ABS( C11I )+ABS( C12 ).GT.ABS( C21 )+
00962      $          ABS( C22R )+ABS( C22I ) ) THEN
00963                T1 = DLAPY3( C12, C11R, C11I )
00964                CZ = C12 / T1
00965                SZR = -C11R / T1
00966                SZI = -C11I / T1
00967             ELSE
00968                CZ = DLAPY2( C22R, C22I )
00969                IF( CZ.LE.SAFMIN ) THEN
00970                   CZ = ZERO
00971                   SZR = ONE
00972                   SZI = ZERO
00973                ELSE
00974                   TEMPR = C22R / CZ
00975                   TEMPI = C22I / CZ
00976                   T1 = DLAPY2( CZ, C21 )
00977                   CZ = CZ / T1
00978                   SZR = -C21*TEMPR / T1
00979                   SZI = C21*TEMPI / T1
00980                END IF
00981             END IF
00982 *
00983 *           Compute Givens rotation on left
00984 *
00985 *           (  CQ   SQ )
00986 *           (  __      )  A or B
00987 *           ( -SQ   CQ )
00988 *
00989             AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 )
00990             BN = ABS( B11 ) + ABS( B22 )
00991             WABS = ABS( WR ) + ABS( WI )
00992             IF( S1*AN.GT.WABS*BN ) THEN
00993                CQ = CZ*B11
00994                SQR = SZR*B22
00995                SQI = -SZI*B22
00996             ELSE
00997                A1R = CZ*A11 + SZR*A12
00998                A1I = SZI*A12
00999                A2R = CZ*A21 + SZR*A22
01000                A2I = SZI*A22
01001                CQ = DLAPY2( A1R, A1I )
01002                IF( CQ.LE.SAFMIN ) THEN
01003                   CQ = ZERO
01004                   SQR = ONE
01005                   SQI = ZERO
01006                ELSE
01007                   TEMPR = A1R / CQ
01008                   TEMPI = A1I / CQ
01009                   SQR = TEMPR*A2R + TEMPI*A2I
01010                   SQI = TEMPI*A2R - TEMPR*A2I
01011                END IF
01012             END IF
01013             T1 = DLAPY3( CQ, SQR, SQI )
01014             CQ = CQ / T1
01015             SQR = SQR / T1
01016             SQI = SQI / T1
01017 *
01018 *           Compute diagonal elements of QBZ
01019 *
01020             TEMPR = SQR*SZR - SQI*SZI
01021             TEMPI = SQR*SZI + SQI*SZR
01022             B1R = CQ*CZ*B11 + TEMPR*B22
01023             B1I = TEMPI*B22
01024             B1A = DLAPY2( B1R, B1I )
01025             B2R = CQ*CZ*B22 + TEMPR*B11
01026             B2I = -TEMPI*B11
01027             B2A = DLAPY2( B2R, B2I )
01028 *
01029 *           Normalize so beta > 0, and Im( alpha1 ) > 0
01030 *
01031             BETA( ILAST-1 ) = B1A
01032             BETA( ILAST ) = B2A
01033             ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV
01034             ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV
01035             ALPHAR( ILAST ) = ( WR*B2A )*S1INV
01036             ALPHAI( ILAST ) = -( WI*B2A )*S1INV
01037 *
01038 *           Step 3: Go to next block -- exit if finished.
01039 *
01040             ILAST = IFIRST - 1
01041             IF( ILAST.LT.ILO )
01042      $         GO TO 380
01043 *
01044 *           Reset counters
01045 *
01046             IITER = 0
01047             ESHIFT = ZERO
01048             IF( .NOT.ILSCHR ) THEN
01049                ILASTM = ILAST
01050                IF( IFRSTM.GT.ILAST )
01051      $            IFRSTM = ILO
01052             END IF
01053             GO TO 350
01054          ELSE
01055 *
01056 *           Usual case: 3x3 or larger block, using Francis implicit
01057 *                       double-shift
01058 *
01059 *                                    2
01060 *           Eigenvalue equation is  w  - c w + d = 0,
01061 *
01062 *                                         -1 2        -1
01063 *           so compute 1st column of  (A B  )  - c A B   + d
01064 *           using the formula in QZIT (from EISPACK)
01065 *
01066 *           We assume that the block is at least 3x3
01067 *
01068             AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
01069      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
01070             AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
01071      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
01072             AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
01073      $             ( BSCALE*T( ILAST, ILAST ) )
01074             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
01075      $             ( BSCALE*T( ILAST, ILAST ) )
01076             U12 = T( ILAST-1, ILAST ) / T( ILAST, ILAST )
01077             AD11L = ( ASCALE*H( IFIRST, IFIRST ) ) /
01078      $              ( BSCALE*T( IFIRST, IFIRST ) )
01079             AD21L = ( ASCALE*H( IFIRST+1, IFIRST ) ) /
01080      $              ( BSCALE*T( IFIRST, IFIRST ) )
01081             AD12L = ( ASCALE*H( IFIRST, IFIRST+1 ) ) /
01082      $              ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
01083             AD22L = ( ASCALE*H( IFIRST+1, IFIRST+1 ) ) /
01084      $              ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
01085             AD32L = ( ASCALE*H( IFIRST+2, IFIRST+1 ) ) /
01086      $              ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
01087             U12L = T( IFIRST, IFIRST+1 ) / T( IFIRST+1, IFIRST+1 )
01088 *
01089             V( 1 ) = ( AD11-AD11L )*( AD22-AD11L ) - AD12*AD21 +
01090      $               AD21*U12*AD11L + ( AD12L-AD11L*U12L )*AD21L
01091             V( 2 ) = ( ( AD22L-AD11L )-AD21L*U12L-( AD11-AD11L )-
01092      $               ( AD22-AD11L )+AD21*U12 )*AD21L
01093             V( 3 ) = AD32L*AD21L
01094 *
01095             ISTART = IFIRST
01096 *
01097             CALL DLARFG( 3, V( 1 ), V( 2 ), 1, TAU )
01098             V( 1 ) = ONE
01099 *
01100 *           Sweep
01101 *
01102             DO 290 J = ISTART, ILAST - 2
01103 *
01104 *              All but last elements: use 3x3 Householder transforms.
01105 *
01106 *              Zero (j-1)st column of A
01107 *
01108                IF( J.GT.ISTART ) THEN
01109                   V( 1 ) = H( J, J-1 )
01110                   V( 2 ) = H( J+1, J-1 )
01111                   V( 3 ) = H( J+2, J-1 )
01112 *
01113                   CALL DLARFG( 3, H( J, J-1 ), V( 2 ), 1, TAU )
01114                   V( 1 ) = ONE
01115                   H( J+1, J-1 ) = ZERO
01116                   H( J+2, J-1 ) = ZERO
01117                END IF
01118 *
01119                DO 230 JC = J, ILASTM
01120                   TEMP = TAU*( H( J, JC )+V( 2 )*H( J+1, JC )+V( 3 )*
01121      $                   H( J+2, JC ) )
01122                   H( J, JC ) = H( J, JC ) - TEMP
01123                   H( J+1, JC ) = H( J+1, JC ) - TEMP*V( 2 )
01124                   H( J+2, JC ) = H( J+2, JC ) - TEMP*V( 3 )
01125                   TEMP2 = TAU*( T( J, JC )+V( 2 )*T( J+1, JC )+V( 3 )*
01126      $                    T( J+2, JC ) )
01127                   T( J, JC ) = T( J, JC ) - TEMP2
01128                   T( J+1, JC ) = T( J+1, JC ) - TEMP2*V( 2 )
01129                   T( J+2, JC ) = T( J+2, JC ) - TEMP2*V( 3 )
01130   230          CONTINUE
01131                IF( ILQ ) THEN
01132                   DO 240 JR = 1, N
01133                      TEMP = TAU*( Q( JR, J )+V( 2 )*Q( JR, J+1 )+V( 3 )*
01134      $                      Q( JR, J+2 ) )
01135                      Q( JR, J ) = Q( JR, J ) - TEMP
01136                      Q( JR, J+1 ) = Q( JR, J+1 ) - TEMP*V( 2 )
01137                      Q( JR, J+2 ) = Q( JR, J+2 ) - TEMP*V( 3 )
01138   240             CONTINUE
01139                END IF
01140 *
01141 *              Zero j-th column of B (see DLAGBC for details)
01142 *
01143 *              Swap rows to pivot
01144 *
01145                ILPIVT = .FALSE.
01146                TEMP = MAX( ABS( T( J+1, J+1 ) ), ABS( T( J+1, J+2 ) ) )
01147                TEMP2 = MAX( ABS( T( J+2, J+1 ) ), ABS( T( J+2, J+2 ) ) )
01148                IF( MAX( TEMP, TEMP2 ).LT.SAFMIN ) THEN
01149                   SCALE = ZERO
01150                   U1 = ONE
01151                   U2 = ZERO
01152                   GO TO 250
01153                ELSE IF( TEMP.GE.TEMP2 ) THEN
01154                   W11 = T( J+1, J+1 )
01155                   W21 = T( J+2, J+1 )
01156                   W12 = T( J+1, J+2 )
01157                   W22 = T( J+2, J+2 )
01158                   U1 = T( J+1, J )
01159                   U2 = T( J+2, J )
01160                ELSE
01161                   W21 = T( J+1, J+1 )
01162                   W11 = T( J+2, J+1 )
01163                   W22 = T( J+1, J+2 )
01164                   W12 = T( J+2, J+2 )
01165                   U2 = T( J+1, J )
01166                   U1 = T( J+2, J )
01167                END IF
01168 *
01169 *              Swap columns if nec.
01170 *
01171                IF( ABS( W12 ).GT.ABS( W11 ) ) THEN
01172                   ILPIVT = .TRUE.
01173                   TEMP = W12
01174                   TEMP2 = W22
01175                   W12 = W11
01176                   W22 = W21
01177                   W11 = TEMP
01178                   W21 = TEMP2
01179                END IF
01180 *
01181 *              LU-factor
01182 *
01183                TEMP = W21 / W11
01184                U2 = U2 - TEMP*U1
01185                W22 = W22 - TEMP*W12
01186                W21 = ZERO
01187 *
01188 *              Compute SCALE
01189 *
01190                SCALE = ONE
01191                IF( ABS( W22 ).LT.SAFMIN ) THEN
01192                   SCALE = ZERO
01193                   U2 = ONE
01194                   U1 = -W12 / W11
01195                   GO TO 250
01196                END IF
01197                IF( ABS( W22 ).LT.ABS( U2 ) )
01198      $            SCALE = ABS( W22 / U2 )
01199                IF( ABS( W11 ).LT.ABS( U1 ) )
01200      $            SCALE = MIN( SCALE, ABS( W11 / U1 ) )
01201 *
01202 *              Solve
01203 *
01204                U2 = ( SCALE*U2 ) / W22
01205                U1 = ( SCALE*U1-W12*U2 ) / W11
01206 *
01207   250          CONTINUE
01208                IF( ILPIVT ) THEN
01209                   TEMP = U2
01210                   U2 = U1
01211                   U1 = TEMP
01212                END IF
01213 *
01214 *              Compute Householder Vector
01215 *
01216                T1 = SQRT( SCALE**2+U1**2+U2**2 )
01217                TAU = ONE + SCALE / T1
01218                VS = -ONE / ( SCALE+T1 )
01219                V( 1 ) = ONE
01220                V( 2 ) = VS*U1
01221                V( 3 ) = VS*U2
01222 *
01223 *              Apply transformations from the right.
01224 *
01225                DO 260 JR = IFRSTM, MIN( J+3, ILAST )
01226                   TEMP = TAU*( H( JR, J )+V( 2 )*H( JR, J+1 )+V( 3 )*
01227      $                   H( JR, J+2 ) )
01228                   H( JR, J ) = H( JR, J ) - TEMP
01229                   H( JR, J+1 ) = H( JR, J+1 ) - TEMP*V( 2 )
01230                   H( JR, J+2 ) = H( JR, J+2 ) - TEMP*V( 3 )
01231   260          CONTINUE
01232                DO 270 JR = IFRSTM, J + 2
01233                   TEMP = TAU*( T( JR, J )+V( 2 )*T( JR, J+1 )+V( 3 )*
01234      $                   T( JR, J+2 ) )
01235                   T( JR, J ) = T( JR, J ) - TEMP
01236                   T( JR, J+1 ) = T( JR, J+1 ) - TEMP*V( 2 )
01237                   T( JR, J+2 ) = T( JR, J+2 ) - TEMP*V( 3 )
01238   270          CONTINUE
01239                IF( ILZ ) THEN
01240                   DO 280 JR = 1, N
01241                      TEMP = TAU*( Z( JR, J )+V( 2 )*Z( JR, J+1 )+V( 3 )*
01242      $                      Z( JR, J+2 ) )
01243                      Z( JR, J ) = Z( JR, J ) - TEMP
01244                      Z( JR, J+1 ) = Z( JR, J+1 ) - TEMP*V( 2 )
01245                      Z( JR, J+2 ) = Z( JR, J+2 ) - TEMP*V( 3 )
01246   280             CONTINUE
01247                END IF
01248                T( J+1, J ) = ZERO
01249                T( J+2, J ) = ZERO
01250   290       CONTINUE
01251 *
01252 *           Last elements: Use Givens rotations
01253 *
01254 *           Rotations from the left
01255 *
01256             J = ILAST - 1
01257             TEMP = H( J, J-1 )
01258             CALL DLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
01259             H( J+1, J-1 ) = ZERO
01260 *
01261             DO 300 JC = J, ILASTM
01262                TEMP = C*H( J, JC ) + S*H( J+1, JC )
01263                H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
01264                H( J, JC ) = TEMP
01265                TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
01266                T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
01267                T( J, JC ) = TEMP2
01268   300       CONTINUE
01269             IF( ILQ ) THEN
01270                DO 310 JR = 1, N
01271                   TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
01272                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
01273                   Q( JR, J ) = TEMP
01274   310          CONTINUE
01275             END IF
01276 *
01277 *           Rotations from the right.
01278 *
01279             TEMP = T( J+1, J+1 )
01280             CALL DLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
01281             T( J+1, J ) = ZERO
01282 *
01283             DO 320 JR = IFRSTM, ILAST
01284                TEMP = C*H( JR, J+1 ) + S*H( JR, J )
01285                H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
01286                H( JR, J+1 ) = TEMP
01287   320       CONTINUE
01288             DO 330 JR = IFRSTM, ILAST - 1
01289                TEMP = C*T( JR, J+1 ) + S*T( JR, J )
01290                T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
01291                T( JR, J+1 ) = TEMP
01292   330       CONTINUE
01293             IF( ILZ ) THEN
01294                DO 340 JR = 1, N
01295                   TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
01296                   Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
01297                   Z( JR, J+1 ) = TEMP
01298   340          CONTINUE
01299             END IF
01300 *
01301 *           End of Double-Shift code
01302 *
01303          END IF
01304 *
01305          GO TO 350
01306 *
01307 *        End of iteration loop
01308 *
01309   350    CONTINUE
01310   360 CONTINUE
01311 *
01312 *     Drop-through = non-convergence
01313 *
01314       INFO = ILAST
01315       GO TO 420
01316 *
01317 *     Successful completion of all QZ steps
01318 *
01319   380 CONTINUE
01320 *
01321 *     Set Eigenvalues 1:ILO-1
01322 *
01323       DO 410 J = 1, ILO - 1
01324          IF( T( J, J ).LT.ZERO ) THEN
01325             IF( ILSCHR ) THEN
01326                DO 390 JR = 1, J
01327                   H( JR, J ) = -H( JR, J )
01328                   T( JR, J ) = -T( JR, J )
01329   390          CONTINUE
01330             ELSE
01331                H( J, J ) = -H( J, J )
01332                T( J, J ) = -T( J, J )
01333             END IF
01334             IF( ILZ ) THEN
01335                DO 400 JR = 1, N
01336                   Z( JR, J ) = -Z( JR, J )
01337   400          CONTINUE
01338             END IF
01339          END IF
01340          ALPHAR( J ) = H( J, J )
01341          ALPHAI( J ) = ZERO
01342          BETA( J ) = T( J, J )
01343   410 CONTINUE
01344 *
01345 *     Normal Termination
01346 *
01347       INFO = 0
01348 *
01349 *     Exit (other than argument error) -- return optimal workspace size
01350 *
01351   420 CONTINUE
01352       WORK( 1 ) = DBLE( N )
01353       RETURN
01354 *
01355 *     End of DHGEQZ
01356 *
01357       END
 All Files Functions