 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, 00002 $ 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 COMPLEX H( LDH, * ), W( * ), Z( LDZ, * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * CLAHQR is an auxiliary routine called by CHSEQR to update the 00020 * eigenvalues and Schur decomposition already computed by CHSEQR, 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 triangular in rows and 00041 * columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). 00042 * CLAHQR works primarily with the Hessenberg submatrix in rows 00043 * and columns ILO to IHI, but applies transformations to all of 00044 * H if WANTT is .TRUE.. 00045 * 1 <= ILO <= max(1,IHI); IHI <= N. 00046 * 00047 * H (input/output) COMPLEX array, dimension (LDH,N) 00048 * On entry, the upper Hessenberg matrix H. 00049 * On exit, if INFO is zero and if WANTT is .TRUE., then H 00050 * is upper triangular in rows and columns ILO:IHI. If INFO 00051 * is zero and if WANTT is .FALSE., then the contents of H 00052 * are unspecified on exit. The output state of H in case 00053 * INF is positive is below under the description of INFO. 00054 * 00055 * LDH (input) INTEGER 00056 * The leading dimension of the array H. LDH >= max(1,N). 00057 * 00058 * W (output) COMPLEX array, dimension (N) 00059 * The computed eigenvalues ILO to IHI are stored in the 00060 * corresponding elements of W. If WANTT is .TRUE., the 00061 * eigenvalues are stored in the same order as on the diagonal 00062 * of the Schur form returned in H, with W(i) = H(i,i). 00063 * 00064 * ILOZ (input) INTEGER 00065 * IHIZ (input) INTEGER 00066 * Specify the rows of Z to which transformations must be 00067 * applied if WANTZ is .TRUE.. 00068 * 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. 00069 * 00070 * Z (input/output) COMPLEX array, dimension (LDZ,N) 00071 * If WANTZ is .TRUE., on entry Z must contain the current 00072 * matrix Z of transformations accumulated by CHSEQR, and on 00073 * exit Z has been updated; transformations are applied only to 00074 * the submatrix Z(ILOZ:IHIZ,ILO:IHI). 00075 * If WANTZ is .FALSE., Z is not referenced. 00076 * 00077 * LDZ (input) INTEGER 00078 * The leading dimension of the array Z. LDZ >= max(1,N). 00079 * 00080 * INFO (output) INTEGER 00081 * = 0: successful exit 00082 * .GT. 0: if INFO = i, CLAHQR failed to compute all the 00083 * eigenvalues ILO to IHI in a total of 30 iterations 00084 * per eigenvalue; elements i+1:ihi of W contain 00085 * those eigenvalues which have been successfully 00086 * computed. 00087 * 00088 * If INFO .GT. 0 and WANTT is .FALSE., then on exit, 00089 * the remaining unconverged eigenvalues are the 00090 * eigenvalues of the upper Hessenberg matrix 00091 * rows and columns ILO thorugh INFO of the final, 00092 * output value of H. 00093 * 00094 * If INFO .GT. 0 and WANTT is .TRUE., then on exit 00095 * (*) (initial value of H)*U = U*(final value of H) 00096 * where U is an orthognal matrix. The final 00097 * value of H is upper Hessenberg and triangular in 00098 * rows and columns INFO+1 through IHI. 00099 * 00100 * If INFO .GT. 0 and WANTZ is .TRUE., then on exit 00101 * (final value of Z) = (initial value of Z)*U 00102 * where U is the orthogonal matrix in (*) 00103 * (regardless of the value of WANTT.) 00104 * 00105 * Further Details 00106 * =============== 00107 * 00108 * 02-96 Based on modifications by 00109 * David Day, Sandia National Laboratory, USA 00110 * 00111 * 12-04 Further modifications by 00112 * Ralph Byers, University of Kansas, USA 00113 * This is a modified version of CLAHQR from LAPACK version 3.0. 00114 * It is (1) more robust against overflow and underflow and 00115 * (2) adopts the more conservative Ahues & Tisseur stopping 00116 * criterion (LAWN 122, 1997). 00117 * 00118 * ========================================================= 00119 * 00120 * .. Parameters .. 00121 INTEGER ITMAX 00122 PARAMETER ( ITMAX = 30 ) 00123 COMPLEX ZERO, ONE 00124 PARAMETER ( ZERO = ( 0.0e0, 0.0e0 ), 00125 $ ONE = ( 1.0e0, 0.0e0 ) ) 00126 REAL RZERO, RONE, HALF 00127 PARAMETER ( RZERO = 0.0e0, RONE = 1.0e0, HALF = 0.5e0 ) 00128 REAL DAT1 00129 PARAMETER ( DAT1 = 3.0e0 / 4.0e0 ) 00130 * .. 00131 * .. Local Scalars .. 00132 COMPLEX CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U, 00133 $ V2, X, Y 00134 REAL AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX, 00135 $ SAFMIN, SMLNUM, SX, T2, TST, ULP 00136 INTEGER I, I1, I2, ITS, J, JHI, JLO, K, L, M, NH, NZ 00137 * .. 00138 * .. Local Arrays .. 00139 COMPLEX V( 2 ) 00140 * .. 00141 * .. External Functions .. 00142 COMPLEX CLADIV 00143 REAL SLAMCH 00144 EXTERNAL CLADIV, SLAMCH 00145 * .. 00146 * .. External Subroutines .. 00147 EXTERNAL CCOPY, CLARFG, CSCAL, SLABAD 00148 * .. 00149 * .. Statement Functions .. 00150 REAL CABS1 00151 * .. 00152 * .. Intrinsic Functions .. 00153 INTRINSIC ABS, AIMAG, CONJG, MAX, MIN, REAL, SQRT 00154 * .. 00155 * .. Statement Function definitions .. 00156 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) 00157 * .. 00158 * .. Executable Statements .. 00159 * 00160 INFO = 0 00161 * 00162 * Quick return if possible 00163 * 00164 IF( N.EQ.0 ) 00165 $ RETURN 00166 IF( ILO.EQ.IHI ) THEN 00167 W( ILO ) = H( ILO, ILO ) 00168 RETURN 00169 END IF 00170 * 00171 * ==== clear out the trash ==== 00172 DO 10 J = ILO, IHI - 3 00173 H( J+2, J ) = ZERO 00174 H( J+3, J ) = ZERO 00175 10 CONTINUE 00176 IF( ILO.LE.IHI-2 ) 00177 $ H( IHI, IHI-2 ) = ZERO 00178 * ==== ensure that subdiagonal entries are real ==== 00179 IF( WANTT ) THEN 00180 JLO = 1 00181 JHI = N 00182 ELSE 00183 JLO = ILO 00184 JHI = IHI 00185 END IF 00186 DO 20 I = ILO + 1, IHI 00187 IF( AIMAG( H( I, I-1 ) ).NE.RZERO ) THEN 00188 * ==== The following redundant normalization 00189 * . avoids problems with both gradual and 00190 * . sudden underflow in ABS(H(I,I-1)) ==== 00191 SC = H( I, I-1 ) / CABS1( H( I, I-1 ) ) 00192 SC = CONJG( SC ) / ABS( SC ) 00193 H( I, I-1 ) = ABS( H( I, I-1 ) ) 00194 CALL CSCAL( JHI-I+1, SC, H( I, I ), LDH ) 00195 CALL CSCAL( MIN( JHI, I+1 )-JLO+1, CONJG( SC ), H( JLO, I ), 00196 $ 1 ) 00197 IF( WANTZ ) 00198 $ CALL CSCAL( IHIZ-ILOZ+1, CONJG( SC ), Z( ILOZ, I ), 1 ) 00199 END IF 00200 20 CONTINUE 00201 * 00202 NH = IHI - ILO + 1 00203 NZ = IHIZ - ILOZ + 1 00204 * 00205 * Set machine-dependent constants for the stopping criterion. 00206 * 00207 SAFMIN = SLAMCH( 'SAFE MINIMUM' ) 00208 SAFMAX = RONE / SAFMIN 00209 CALL SLABAD( SAFMIN, SAFMAX ) 00210 ULP = SLAMCH( 'PRECISION' ) 00211 SMLNUM = SAFMIN*( REAL( NH ) / ULP ) 00212 * 00213 * I1 and I2 are the indices of the first row and last column of H 00214 * to which transformations must be applied. If eigenvalues only are 00215 * being computed, I1 and I2 are set inside the main loop. 00216 * 00217 IF( WANTT ) THEN 00218 I1 = 1 00219 I2 = N 00220 END IF 00221 * 00222 * The main loop begins here. I is the loop index and decreases from 00223 * IHI to ILO in steps of 1. Each iteration of the loop works 00224 * with the active submatrix in rows and columns L to I. 00225 * Eigenvalues I+1 to IHI have already converged. Either L = ILO, or 00226 * H(L,L-1) is negligible so that the matrix splits. 00227 * 00228 I = IHI 00229 30 CONTINUE 00230 IF( I.LT.ILO ) 00231 $ GO TO 150 00232 * 00233 * Perform QR iterations on rows and columns ILO to I until a 00234 * submatrix of order 1 splits off at the bottom because a 00235 * subdiagonal element has become negligible. 00236 * 00237 L = ILO 00238 DO 130 ITS = 0, ITMAX 00239 * 00240 * Look for a single small subdiagonal element. 00241 * 00242 DO 40 K = I, L + 1, -1 00243 IF( CABS1( H( K, K-1 ) ).LE.SMLNUM ) 00244 $ GO TO 50 00245 TST = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) ) 00246 IF( TST.EQ.ZERO ) THEN 00247 IF( K-2.GE.ILO ) 00248 $ TST = TST + ABS( REAL( H( K-1, K-2 ) ) ) 00249 IF( K+1.LE.IHI ) 00250 $ TST = TST + ABS( REAL( H( K+1, K ) ) ) 00251 END IF 00252 * ==== The following is a conservative small subdiagonal 00253 * . deflation criterion due to Ahues & Tisseur (LAWN 122, 00254 * . 1997). It has better mathematical foundation and 00255 * . improves accuracy in some examples. ==== 00256 IF( ABS( REAL( H( K, K-1 ) ) ).LE.ULP*TST ) THEN 00257 AB = MAX( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) ) 00258 BA = MIN( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) ) 00259 AA = MAX( CABS1( H( K, K ) ), 00260 $ CABS1( H( K-1, K-1 )-H( K, K ) ) ) 00261 BB = MIN( CABS1( H( K, K ) ), 00262 $ CABS1( H( K-1, K-1 )-H( K, K ) ) ) 00263 S = AA + AB 00264 IF( BA*( AB / S ).LE.MAX( SMLNUM, 00265 $ ULP*( BB*( AA / S ) ) ) )GO TO 50 00266 END IF 00267 40 CONTINUE 00268 50 CONTINUE 00269 L = K 00270 IF( L.GT.ILO ) THEN 00271 * 00272 * H(L,L-1) is negligible 00273 * 00274 H( L, L-1 ) = ZERO 00275 END IF 00276 * 00277 * Exit from loop if a submatrix of order 1 has split off. 00278 * 00279 IF( L.GE.I ) 00280 $ GO TO 140 00281 * 00282 * Now the active submatrix is in rows and columns L to I. If 00283 * eigenvalues only are being computed, only the active submatrix 00284 * need be transformed. 00285 * 00286 IF( .NOT.WANTT ) THEN 00287 I1 = L 00288 I2 = I 00289 END IF 00290 * 00291 IF( ITS.EQ.10 ) THEN 00292 * 00293 * Exceptional shift. 00294 * 00295 S = DAT1*ABS( REAL( H( L+1, L ) ) ) 00296 T = S + H( L, L ) 00297 ELSE IF( ITS.EQ.20 ) THEN 00298 * 00299 * Exceptional shift. 00300 * 00301 S = DAT1*ABS( REAL( H( I, I-1 ) ) ) 00302 T = S + H( I, I ) 00303 ELSE 00304 * 00305 * Wilkinson's shift. 00306 * 00307 T = H( I, I ) 00308 U = SQRT( H( I-1, I ) )*SQRT( H( I, I-1 ) ) 00309 S = CABS1( U ) 00310 IF( S.NE.RZERO ) THEN 00311 X = HALF*( H( I-1, I-1 )-T ) 00312 SX = CABS1( X ) 00313 S = MAX( S, CABS1( X ) ) 00314 Y = S*SQRT( ( X / S )**2+( U / S )**2 ) 00315 IF( SX.GT.RZERO ) THEN 00316 IF( REAL( X / SX )*REAL( Y )+AIMAG( X / SX )* 00317 $ AIMAG( Y ).LT.RZERO )Y = -Y 00318 END IF 00319 T = T - U*CLADIV( U, ( X+Y ) ) 00320 END IF 00321 END IF 00322 * 00323 * Look for two consecutive small subdiagonal elements. 00324 * 00325 DO 60 M = I - 1, L + 1, -1 00326 * 00327 * Determine the effect of starting the single-shift QR 00328 * iteration at row M, and see if this would make H(M,M-1) 00329 * negligible. 00330 * 00331 H11 = H( M, M ) 00332 H22 = H( M+1, M+1 ) 00333 H11S = H11 - T 00334 H21 = REAL( H( M+1, M ) ) 00335 S = CABS1( H11S ) + ABS( H21 ) 00336 H11S = H11S / S 00337 H21 = H21 / S 00338 V( 1 ) = H11S 00339 V( 2 ) = H21 00340 H10 = REAL( H( M, M-1 ) ) 00341 IF( ABS( H10 )*ABS( H21 ).LE.ULP* 00342 $ ( CABS1( H11S )*( CABS1( H11 )+CABS1( H22 ) ) ) ) 00343 $ GO TO 70 00344 60 CONTINUE 00345 H11 = H( L, L ) 00346 H22 = H( L+1, L+1 ) 00347 H11S = H11 - T 00348 H21 = REAL( H( L+1, L ) ) 00349 S = CABS1( H11S ) + ABS( H21 ) 00350 H11S = H11S / S 00351 H21 = H21 / S 00352 V( 1 ) = H11S 00353 V( 2 ) = H21 00354 70 CONTINUE 00355 * 00356 * Single-shift QR step 00357 * 00358 DO 120 K = M, I - 1 00359 * 00360 * The first iteration of this loop determines a reflection G 00361 * from the vector V and applies it from left and right to H, 00362 * thus creating a nonzero bulge below the subdiagonal. 00363 * 00364 * Each subsequent iteration determines a reflection G to 00365 * restore the Hessenberg form in the (K-1)th column, and thus 00366 * chases the bulge one step toward the bottom of the active 00367 * submatrix. 00368 * 00369 * V(2) is always real before the call to CLARFG, and hence 00370 * after the call T2 ( = T1*V(2) ) is also real. 00371 * 00372 IF( K.GT.M ) 00373 $ CALL CCOPY( 2, H( K, K-1 ), 1, V, 1 ) 00374 CALL CLARFG( 2, V( 1 ), V( 2 ), 1, T1 ) 00375 IF( K.GT.M ) THEN 00376 H( K, K-1 ) = V( 1 ) 00377 H( K+1, K-1 ) = ZERO 00378 END IF 00379 V2 = V( 2 ) 00380 T2 = REAL( T1*V2 ) 00381 * 00382 * Apply G from the left to transform the rows of the matrix 00383 * in columns K to I2. 00384 * 00385 DO 80 J = K, I2 00386 SUM = CONJG( T1 )*H( K, J ) + T2*H( K+1, J ) 00387 H( K, J ) = H( K, J ) - SUM 00388 H( K+1, J ) = H( K+1, J ) - SUM*V2 00389 80 CONTINUE 00390 * 00391 * Apply G from the right to transform the columns of the 00392 * matrix in rows I1 to min(K+2,I). 00393 * 00394 DO 90 J = I1, MIN( K+2, I ) 00395 SUM = T1*H( J, K ) + T2*H( J, K+1 ) 00396 H( J, K ) = H( J, K ) - SUM 00397 H( J, K+1 ) = H( J, K+1 ) - SUM*CONJG( V2 ) 00398 90 CONTINUE 00399 * 00400 IF( WANTZ ) THEN 00401 * 00402 * Accumulate transformations in the matrix Z 00403 * 00404 DO 100 J = ILOZ, IHIZ 00405 SUM = T1*Z( J, K ) + T2*Z( J, K+1 ) 00406 Z( J, K ) = Z( J, K ) - SUM 00407 Z( J, K+1 ) = Z( J, K+1 ) - SUM*CONJG( V2 ) 00408 100 CONTINUE 00409 END IF 00410 * 00411 IF( K.EQ.M .AND. M.GT.L ) THEN 00412 * 00413 * If the QR step was started at row M > L because two 00414 * consecutive small subdiagonals were found, then extra 00415 * scaling must be performed to ensure that H(M,M-1) remains 00416 * real. 00417 * 00418 TEMP = ONE - T1 00419 TEMP = TEMP / ABS( TEMP ) 00420 H( M+1, M ) = H( M+1, M )*CONJG( TEMP ) 00421 IF( M+2.LE.I ) 00422 $ H( M+2, M+1 ) = H( M+2, M+1 )*TEMP 00423 DO 110 J = M, I 00424 IF( J.NE.M+1 ) THEN 00425 IF( I2.GT.J ) 00426 $ CALL CSCAL( I2-J, TEMP, H( J, J+1 ), LDH ) 00427 CALL CSCAL( J-I1, CONJG( TEMP ), H( I1, J ), 1 ) 00428 IF( WANTZ ) THEN 00429 CALL CSCAL( NZ, CONJG( TEMP ), Z( ILOZ, J ), 1 ) 00430 END IF 00431 END IF 00432 110 CONTINUE 00433 END IF 00434 120 CONTINUE 00435 * 00436 * Ensure that H(I,I-1) is real. 00437 * 00438 TEMP = H( I, I-1 ) 00439 IF( AIMAG( TEMP ).NE.RZERO ) THEN 00440 RTEMP = ABS( TEMP ) 00441 H( I, I-1 ) = RTEMP 00442 TEMP = TEMP / RTEMP 00443 IF( I2.GT.I ) 00444 $ CALL CSCAL( I2-I, CONJG( TEMP ), H( I, I+1 ), LDH ) 00445 CALL CSCAL( I-I1, TEMP, H( I1, I ), 1 ) 00446 IF( WANTZ ) THEN 00447 CALL CSCAL( NZ, TEMP, Z( ILOZ, I ), 1 ) 00448 END IF 00449 END IF 00450 * 00451 130 CONTINUE 00452 * 00453 * Failure to converge in remaining number of iterations 00454 * 00455 INFO = I 00456 RETURN 00457 * 00458 140 CONTINUE 00459 * 00460 * H(I,I-1) is negligible: one eigenvalue has converged. 00461 * 00462 W( I ) = H( I, I ) 00463 * 00464 * return to start of the main loop with new value of I. 00465 * 00466 I = L - 1 00467 GO TO 30 00468 * 00469 150 CONTINUE 00470 RETURN 00471 * 00472 * End of CLAHQR 00473 * 00474 END
1.7.3 
