 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, 00002 $ ILOZ, IHIZ, Z, LDZ, INFO ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.2) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N 00010 LOGICAL WANTT, WANTZ 00011 * .. 00012 * .. Array Arguments .. 00013 REAL H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * SLAHQR is an auxiliary routine called by SHSEQR to update the 00020 * eigenvalues and Schur decomposition already computed by SHSEQR, by 00021 * dealing with the Hessenberg submatrix in rows and columns ILO to 00022 * IHI. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * WANTT (input) LOGICAL 00028 * = .TRUE. : the full Schur form T is required; 00029 * = .FALSE.: only eigenvalues are required. 00030 * 00031 * WANTZ (input) LOGICAL 00032 * = .TRUE. : the matrix of Schur vectors Z is required; 00033 * = .FALSE.: Schur vectors are not required. 00034 * 00035 * N (input) INTEGER 00036 * The order of the matrix H. N >= 0. 00037 * 00038 * ILO (input) INTEGER 00039 * IHI (input) INTEGER 00040 * It is assumed that H is already upper quasi-triangular in 00041 * rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless 00042 * ILO = 1). SLAHQR works primarily with the Hessenberg 00043 * submatrix in rows and columns ILO to IHI, but applies 00044 * transformations to all of H if WANTT is .TRUE.. 00045 * 1 <= ILO <= max(1,IHI); IHI <= N. 00046 * 00047 * H (input/output) REAL array, dimension (LDH,N) 00048 * On entry, the upper Hessenberg matrix H. 00049 * On exit, if INFO is zero and if WANTT is .TRUE., H is upper 00050 * quasi-triangular in rows and columns ILO:IHI, with any 00051 * 2-by-2 diagonal blocks in standard form. If INFO is zero 00052 * and WANTT is .FALSE., the contents of H are unspecified on 00053 * exit. The output state of H if INFO is nonzero is given 00054 * below under the description of INFO. 00055 * 00056 * LDH (input) INTEGER 00057 * The leading dimension of the array H. LDH >= max(1,N). 00058 * 00059 * WR (output) REAL array, dimension (N) 00060 * WI (output) REAL array, dimension (N) 00061 * The real and imaginary parts, respectively, of the computed 00062 * eigenvalues ILO to IHI are stored in the corresponding 00063 * elements of WR and WI. If two eigenvalues are computed as a 00064 * complex conjugate pair, they are stored in consecutive 00065 * elements of WR and WI, say the i-th and (i+1)th, with 00066 * WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the 00067 * eigenvalues are stored in the same order as on the diagonal 00068 * of the Schur form returned in H, with WR(i) = H(i,i), and, if 00069 * H(i:i+1,i:i+1) is a 2-by-2 diagonal block, 00070 * WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). 00071 * 00072 * ILOZ (input) INTEGER 00073 * IHIZ (input) INTEGER 00074 * Specify the rows of Z to which transformations must be 00075 * applied if WANTZ is .TRUE.. 00076 * 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. 00077 * 00078 * Z (input/output) REAL array, dimension (LDZ,N) 00079 * If WANTZ is .TRUE., on entry Z must contain the current 00080 * matrix Z of transformations accumulated by SHSEQR, and on 00081 * exit Z has been updated; transformations are applied only to 00082 * the submatrix Z(ILOZ:IHIZ,ILO:IHI). 00083 * If WANTZ is .FALSE., Z is not referenced. 00084 * 00085 * LDZ (input) INTEGER 00086 * The leading dimension of the array Z. LDZ >= max(1,N). 00087 * 00088 * INFO (output) INTEGER 00089 * = 0: successful exit 00090 * .GT. 0: If INFO = i, SLAHQR failed to compute all the 00091 * eigenvalues ILO to IHI in a total of 30 iterations 00092 * per eigenvalue; elements i+1:ihi of WR and WI 00093 * contain those eigenvalues which have been 00094 * successfully computed. 00095 * 00096 * If INFO .GT. 0 and WANTT is .FALSE., then on exit, 00097 * the remaining unconverged eigenvalues are the 00098 * eigenvalues of the upper Hessenberg matrix rows 00099 * and columns ILO thorugh INFO of the final, output 00100 * value of H. 00101 * 00102 * If INFO .GT. 0 and WANTT is .TRUE., then on exit 00103 * (*) (initial value of H)*U = U*(final value of H) 00104 * where U is an orthognal matrix. The final 00105 * value of H is upper Hessenberg and triangular in 00106 * rows and columns INFO+1 through IHI. 00107 * 00108 * If INFO .GT. 0 and WANTZ is .TRUE., then on exit 00109 * (final value of Z) = (initial value of Z)*U 00110 * where U is the orthogonal matrix in (*) 00111 * (regardless of the value of WANTT.) 00112 * 00113 * Further Details 00114 * =============== 00115 * 00116 * 02-96 Based on modifications by 00117 * David Day, Sandia National Laboratory, USA 00118 * 00119 * 12-04 Further modifications by 00120 * Ralph Byers, University of Kansas, USA 00121 * This is a modified version of SLAHQR from LAPACK version 3.0. 00122 * It is (1) more robust against overflow and underflow and 00123 * (2) adopts the more conservative Ahues & Tisseur stopping 00124 * criterion (LAWN 122, 1997). 00125 * 00126 * ========================================================= 00127 * 00128 * .. Parameters .. 00129 INTEGER ITMAX 00130 PARAMETER ( ITMAX = 30 ) 00131 REAL ZERO, ONE, TWO 00132 PARAMETER ( ZERO = 0.0e0, ONE = 1.0e0, TWO = 2.0e0 ) 00133 REAL DAT1, DAT2 00134 PARAMETER ( DAT1 = 3.0e0 / 4.0e0, DAT2 = -0.4375e0 ) 00135 * .. 00136 * .. Local Scalars .. 00137 REAL AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S, 00138 $ H22, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SAFMAX, 00139 $ SAFMIN, SMLNUM, SN, SUM, T1, T2, T3, TR, TST, 00140 $ ULP, V2, V3 00141 INTEGER I, I1, I2, ITS, J, K, L, M, NH, NR, NZ 00142 * .. 00143 * .. Local Arrays .. 00144 REAL V( 3 ) 00145 * .. 00146 * .. External Functions .. 00147 REAL SLAMCH 00148 EXTERNAL SLAMCH 00149 * .. 00150 * .. External Subroutines .. 00151 EXTERNAL SCOPY, SLABAD, SLANV2, SLARFG, SROT 00152 * .. 00153 * .. Intrinsic Functions .. 00154 INTRINSIC ABS, MAX, MIN, REAL, SQRT 00155 * .. 00156 * .. Executable Statements .. 00157 * 00158 INFO = 0 00159 * 00160 * Quick return if possible 00161 * 00162 IF( N.EQ.0 ) 00163 $ RETURN 00164 IF( ILO.EQ.IHI ) THEN 00165 WR( ILO ) = H( ILO, ILO ) 00166 WI( ILO ) = ZERO 00167 RETURN 00168 END IF 00169 * 00170 * ==== clear out the trash ==== 00171 DO 10 J = ILO, IHI - 3 00172 H( J+2, J ) = ZERO 00173 H( J+3, J ) = ZERO 00174 10 CONTINUE 00175 IF( ILO.LE.IHI-2 ) 00176 $ H( IHI, IHI-2 ) = ZERO 00177 * 00178 NH = IHI - ILO + 1 00179 NZ = IHIZ - ILOZ + 1 00180 * 00181 * Set machine-dependent constants for the stopping criterion. 00182 * 00183 SAFMIN = SLAMCH( 'SAFE MINIMUM' ) 00184 SAFMAX = ONE / SAFMIN 00185 CALL SLABAD( SAFMIN, SAFMAX ) 00186 ULP = SLAMCH( 'PRECISION' ) 00187 SMLNUM = SAFMIN*( REAL( NH ) / ULP ) 00188 * 00189 * I1 and I2 are the indices of the first row and last column of H 00190 * to which transformations must be applied. If eigenvalues only are 00191 * being computed, I1 and I2 are set inside the main loop. 00192 * 00193 IF( WANTT ) THEN 00194 I1 = 1 00195 I2 = N 00196 END IF 00197 * 00198 * The main loop begins here. I is the loop index and decreases from 00199 * IHI to ILO in steps of 1 or 2. Each iteration of the loop works 00200 * with the active submatrix in rows and columns L to I. 00201 * Eigenvalues I+1 to IHI have already converged. Either L = ILO or 00202 * H(L,L-1) is negligible so that the matrix splits. 00203 * 00204 I = IHI 00205 20 CONTINUE 00206 L = ILO 00207 IF( I.LT.ILO ) 00208 $ GO TO 160 00209 * 00210 * Perform QR iterations on rows and columns ILO to I until a 00211 * submatrix of order 1 or 2 splits off at the bottom because a 00212 * subdiagonal element has become negligible. 00213 * 00214 DO 140 ITS = 0, ITMAX 00215 * 00216 * Look for a single small subdiagonal element. 00217 * 00218 DO 30 K = I, L + 1, -1 00219 IF( ABS( H( K, K-1 ) ).LE.SMLNUM ) 00220 $ GO TO 40 00221 TST = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) ) 00222 IF( TST.EQ.ZERO ) THEN 00223 IF( K-2.GE.ILO ) 00224 $ TST = TST + ABS( H( K-1, K-2 ) ) 00225 IF( K+1.LE.IHI ) 00226 $ TST = TST + ABS( H( K+1, K ) ) 00227 END IF 00228 * ==== The following is a conservative small subdiagonal 00229 * . deflation criterion due to Ahues & Tisseur (LAWN 122, 00230 * . 1997). It has better mathematical foundation and 00231 * . improves accuracy in some cases. ==== 00232 IF( ABS( H( K, K-1 ) ).LE.ULP*TST ) THEN 00233 AB = MAX( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) ) 00234 BA = MIN( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) ) 00235 AA = MAX( ABS( H( K, K ) ), 00236 $ ABS( H( K-1, K-1 )-H( K, K ) ) ) 00237 BB = MIN( ABS( H( K, K ) ), 00238 $ ABS( H( K-1, K-1 )-H( K, K ) ) ) 00239 S = AA + AB 00240 IF( BA*( AB / S ).LE.MAX( SMLNUM, 00241 $ ULP*( BB*( AA / S ) ) ) )GO TO 40 00242 END IF 00243 30 CONTINUE 00244 40 CONTINUE 00245 L = K 00246 IF( L.GT.ILO ) THEN 00247 * 00248 * H(L,L-1) is negligible 00249 * 00250 H( L, L-1 ) = ZERO 00251 END IF 00252 * 00253 * Exit from loop if a submatrix of order 1 or 2 has split off. 00254 * 00255 IF( L.GE.I-1 ) 00256 $ GO TO 150 00257 * 00258 * Now the active submatrix is in rows and columns L to I. If 00259 * eigenvalues only are being computed, only the active submatrix 00260 * need be transformed. 00261 * 00262 IF( .NOT.WANTT ) THEN 00263 I1 = L 00264 I2 = I 00265 END IF 00266 * 00267 IF( ITS.EQ.10 ) THEN 00268 * 00269 * Exceptional shift. 00270 * 00271 S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) ) 00272 H11 = DAT1*S + H( L, L ) 00273 H12 = DAT2*S 00274 H21 = S 00275 H22 = H11 00276 ELSE IF( ITS.EQ.20 ) THEN 00277 * 00278 * Exceptional shift. 00279 * 00280 S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) ) 00281 H11 = DAT1*S + H( I, I ) 00282 H12 = DAT2*S 00283 H21 = S 00284 H22 = H11 00285 ELSE 00286 * 00287 * Prepare to use Francis' double shift 00288 * (i.e. 2nd degree generalized Rayleigh quotient) 00289 * 00290 H11 = H( I-1, I-1 ) 00291 H21 = H( I, I-1 ) 00292 H12 = H( I-1, I ) 00293 H22 = H( I, I ) 00294 END IF 00295 S = ABS( H11 ) + ABS( H12 ) + ABS( H21 ) + ABS( H22 ) 00296 IF( S.EQ.ZERO ) THEN 00297 RT1R = ZERO 00298 RT1I = ZERO 00299 RT2R = ZERO 00300 RT2I = ZERO 00301 ELSE 00302 H11 = H11 / S 00303 H21 = H21 / S 00304 H12 = H12 / S 00305 H22 = H22 / S 00306 TR = ( H11+H22 ) / TWO 00307 DET = ( H11-TR )*( H22-TR ) - H12*H21 00308 RTDISC = SQRT( ABS( DET ) ) 00309 IF( DET.GE.ZERO ) THEN 00310 * 00311 * ==== complex conjugate shifts ==== 00312 * 00313 RT1R = TR*S 00314 RT2R = RT1R 00315 RT1I = RTDISC*S 00316 RT2I = -RT1I 00317 ELSE 00318 * 00319 * ==== real shifts (use only one of them) ==== 00320 * 00321 RT1R = TR + RTDISC 00322 RT2R = TR - RTDISC 00323 IF( ABS( RT1R-H22 ).LE.ABS( RT2R-H22 ) ) THEN 00324 RT1R = RT1R*S 00325 RT2R = RT1R 00326 ELSE 00327 RT2R = RT2R*S 00328 RT1R = RT2R 00329 END IF 00330 RT1I = ZERO 00331 RT2I = ZERO 00332 END IF 00333 END IF 00334 * 00335 * Look for two consecutive small subdiagonal elements. 00336 * 00337 DO 50 M = I - 2, L, -1 00338 * Determine the effect of starting the double-shift QR 00339 * iteration at row M, and see if this would make H(M,M-1) 00340 * negligible. (The following uses scaling to avoid 00341 * overflows and most underflows.) 00342 * 00343 H21S = H( M+1, M ) 00344 S = ABS( H( M, M )-RT2R ) + ABS( RT2I ) + ABS( H21S ) 00345 H21S = H( M+1, M ) / S 00346 V( 1 ) = H21S*H( M, M+1 ) + ( H( M, M )-RT1R )* 00347 $ ( ( H( M, M )-RT2R ) / S ) - RT1I*( RT2I / S ) 00348 V( 2 ) = H21S*( H( M, M )+H( M+1, M+1 )-RT1R-RT2R ) 00349 V( 3 ) = H21S*H( M+2, M+1 ) 00350 S = ABS( V( 1 ) ) + ABS( V( 2 ) ) + ABS( V( 3 ) ) 00351 V( 1 ) = V( 1 ) / S 00352 V( 2 ) = V( 2 ) / S 00353 V( 3 ) = V( 3 ) / S 00354 IF( M.EQ.L ) 00355 $ GO TO 60 00356 IF( ABS( H( M, M-1 ) )*( ABS( V( 2 ) )+ABS( V( 3 ) ) ).LE. 00357 $ ULP*ABS( V( 1 ) )*( ABS( H( M-1, M-1 ) )+ABS( H( M, 00358 $ M ) )+ABS( H( M+1, M+1 ) ) ) )GO TO 60 00359 50 CONTINUE 00360 60 CONTINUE 00361 * 00362 * Double-shift QR step 00363 * 00364 DO 130 K = M, I - 1 00365 * 00366 * The first iteration of this loop determines a reflection G 00367 * from the vector V and applies it from left and right to H, 00368 * thus creating a nonzero bulge below the subdiagonal. 00369 * 00370 * Each subsequent iteration determines a reflection G to 00371 * restore the Hessenberg form in the (K-1)th column, and thus 00372 * chases the bulge one step toward the bottom of the active 00373 * submatrix. NR is the order of G. 00374 * 00375 NR = MIN( 3, I-K+1 ) 00376 IF( K.GT.M ) 00377 $ CALL SCOPY( NR, H( K, K-1 ), 1, V, 1 ) 00378 CALL SLARFG( NR, V( 1 ), V( 2 ), 1, T1 ) 00379 IF( K.GT.M ) THEN 00380 H( K, K-1 ) = V( 1 ) 00381 H( K+1, K-1 ) = ZERO 00382 IF( K.LT.I-1 ) 00383 $ H( K+2, K-1 ) = ZERO 00384 ELSE IF( M.GT.L ) THEN 00385 * ==== Use the following instead of 00386 * . H( K, K-1 ) = -H( K, K-1 ) to 00387 * . avoid a bug when v(2) and v(3) 00388 * . underflow. ==== 00389 H( K, K-1 ) = H( K, K-1 )*( ONE-T1 ) 00390 END IF 00391 V2 = V( 2 ) 00392 T2 = T1*V2 00393 IF( NR.EQ.3 ) THEN 00394 V3 = V( 3 ) 00395 T3 = T1*V3 00396 * 00397 * Apply G from the left to transform the rows of the matrix 00398 * in columns K to I2. 00399 * 00400 DO 70 J = K, I2 00401 SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J ) 00402 H( K, J ) = H( K, J ) - SUM*T1 00403 H( K+1, J ) = H( K+1, J ) - SUM*T2 00404 H( K+2, J ) = H( K+2, J ) - SUM*T3 00405 70 CONTINUE 00406 * 00407 * Apply G from the right to transform the columns of the 00408 * matrix in rows I1 to min(K+3,I). 00409 * 00410 DO 80 J = I1, MIN( K+3, I ) 00411 SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 ) 00412 H( J, K ) = H( J, K ) - SUM*T1 00413 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 00414 H( J, K+2 ) = H( J, K+2 ) - SUM*T3 00415 80 CONTINUE 00416 * 00417 IF( WANTZ ) THEN 00418 * 00419 * Accumulate transformations in the matrix Z 00420 * 00421 DO 90 J = ILOZ, IHIZ 00422 SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 ) 00423 Z( J, K ) = Z( J, K ) - SUM*T1 00424 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 00425 Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3 00426 90 CONTINUE 00427 END IF 00428 ELSE IF( NR.EQ.2 ) THEN 00429 * 00430 * Apply G from the left to transform the rows of the matrix 00431 * in columns K to I2. 00432 * 00433 DO 100 J = K, I2 00434 SUM = H( K, J ) + V2*H( K+1, J ) 00435 H( K, J ) = H( K, J ) - SUM*T1 00436 H( K+1, J ) = H( K+1, J ) - SUM*T2 00437 100 CONTINUE 00438 * 00439 * Apply G from the right to transform the columns of the 00440 * matrix in rows I1 to min(K+3,I). 00441 * 00442 DO 110 J = I1, I 00443 SUM = H( J, K ) + V2*H( J, K+1 ) 00444 H( J, K ) = H( J, K ) - SUM*T1 00445 H( J, K+1 ) = H( J, K+1 ) - SUM*T2 00446 110 CONTINUE 00447 * 00448 IF( WANTZ ) THEN 00449 * 00450 * Accumulate transformations in the matrix Z 00451 * 00452 DO 120 J = ILOZ, IHIZ 00453 SUM = Z( J, K ) + V2*Z( J, K+1 ) 00454 Z( J, K ) = Z( J, K ) - SUM*T1 00455 Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2 00456 120 CONTINUE 00457 END IF 00458 END IF 00459 130 CONTINUE 00460 * 00461 140 CONTINUE 00462 * 00463 * Failure to converge in remaining number of iterations 00464 * 00465 INFO = I 00466 RETURN 00467 * 00468 150 CONTINUE 00469 * 00470 IF( L.EQ.I ) THEN 00471 * 00472 * H(I,I-1) is negligible: one eigenvalue has converged. 00473 * 00474 WR( I ) = H( I, I ) 00475 WI( I ) = ZERO 00476 ELSE IF( L.EQ.I-1 ) THEN 00477 * 00478 * H(I-1,I-2) is negligible: a pair of eigenvalues have converged. 00479 * 00480 * Transform the 2-by-2 submatrix to standard Schur form, 00481 * and compute and store the eigenvalues. 00482 * 00483 CALL SLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ), 00484 $ H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ), 00485 $ CS, SN ) 00486 * 00487 IF( WANTT ) THEN 00488 * 00489 * Apply the transformation to the rest of H. 00490 * 00491 IF( I2.GT.I ) 00492 $ CALL SROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH, 00493 $ CS, SN ) 00494 CALL SROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN ) 00495 END IF 00496 IF( WANTZ ) THEN 00497 * 00498 * Apply the transformation to Z. 00499 * 00500 CALL SROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN ) 00501 END IF 00502 END IF 00503 * 00504 * return to start of the main loop with new value of I. 00505 * 00506 I = L - 1 00507 GO TO 20 00508 * 00509 160 CONTINUE 00510 RETURN 00511 * 00512 * End of SLAHQR 00513 * 00514 END
1.7.3 
