 |
LAPACK
3.4.0
LAPACK: Linear Algebra PACKage
|
|
LAPACK
3.4.0
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 November 2011 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.0) -- 00308 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00309 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00310 * November 2011 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-1, ILAST ) / 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 END IF 00924 * 00925 * Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.) 00926 * 00927 * Recompute shift 00928 * 00929 CALL SLAG2( H( ILAST-1, ILAST-1 ), LDH, 00930 $ T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1, 00931 $ TEMP, WR, TEMP2, WI ) 00932 * 00933 * If standardization has perturbed the shift onto real line, 00934 * do another (real single-shift) QR step. 00935 * 00936 IF( WI.EQ.ZERO ) 00937 $ GO TO 350 00938 S1INV = ONE / S1 00939 * 00940 * Do EISPACK (QZVAL) computation of alpha and beta 00941 * 00942 A11 = H( ILAST-1, ILAST-1 ) 00943 A21 = H( ILAST, ILAST-1 ) 00944 A12 = H( ILAST-1, ILAST ) 00945 A22 = H( ILAST, ILAST ) 00946 * 00947 * Compute complex Givens rotation on right 00948 * (Assume some element of C = (sA - wB) > unfl ) 00949 * __ 00950 * (sA - wB) ( CZ -SZ ) 00951 * ( SZ CZ ) 00952 * 00953 C11R = S1*A11 - WR*B11 00954 C11I = -WI*B11 00955 C12 = S1*A12 00956 C21 = S1*A21 00957 C22R = S1*A22 - WR*B22 00958 C22I = -WI*B22 00959 * 00960 IF( ABS( C11R )+ABS( C11I )+ABS( C12 ).GT.ABS( C21 )+ 00961 $ ABS( C22R )+ABS( C22I ) ) THEN 00962 T1 = SLAPY3( C12, C11R, C11I ) 00963 CZ = C12 / T1 00964 SZR = -C11R / T1 00965 SZI = -C11I / T1 00966 ELSE 00967 CZ = SLAPY2( C22R, C22I ) 00968 IF( CZ.LE.SAFMIN ) THEN 00969 CZ = ZERO 00970 SZR = ONE 00971 SZI = ZERO 00972 ELSE 00973 TEMPR = C22R / CZ 00974 TEMPI = C22I / CZ 00975 T1 = SLAPY2( CZ, C21 ) 00976 CZ = CZ / T1 00977 SZR = -C21*TEMPR / T1 00978 SZI = C21*TEMPI / T1 00979 END IF 00980 END IF 00981 * 00982 * Compute Givens rotation on left 00983 * 00984 * ( CQ SQ ) 00985 * ( __ ) A or B 00986 * ( -SQ CQ ) 00987 * 00988 AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 ) 00989 BN = ABS( B11 ) + ABS( B22 ) 00990 WABS = ABS( WR ) + ABS( WI ) 00991 IF( S1*AN.GT.WABS*BN ) THEN 00992 CQ = CZ*B11 00993 SQR = SZR*B22 00994 SQI = -SZI*B22 00995 ELSE 00996 A1R = CZ*A11 + SZR*A12 00997 A1I = SZI*A12 00998 A2R = CZ*A21 + SZR*A22 00999 A2I = SZI*A22 01000 CQ = SLAPY2( A1R, A1I ) 01001 IF( CQ.LE.SAFMIN ) THEN 01002 CQ = ZERO 01003 SQR = ONE 01004 SQI = ZERO 01005 ELSE 01006 TEMPR = A1R / CQ 01007 TEMPI = A1I / CQ 01008 SQR = TEMPR*A2R + TEMPI*A2I 01009 SQI = TEMPI*A2R - TEMPR*A2I 01010 END IF 01011 END IF 01012 T1 = SLAPY3( CQ, SQR, SQI ) 01013 CQ = CQ / T1 01014 SQR = SQR / T1 01015 SQI = SQI / T1 01016 * 01017 * Compute diagonal elements of QBZ 01018 * 01019 TEMPR = SQR*SZR - SQI*SZI 01020 TEMPI = SQR*SZI + SQI*SZR 01021 B1R = CQ*CZ*B11 + TEMPR*B22 01022 B1I = TEMPI*B22 01023 B1A = SLAPY2( B1R, B1I ) 01024 B2R = CQ*CZ*B22 + TEMPR*B11 01025 B2I = -TEMPI*B11 01026 B2A = SLAPY2( B2R, B2I ) 01027 * 01028 * Normalize so beta > 0, and Im( alpha1 ) > 0 01029 * 01030 BETA( ILAST-1 ) = B1A 01031 BETA( ILAST ) = B2A 01032 ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV 01033 ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV 01034 ALPHAR( ILAST ) = ( WR*B2A )*S1INV 01035 ALPHAI( ILAST ) = -( WI*B2A )*S1INV 01036 * 01037 * Step 3: Go to next block -- exit if finished. 01038 * 01039 ILAST = IFIRST - 1 01040 IF( ILAST.LT.ILO ) 01041 $ GO TO 380 01042 * 01043 * Reset counters 01044 * 01045 IITER = 0 01046 ESHIFT = ZERO 01047 IF( .NOT.ILSCHR ) THEN 01048 ILASTM = ILAST 01049 IF( IFRSTM.GT.ILAST ) 01050 $ IFRSTM = ILO 01051 END IF 01052 GO TO 350 01053 ELSE 01054 * 01055 * Usual case: 3x3 or larger block, using Francis implicit 01056 * double-shift 01057 * 01058 * 2 01059 * Eigenvalue equation is w - c w + d = 0, 01060 * 01061 * -1 2 -1 01062 * so compute 1st column of (A B ) - c A B + d 01063 * using the formula in QZIT (from EISPACK) 01064 * 01065 * We assume that the block is at least 3x3 01066 * 01067 AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) / 01068 $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) 01069 AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) / 01070 $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) 01071 AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) / 01072 $ ( BSCALE*T( ILAST, ILAST ) ) 01073 AD22 = ( ASCALE*H( ILAST, ILAST ) ) / 01074 $ ( BSCALE*T( ILAST, ILAST ) ) 01075 U12 = T( ILAST-1, ILAST ) / T( ILAST, ILAST ) 01076 AD11L = ( ASCALE*H( IFIRST, IFIRST ) ) / 01077 $ ( BSCALE*T( IFIRST, IFIRST ) ) 01078 AD21L = ( ASCALE*H( IFIRST+1, IFIRST ) ) / 01079 $ ( BSCALE*T( IFIRST, IFIRST ) ) 01080 AD12L = ( ASCALE*H( IFIRST, IFIRST+1 ) ) / 01081 $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) ) 01082 AD22L = ( ASCALE*H( IFIRST+1, IFIRST+1 ) ) / 01083 $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) ) 01084 AD32L = ( ASCALE*H( IFIRST+2, IFIRST+1 ) ) / 01085 $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) ) 01086 U12L = T( IFIRST, IFIRST+1 ) / T( IFIRST+1, IFIRST+1 ) 01087 * 01088 V( 1 ) = ( AD11-AD11L )*( AD22-AD11L ) - AD12*AD21 + 01089 $ AD21*U12*AD11L + ( AD12L-AD11L*U12L )*AD21L 01090 V( 2 ) = ( ( AD22L-AD11L )-AD21L*U12L-( AD11-AD11L )- 01091 $ ( AD22-AD11L )+AD21*U12 )*AD21L 01092 V( 3 ) = AD32L*AD21L 01093 * 01094 ISTART = IFIRST 01095 * 01096 CALL SLARFG( 3, V( 1 ), V( 2 ), 1, TAU ) 01097 V( 1 ) = ONE 01098 * 01099 * Sweep 01100 * 01101 DO 290 J = ISTART, ILAST - 2 01102 * 01103 * All but last elements: use 3x3 Householder transforms. 01104 * 01105 * Zero (j-1)st column of A 01106 * 01107 IF( J.GT.ISTART ) THEN 01108 V( 1 ) = H( J, J-1 ) 01109 V( 2 ) = H( J+1, J-1 ) 01110 V( 3 ) = H( J+2, J-1 ) 01111 * 01112 CALL SLARFG( 3, H( J, J-1 ), V( 2 ), 1, TAU ) 01113 V( 1 ) = ONE 01114 H( J+1, J-1 ) = ZERO 01115 H( J+2, J-1 ) = ZERO 01116 END IF 01117 * 01118 DO 230 JC = J, ILASTM 01119 TEMP = TAU*( H( J, JC )+V( 2 )*H( J+1, JC )+V( 3 )* 01120 $ H( J+2, JC ) ) 01121 H( J, JC ) = H( J, JC ) - TEMP 01122 H( J+1, JC ) = H( J+1, JC ) - TEMP*V( 2 ) 01123 H( J+2, JC ) = H( J+2, JC ) - TEMP*V( 3 ) 01124 TEMP2 = TAU*( T( J, JC )+V( 2 )*T( J+1, JC )+V( 3 )* 01125 $ T( J+2, JC ) ) 01126 T( J, JC ) = T( J, JC ) - TEMP2 01127 T( J+1, JC ) = T( J+1, JC ) - TEMP2*V( 2 ) 01128 T( J+2, JC ) = T( J+2, JC ) - TEMP2*V( 3 ) 01129 230 CONTINUE 01130 IF( ILQ ) THEN 01131 DO 240 JR = 1, N 01132 TEMP = TAU*( Q( JR, J )+V( 2 )*Q( JR, J+1 )+V( 3 )* 01133 $ Q( JR, J+2 ) ) 01134 Q( JR, J ) = Q( JR, J ) - TEMP 01135 Q( JR, J+1 ) = Q( JR, J+1 ) - TEMP*V( 2 ) 01136 Q( JR, J+2 ) = Q( JR, J+2 ) - TEMP*V( 3 ) 01137 240 CONTINUE 01138 END IF 01139 * 01140 * Zero j-th column of B (see SLAGBC for details) 01141 * 01142 * Swap rows to pivot 01143 * 01144 ILPIVT = .FALSE. 01145 TEMP = MAX( ABS( T( J+1, J+1 ) ), ABS( T( J+1, J+2 ) ) ) 01146 TEMP2 = MAX( ABS( T( J+2, J+1 ) ), ABS( T( J+2, J+2 ) ) ) 01147 IF( MAX( TEMP, TEMP2 ).LT.SAFMIN ) THEN 01148 SCALE = ZERO 01149 U1 = ONE 01150 U2 = ZERO 01151 GO TO 250 01152 ELSE IF( TEMP.GE.TEMP2 ) THEN 01153 W11 = T( J+1, J+1 ) 01154 W21 = T( J+2, J+1 ) 01155 W12 = T( J+1, J+2 ) 01156 W22 = T( J+2, J+2 ) 01157 U1 = T( J+1, J ) 01158 U2 = T( J+2, J ) 01159 ELSE 01160 W21 = T( J+1, J+1 ) 01161 W11 = T( J+2, J+1 ) 01162 W22 = T( J+1, J+2 ) 01163 W12 = T( J+2, J+2 ) 01164 U2 = T( J+1, J ) 01165 U1 = T( J+2, J ) 01166 END IF 01167 * 01168 * Swap columns if nec. 01169 * 01170 IF( ABS( W12 ).GT.ABS( W11 ) ) THEN 01171 ILPIVT = .TRUE. 01172 TEMP = W12 01173 TEMP2 = W22 01174 W12 = W11 01175 W22 = W21 01176 W11 = TEMP 01177 W21 = TEMP2 01178 END IF 01179 * 01180 * LU-factor 01181 * 01182 TEMP = W21 / W11 01183 U2 = U2 - TEMP*U1 01184 W22 = W22 - TEMP*W12 01185 W21 = ZERO 01186 * 01187 * Compute SCALE 01188 * 01189 SCALE = ONE 01190 IF( ABS( W22 ).LT.SAFMIN ) THEN 01191 SCALE = ZERO 01192 U2 = ONE 01193 U1 = -W12 / W11 01194 GO TO 250 01195 END IF 01196 IF( ABS( W22 ).LT.ABS( U2 ) ) 01197 $ SCALE = ABS( W22 / U2 ) 01198 IF( ABS( W11 ).LT.ABS( U1 ) ) 01199 $ SCALE = MIN( SCALE, ABS( W11 / U1 ) ) 01200 * 01201 * Solve 01202 * 01203 U2 = ( SCALE*U2 ) / W22 01204 U1 = ( SCALE*U1-W12*U2 ) / W11 01205 * 01206 250 CONTINUE 01207 IF( ILPIVT ) THEN 01208 TEMP = U2 01209 U2 = U1 01210 U1 = TEMP 01211 END IF 01212 * 01213 * Compute Householder Vector 01214 * 01215 T1 = SQRT( SCALE**2+U1**2+U2**2 ) 01216 TAU = ONE + SCALE / T1 01217 VS = -ONE / ( SCALE+T1 ) 01218 V( 1 ) = ONE 01219 V( 2 ) = VS*U1 01220 V( 3 ) = VS*U2 01221 * 01222 * Apply transformations from the right. 01223 * 01224 DO 260 JR = IFRSTM, MIN( J+3, ILAST ) 01225 TEMP = TAU*( H( JR, J )+V( 2 )*H( JR, J+1 )+V( 3 )* 01226 $ H( JR, J+2 ) ) 01227 H( JR, J ) = H( JR, J ) - TEMP 01228 H( JR, J+1 ) = H( JR, J+1 ) - TEMP*V( 2 ) 01229 H( JR, J+2 ) = H( JR, J+2 ) - TEMP*V( 3 ) 01230 260 CONTINUE 01231 DO 270 JR = IFRSTM, J + 2 01232 TEMP = TAU*( T( JR, J )+V( 2 )*T( JR, J+1 )+V( 3 )* 01233 $ T( JR, J+2 ) ) 01234 T( JR, J ) = T( JR, J ) - TEMP 01235 T( JR, J+1 ) = T( JR, J+1 ) - TEMP*V( 2 ) 01236 T( JR, J+2 ) = T( JR, J+2 ) - TEMP*V( 3 ) 01237 270 CONTINUE 01238 IF( ILZ ) THEN 01239 DO 280 JR = 1, N 01240 TEMP = TAU*( Z( JR, J )+V( 2 )*Z( JR, J+1 )+V( 3 )* 01241 $ Z( JR, J+2 ) ) 01242 Z( JR, J ) = Z( JR, J ) - TEMP 01243 Z( JR, J+1 ) = Z( JR, J+1 ) - TEMP*V( 2 ) 01244 Z( JR, J+2 ) = Z( JR, J+2 ) - TEMP*V( 3 ) 01245 280 CONTINUE 01246 END IF 01247 T( J+1, J ) = ZERO 01248 T( J+2, J ) = ZERO 01249 290 CONTINUE 01250 * 01251 * Last elements: Use Givens rotations 01252 * 01253 * Rotations from the left 01254 * 01255 J = ILAST - 1 01256 TEMP = H( J, J-1 ) 01257 CALL SLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) ) 01258 H( J+1, J-1 ) = ZERO 01259 * 01260 DO 300 JC = J, ILASTM 01261 TEMP = C*H( J, JC ) + S*H( J+1, JC ) 01262 H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC ) 01263 H( J, JC ) = TEMP 01264 TEMP2 = C*T( J, JC ) + S*T( J+1, JC ) 01265 T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC ) 01266 T( J, JC ) = TEMP2 01267 300 CONTINUE 01268 IF( ILQ ) THEN 01269 DO 310 JR = 1, N 01270 TEMP = C*Q( JR, J ) + S*Q( JR, J+1 ) 01271 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) 01272 Q( JR, J ) = TEMP 01273 310 CONTINUE 01274 END IF 01275 * 01276 * Rotations from the right. 01277 * 01278 TEMP = T( J+1, J+1 ) 01279 CALL SLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) ) 01280 T( J+1, J ) = ZERO 01281 * 01282 DO 320 JR = IFRSTM, ILAST 01283 TEMP = C*H( JR, J+1 ) + S*H( JR, J ) 01284 H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J ) 01285 H( JR, J+1 ) = TEMP 01286 320 CONTINUE 01287 DO 330 JR = IFRSTM, ILAST - 1 01288 TEMP = C*T( JR, J+1 ) + S*T( JR, J ) 01289 T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J ) 01290 T( JR, J+1 ) = TEMP 01291 330 CONTINUE 01292 IF( ILZ ) THEN 01293 DO 340 JR = 1, N 01294 TEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) 01295 Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J ) 01296 Z( JR, J+1 ) = TEMP 01297 340 CONTINUE 01298 END IF 01299 * 01300 * End of Double-Shift code 01301 * 01302 END IF 01303 * 01304 GO TO 350 01305 * 01306 * End of iteration loop 01307 * 01308 350 CONTINUE 01309 360 CONTINUE 01310 * 01311 * Drop-through = non-convergence 01312 * 01313 INFO = ILAST 01314 GO TO 420 01315 * 01316 * Successful completion of all QZ steps 01317 * 01318 380 CONTINUE 01319 * 01320 * Set Eigenvalues 1:ILO-1 01321 * 01322 DO 410 J = 1, ILO - 1 01323 IF( T( J, J ).LT.ZERO ) THEN 01324 IF( ILSCHR ) THEN 01325 DO 390 JR = 1, J 01326 H( JR, J ) = -H( JR, J ) 01327 T( JR, J ) = -T( JR, J ) 01328 390 CONTINUE 01329 ELSE 01330 H( J, J ) = -H( J, J ) 01331 T( J, J ) = -T( J, J ) 01332 END IF 01333 IF( ILZ ) THEN 01334 DO 400 JR = 1, N 01335 Z( JR, J ) = -Z( JR, J ) 01336 400 CONTINUE 01337 END IF 01338 END IF 01339 ALPHAR( J ) = H( J, J ) 01340 ALPHAI( J ) = ZERO 01341 BETA( J ) = T( J, J ) 01342 410 CONTINUE 01343 * 01344 * Normal Termination 01345 * 01346 INFO = 0 01347 * 01348 * Exit (other than argument error) -- return optimal workspace size 01349 * 01350 420 CONTINUE 01351 WORK( 1 ) = REAL( N ) 01352 RETURN 01353 * 01354 * End of SHGEQZ 01355 * 01356 END
1.7.6.1 
