1 SUBROUTINE CLAHQR2( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, 2 $ IHIZ, Z, LDZ, INFO ) 3* 4* -- ScaLAPACK routine (version 1.7) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 22, 2000 8* 9* .. Scalar Arguments .. 10 LOGICAL WANTT, WANTZ 11 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N 12* .. 13* .. Array Arguments .. 14 COMPLEX H( LDH, * ), W( * ), Z( LDZ, * ) 15* .. 16* 17* Purpose 18* ======= 19* 20* CLAHQR2 is an auxiliary routine called by CHSEQR to update the 21* eigenvalues and Schur decomposition already computed by CHSEQR, by 22* dealing with the Hessenberg submatrix in rows and columns ILO to IHI. 23* This version of CLAHQR (not the standard LAPACK version) uses a 24* double-shift algorithm (like LAPACK's SLAHQR). 25* Unlike the standard LAPACK convention, this does not assume the 26* subdiagonal is real, nor does it work to preserve this quality if 27* given. 28* 29* Arguments 30* ========= 31* 32* WANTT (input) LOGICAL 33* = .TRUE. : the full Schur form T is required; 34* = .FALSE.: only eigenvalues are required. 35* 36* WANTZ (input) LOGICAL 37* = .TRUE. : the matrix of Schur vectors Z is required; 38* = .FALSE.: Schur vectors are not required. 39* 40* N (input) INTEGER 41* The order of the matrix H. N >= 0. 42* 43* ILO (input) INTEGER 44* IHI (input) INTEGER 45* It is assumed that H is already upper triangular in rows and 46* columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). 47* CLAHQR works primarily with the Hessenberg submatrix in rows 48* and columns ILO to IHI, but applies transformations to all of 49* H if WANTT is .TRUE.. 50* 1 <= ILO <= max(1,IHI); IHI <= N. 51* 52* H (input/output) COMPLEX array, dimension (LDH,N) 53* On entry, the upper Hessenberg matrix H. 54* On exit, if WANTT is .TRUE., H is upper triangular in rows 55* and columns ILO:IHI. If WANTT is .FALSE., the contents of H 56* are unspecified on exit. 57* 58* LDH (input) INTEGER 59* The leading dimension of the array H. LDH >= max(1,N). 60* 61* W (output) COMPLEX array, dimension (N) 62* The computed eigenvalues ILO to IHI are stored in the 63* corresponding elements of W. If WANTT is .TRUE., the 64* eigenvalues are stored in the same order as on the diagonal 65* of the Schur form returned in H, with W(i) = H(i,i). 66* 67* ILOZ (input) INTEGER 68* IHIZ (input) INTEGER 69* Specify the rows of Z to which transformations must be 70* applied if WANTZ is .TRUE.. 71* 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. 72* 73* Z (input/output) COMPLEX array, dimension (LDZ,N) 74* If WANTZ is .TRUE., on entry Z must contain the current 75* matrix Z of transformations, and on exit Z has been updated; 76* transformations are applied only to the submatrix 77* Z(ILOZ:IHIZ,ILO:IHI). If WANTZ is .FALSE., Z is not 78* referenced. 79* 80* LDZ (input) INTEGER 81* The leading dimension of the array Z. LDZ >= max(1,N). 82* 83* INFO (output) INTEGER 84* = 0: successful exit 85* > 0: if INFO = i, CLAHQR failed to compute all the 86* eigenvalues ILO to IHI in a total of 30*(IHI-ILO+1) 87* iterations; elements i+1:ihi of W contain those 88* eigenvalues which have been successfully computed. 89* 90* Further Details 91* =============== 92* 93* Modified by Mark R. Fahey, June, 2000 94* 95* ===================================================================== 96* 97* .. Parameters .. 98 COMPLEX ZERO 99 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) ) 100 REAL RZERO, RONE 101 PARAMETER ( RZERO = 0.0E+0, RONE = 1.0E+0 ) 102 REAL DAT1, DAT2 103 PARAMETER ( DAT1 = 0.75E+0, DAT2 = -0.4375E+0 ) 104* .. 105* .. Local Scalars .. 106 INTEGER I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ 107 REAL CS, OVFL, S, SMLNUM, TST1, ULP, UNFL 108 COMPLEX CDUM, H00, H10, H11, H12, H21, H22, H33, H33S, 109 $ H43H34, H44, H44S, SN, SUM, T1, T2, T3, V1, V2, 110 $ V3 111* .. 112* .. Local Arrays .. 113 REAL RWORK( 1 ) 114 COMPLEX V( 3 ) 115* .. 116* .. External Functions .. 117 REAL SLAMCH, CLANHS 118 EXTERNAL SLAMCH, CLANHS 119* .. 120* .. External Subroutines .. 121 EXTERNAL SLABAD, CCOPY, CLANV2, CLARFG, CROT 122* .. 123* .. Intrinsic Functions .. 124 INTRINSIC ABS, REAL, CONJG, AIMAG, MAX, MIN 125* .. 126* .. Statement Functions .. 127 REAL CABS1 128* .. 129* .. Statement Function definitions .. 130 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) 131* .. 132* .. Executable Statements .. 133* 134 INFO = 0 135* 136* Quick return if possible 137* 138 IF( N.EQ.0 ) 139 $ RETURN 140 IF( ILO.EQ.IHI ) THEN 141 W( ILO ) = H( ILO, ILO ) 142 RETURN 143 END IF 144* 145 NH = IHI - ILO + 1 146 NZ = IHIZ - ILOZ + 1 147* 148* Set machine-dependent constants for the stopping criterion. 149* If norm(H) <= sqrt(OVFL), overflow should not occur. 150* 151 UNFL = SLAMCH( 'Safe minimum' ) 152 OVFL = RONE / UNFL 153 CALL SLABAD( UNFL, OVFL ) 154 ULP = SLAMCH( 'Precision' ) 155 SMLNUM = UNFL*( NH / ULP ) 156* 157* I1 and I2 are the indices of the first row and last column of H 158* to which transformations must be applied. If eigenvalues only are 159* being computed, I1 and I2 are set inside the main loop. 160* 161 IF( WANTT ) THEN 162 I1 = 1 163 I2 = N 164 END IF 165* 166* ITN is the total number of QR iterations allowed. 167* 168 ITN = 30*NH 169* 170* The main loop begins here. I is the loop index and decreases from 171* IHI to ILO in steps of 1 or 2. Each iteration of the loop works 172* with the active submatrix in rows and columns L to I. 173* Eigenvalues I+1 to IHI have already converged. Either L = ILO, or 174* H(L,L-1) is negligible so that the matrix splits. 175* 176 I = IHI 177 10 CONTINUE 178 L = ILO 179 IF( I.LT.ILO ) 180 $ GO TO 150 181* 182* Perform QR iterations on rows and columns ILO to I until a 183* submatrix of order 1 or 2 splits off at the bottom because a 184* subdiagonal element has become negligible. 185* 186 DO 130 ITS = 0, ITN 187* 188* Look for a single small subdiagonal element. 189* 190 DO 20 K = I, L + 1, -1 191 TST1 = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) ) 192 IF( TST1.EQ.RZERO ) 193 $ TST1 = CLANHS( '1', I-L+1, H( L, L ), LDH, RWORK ) 194 IF( CABS1( H( K, K-1 ) ).LE.MAX( ULP*TST1, SMLNUM ) ) 195 $ GO TO 30 196 20 CONTINUE 197 30 CONTINUE 198 L = K 199 IF( L.GT.ILO ) THEN 200* 201* H(L,L-1) is negligible 202* 203 H( L, L-1 ) = ZERO 204 END IF 205* 206* Exit from loop if a submatrix of order 1 or 2 has split off. 207* 208 IF( L.GE.I-1 ) 209 $ GO TO 140 210* 211* Now the active submatrix is in rows and columns L to I. If 212* eigenvalues only are being computed, only the active submatrix 213* need be transformed. 214* 215 IF( .NOT.WANTT ) THEN 216 I1 = L 217 I2 = I 218 END IF 219* 220 IF( ITS.EQ.10 .OR. ITS.EQ.20 ) THEN 221* 222* Exceptional shift. 223* 224* S = ABS( REAL( H( I,I-1 ) ) ) + ABS( REAL( H( I-1,I-2 ) ) ) 225 S = CABS1( H( I, I-1 ) ) + CABS1( H( I-1, I-2 ) ) 226 H44 = DAT1*S 227 H33 = H44 228 H43H34 = DAT2*S*S 229 ELSE 230* 231* Prepare to use Wilkinson's shift. 232* 233 H44 = H( I, I ) 234 H33 = H( I-1, I-1 ) 235 H43H34 = H( I, I-1 )*H( I-1, I ) 236 END IF 237* 238* Look for two consecutive small subdiagonal elements. 239* 240 DO 40 M = I - 2, L, -1 241* 242* Determine the effect of starting the double-shift QR 243* iteration at row M, and see if this would make H(M,M-1) 244* negligible. 245* 246 H11 = H( M, M ) 247 H22 = H( M+1, M+1 ) 248 H21 = H( M+1, M ) 249 H12 = H( M, M+1 ) 250 H44S = H44 - H11 251 H33S = H33 - H11 252 V1 = ( H33S*H44S-H43H34 ) / H21 + H12 253 V2 = H22 - H11 - H33S - H44S 254 V3 = H( M+2, M+1 ) 255 S = CABS1( V1 ) + CABS1( V2 ) + ABS( V3 ) 256 V1 = V1 / S 257 V2 = V2 / S 258 V3 = V3 / S 259 V( 1 ) = V1 260 V( 2 ) = V2 261 V( 3 ) = V3 262 IF( M.EQ.L ) 263 $ GO TO 50 264 H00 = H( M-1, M-1 ) 265 H10 = H( M, M-1 ) 266 TST1 = CABS1( V1 )*( CABS1( H00 )+CABS1( H11 )+ 267 $ CABS1( H22 ) ) 268 IF( CABS1( H10 )*( CABS1( V2 )+CABS1( V3 ) ).LE.ULP*TST1 ) 269 $ GO TO 50 270 40 CONTINUE 271 50 CONTINUE 272* 273* Double-shift QR step 274* 275 DO 120 K = M, I - 1 276* 277* The first iteration of this loop determines a reflection G 278* from the vector V and applies it from left and right to H, 279* thus creating a nonzero bulge below the subdiagonal. 280* 281* Each subsequent iteration determines a reflection G to 282* restore the Hessenberg form in the (K-1)th column, and thus 283* chases the bulge one step toward the bottom of the active 284* submatrix. NR is the order of G 285* 286 NR = MIN( 3, I-K+1 ) 287 IF( K.GT.M ) 288 $ CALL CCOPY( NR, H( K, K-1 ), 1, V, 1 ) 289 CALL CLARFG( NR, V( 1 ), V( 2 ), 1, T1 ) 290 IF( K.GT.M ) THEN 291 H( K, K-1 ) = V( 1 ) 292 H( K+1, K-1 ) = ZERO 293 IF( K.LT.I-1 ) 294 $ H( K+2, K-1 ) = ZERO 295 ELSE IF( M.GT.L ) THEN 296* The real double-shift code uses H( K, K-1 ) = -H( K, K-1 ) 297* instead of the following. 298 H( K, K-1 ) = H( K, K-1 ) - CONJG( T1 )*H( K, K-1 ) 299 END IF 300 V2 = V( 2 ) 301 T2 = T1*V2 302 IF( NR.EQ.3 ) THEN 303 V3 = V( 3 ) 304 T3 = T1*V3 305* 306* Apply G from the left to transform the rows of the matrix 307* in columns K to I2. 308* 309 DO 60 J = K, I2 310 SUM = CONJG( T1 )*H( K, J ) + 311 $ CONJG( T2 )*H( K+1, J ) + 312 $ CONJG( T3 )*H( K+2, J ) 313 H( K, J ) = H( K, J ) - SUM 314 H( K+1, J ) = H( K+1, J ) - SUM*V2 315 H( K+2, J ) = H( K+2, J ) - SUM*V3 316 60 CONTINUE 317* 318* Apply G from the right to transform the columns of the 319* matrix in rows I1 to min(K+3,I). 320* 321 DO 70 J = I1, MIN( K+3, I ) 322 SUM = T1*H( J, K ) + T2*H( J, K+1 ) + T3*H( J, K+2 ) 323 H( J, K ) = H( J, K ) - SUM 324 H( J, K+1 ) = H( J, K+1 ) - SUM*CONJG( V2 ) 325 H( J, K+2 ) = H( J, K+2 ) - SUM*CONJG( V3 ) 326 70 CONTINUE 327* 328 IF( WANTZ ) THEN 329* 330* Accumulate transformations in the matrix Z 331* 332 DO 80 J = ILOZ, IHIZ 333 SUM = T1*Z( J, K ) + T2*Z( J, K+1 ) + 334 $ T3*Z( J, K+2 ) 335 Z( J, K ) = Z( J, K ) - SUM 336 Z( J, K+1 ) = Z( J, K+1 ) - SUM*CONJG( V2 ) 337 Z( J, K+2 ) = Z( J, K+2 ) - SUM*CONJG( V3 ) 338 80 CONTINUE 339 END IF 340 ELSE IF( NR.EQ.2 ) THEN 341* 342* Apply G from the left to transform the rows of the matrix 343* in columns K to I2. 344* 345 DO 90 J = K, I2 346 SUM = CONJG( T1 )*H( K, J ) + 347 $ CONJG( T2 )*H( K+1, J ) 348 H( K, J ) = H( K, J ) - SUM 349 H( K+1, J ) = H( K+1, J ) - SUM*V2 350 90 CONTINUE 351* 352* Apply G from the right to transform the columns of the 353* matrix in rows I1 to min(K+2,I). 354* 355 DO 100 J = I1, MIN( K+2, I ) 356 SUM = T1*H( J, K ) + T2*H( J, K+1 ) 357 H( J, K ) = H( J, K ) - SUM 358 H( J, K+1 ) = H( J, K+1 ) - SUM*CONJG( V2 ) 359 100 CONTINUE 360* 361 IF( WANTZ ) THEN 362* 363* Accumulate transformations in the matrix Z 364* 365 DO 110 J = ILOZ, IHIZ 366 SUM = T1*Z( J, K ) + T2*Z( J, K+1 ) 367 Z( J, K ) = Z( J, K ) - SUM 368 Z( J, K+1 ) = Z( J, K+1 ) - SUM*CONJG( V2 ) 369 110 CONTINUE 370 END IF 371 END IF 372* 373* Since at the start of the QR step we have for M > L 374* H( K, K-1 ) = H( K, K-1 ) - CONJG( T1 )*H( K, K-1 ) 375* then we don't need to do the following 376* IF( K.EQ.M .AND. M.GT.L ) THEN 377* If the QR step was started at row M > L because two 378* consecutive small subdiagonals were found, then H(M,M-1) 379* must also be updated by a factor of (1-T1). 380* TEMP = ONE - T1 381* H( m, m-1 ) = H( m, m-1 )*CONJG( TEMP ) 382* END IF 383 120 CONTINUE 384* 385* Ensure that H(I,I-1) is real. 386* 387 130 CONTINUE 388* 389* Failure to converge in remaining number of iterations 390* 391 INFO = I 392 RETURN 393* 394 140 CONTINUE 395* 396 IF( L.EQ.I ) THEN 397* 398* H(I,I-1) is negligible: one eigenvalue has converged. 399* 400 W( I ) = H( I, I ) 401* 402 ELSE IF( L.EQ.I-1 ) THEN 403* 404* H(I-1,I-2) is negligible: a pair of eigenvalues have converged. 405* 406* Transform the 2-by-2 submatrix to standard Schur form, 407* and compute and store the eigenvalues. 408* 409 CALL CLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ), 410 $ H( I, I ), W( I-1 ), W( I ), CS, SN ) 411* 412 IF( WANTT ) THEN 413* 414* Apply the transformation to the rest of H. 415* 416 IF( I2.GT.I ) 417 $ CALL CROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH, 418 $ CS, SN ) 419 CALL CROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, 420 $ CONJG( SN ) ) 421 END IF 422 IF( WANTZ ) THEN 423* 424* Apply the transformation to Z. 425* 426 CALL CROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, 427 $ CONJG( SN ) ) 428 END IF 429* 430 END IF 431* 432* Decrement number of remaining iterations, and return to start of 433* the main loop with new value of I. 434* 435 ITN = ITN - ITS 436 I = L - 1 437 GO TO 10 438* 439 150 CONTINUE 440 RETURN 441* 442* End of CLAHQR2 443* 444 END 445