![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SHGEQZ 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SHGEQZ + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/shgeqz.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/shgeqz.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/shgeqz.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SHGEQZ( 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 * REAL 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 *> SHGEQZ 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 SGGHRD. 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 SGGHRD 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 realGEcomputational 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 SHGEQZ( 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 REAL 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 REAL HALF, ZERO, ONE, SAFETY 00327 PARAMETER ( HALF = 0.5E+0, ZERO = 0.0E+0, ONE = 1.0E+0, 00328 $ SAFETY = 1.0E+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 REAL 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 REAL V( 3 ) 00349 * .. 00350 * .. External Functions .. 00351 LOGICAL LSAME 00352 REAL SLAMCH, SLANHS, SLAPY2, SLAPY3 00353 EXTERNAL LSAME, SLAMCH, SLANHS, SLAPY2, SLAPY3 00354 * .. 00355 * .. External Subroutines .. 00356 EXTERNAL SLAG2, SLARFG, SLARTG, SLASET, SLASV2, SROT, 00357 $ XERBLA 00358 * .. 00359 * .. Intrinsic Functions .. 00360 INTRINSIC ABS, MAX, MIN, REAL, 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( 'SHGEQZ', -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 ) = REAL( 1 ) 00441 RETURN 00442 END IF 00443 * 00444 * Initialize Q and Z 00445 * 00446 IF( ICOMPQ.EQ.3 ) 00447 $ CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) 00448 IF( ICOMPZ.EQ.3 ) 00449 $ CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) 00450 * 00451 * Machine Constants 00452 * 00453 IN = IHI + 1 - ILO 00454 SAFMIN = SLAMCH( 'S' ) 00455 SAFMAX = ONE / SAFMIN 00456 ULP = SLAMCH( 'E' )*SLAMCH( 'B' ) 00457 ANORM = SLANHS( 'F', IN, H( ILO, ILO ), LDH, WORK ) 00458 BNORM = SLANHS( '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 SLARTG( TEMP, H( JCH+1, JCH ), C, S, 00592 $ H( JCH, JCH ) ) 00593 H( JCH+1, JCH ) = ZERO 00594 CALL SROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH, 00595 $ H( JCH+1, JCH+1 ), LDH, C, S ) 00596 CALL SROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT, 00597 $ T( JCH+1, JCH+1 ), LDT, C, S ) 00598 IF( ILQ ) 00599 $ CALL SROT( 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 SLARTG( 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 SROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT, 00627 $ T( JCH+1, JCH+2 ), LDT, C, S ) 00628 CALL SROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH, 00629 $ H( JCH+1, JCH-1 ), LDH, C, S ) 00630 IF( ILQ ) 00631 $ CALL SROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 00632 $ C, S ) 00633 TEMP = H( JCH+1, JCH ) 00634 CALL SLARTG( TEMP, H( JCH+1, JCH-1 ), C, S, 00635 $ H( JCH+1, JCH ) ) 00636 H( JCH+1, JCH-1 ) = ZERO 00637 CALL SROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1, 00638 $ H( IFRSTM, JCH-1 ), 1, C, S ) 00639 CALL SROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1, 00640 $ T( IFRSTM, JCH-1 ), 1, C, S ) 00641 IF( ILZ ) 00642 $ CALL SROT( 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 SLARTG( TEMP, H( ILAST, ILAST-1 ), C, S, 00670 $ H( ILAST, ILAST ) ) 00671 H( ILAST, ILAST-1 ) = ZERO 00672 CALL SROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1, 00673 $ H( IFRSTM, ILAST-1 ), 1, C, S ) 00674 CALL SROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1, 00675 $ T( IFRSTM, ILAST-1 ), 1, C, S ) 00676 IF( ILZ ) 00677 $ CALL SROT( 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( ( REAL( 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*REAL( 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 SLAG2 is the Wilkinson shift (AEP p.512), 00757 * 00758 CALL SLAG2( 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 SLARTG( 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 SLARTG( 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 SLARTG( 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 SLASV2( 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 SROT( ILASTM+1-IFIRST, H( ILAST-1, ILAST-1 ), LDH, 00887 $ H( ILAST, ILAST-1 ), LDH, CL, SL ) 00888 CALL SROT( ILAST+1-IFRSTM, H( IFRSTM, ILAST-1 ), 1, 00889 $ H( IFRSTM, ILAST ), 1, CR, SR ) 00890 * 00891 IF( ILAST.LT.ILASTM ) 00892 $ CALL SROT( ILASTM-ILAST, T( ILAST-1, ILAST+1 ), LDT, 00893 $ T( ILAST, ILAST+1 ), LDT, CL, SL ) 00894 IF( IFRSTM.LT.ILAST-1 ) 00895 $ CALL SROT( IFIRST-IFRSTM, T( IFRSTM, ILAST-1 ), 1, 00896 $ T( IFRSTM, ILAST ), 1, CR, SR ) 00897 * 00898 IF( ILQ ) 00899 $ CALL SROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL, 00900 $ SL ) 00901 IF( ILZ ) 00902 $ CALL SROT( 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 SLAG2( 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 = SLAPY3( C12, C11R, C11I ) 00964 CZ = C12 / T1 00965 SZR = -C11R / T1 00966 SZI = -C11I / T1 00967 ELSE 00968 CZ = SLAPY2( 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 = SLAPY2( 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 = SLAPY2( 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 = SLAPY3( 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 = SLAPY2( B1R, B1I ) 01025 B2R = CQ*CZ*B22 + TEMPR*B11 01026 B2I = -TEMPI*B11 01027 B2A = SLAPY2( 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 SLARFG( 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 SLARFG( 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 SLAGBC 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 SLARTG( 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 SLARTG( 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 ) = REAL( N ) 01353 RETURN 01354 * 01355 * End of SHGEQZ 01356 * 01357 END