1
SUBROUTINE SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
2
$ ILOZ, IHIZ, Z, LDZ, INFO )
4
* -- LAPACK auxiliary routine (version 2.0) --
5
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6
* Courant Institute, Argonne National Lab, and Rice University
9
* .. Scalar Arguments ..
11
INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
13
* .. Array Arguments ..
14
REAL H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
20
* SLAHQR is an auxiliary routine called by SHSEQR to update the
21
* eigenvalues and Schur decomposition already computed by SHSEQR, by
22
* dealing with the Hessenberg submatrix in rows and columns ILO to IHI.
27
* WANTT (input) LOGICAL
28
* = .TRUE. : the full Schur form T is required;
29
* = .FALSE.: only eigenvalues are required.
31
* WANTZ (input) LOGICAL
32
* = .TRUE. : the matrix of Schur vectors Z is required;
33
* = .FALSE.: Schur vectors are not required.
36
* The order of the matrix H. N >= 0.
40
* It is assumed that H is already upper quasi-triangular in
41
* rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
42
* ILO = 1). SLAHQR works primarily with the Hessenberg
43
* submatrix in rows and columns ILO to IHI, but applies
44
* transformations to all of H if WANTT is .TRUE..
45
* 1 <= ILO <= max(1,IHI); IHI <= N.
47
* H (input/output) REAL array, dimension (LDH,N)
48
* On entry, the upper Hessenberg matrix H.
49
* On exit, if WANTT is .TRUE., H is upper quasi-triangular in
50
* rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in
51
* standard form. If WANTT is .FALSE., the contents of H are
52
* unspecified on exit.
55
* The leading dimension of the array H. LDH >= max(1,N).
57
* WR (output) REAL array, dimension (N)
58
* WI (output) REAL array, dimension (N)
59
* The real and imaginary parts, respectively, of the computed
60
* eigenvalues ILO to IHI are stored in the corresponding
61
* elements of WR and WI. If two eigenvalues are computed as a
62
* complex conjugate pair, they are stored in consecutive
63
* elements of WR and WI, say the i-th and (i+1)th, with
64
* WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
65
* eigenvalues are stored in the same order as on the diagonal
66
* of the Schur form returned in H, with WR(i) = H(i,i), and, if
67
* H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
68
* WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
70
* ILOZ (input) INTEGER
71
* IHIZ (input) INTEGER
72
* Specify the rows of Z to which transformations must be
73
* applied if WANTZ is .TRUE..
74
* 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
76
* Z (input/output) REAL array, dimension (LDZ,N)
77
* If WANTZ is .TRUE., on entry Z must contain the current
78
* matrix Z of transformations accumulated by SHSEQR, and on
79
* exit Z has been updated; transformations are applied only to
80
* the submatrix Z(ILOZ:IHIZ,ILO:IHI).
81
* If WANTZ is .FALSE., Z is not referenced.
84
* The leading dimension of the array Z. LDZ >= max(1,N).
86
* INFO (output) INTEGER
87
* = 0: successful exit
88
* > 0: SLAHQR failed to compute all the eigenvalues ILO to IHI
89
* in a total of 30*(IHI-ILO+1) iterations; if INFO = i,
90
* elements i+1:ihi of WR and WI contain those eigenvalues
91
* which have been successfully computed.
93
* =====================================================================
97
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
99
PARAMETER ( DAT1 = 0.75E+0, DAT2 = -0.4375E+0 )
101
* .. Local Scalars ..
102
INTEGER I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ
103
REAL CS, H00, H10, H11, H12, H21, H22, H33, H33S,
104
$ H43H34, H44, H44S, OVFL, S, SMLNUM, SN, SUM,
105
$ T1, T2, T3, TST1, ULP, UNFL, V1, V2, V3
108
REAL V( 3 ), WORK( 1 )
110
* .. External Functions ..
112
EXTERNAL SLAMCH, SLANHS
114
* .. External Subroutines ..
115
EXTERNAL SCOPY, SLABAD, SLANV2, SLARFG, SROT
117
* .. Intrinsic Functions ..
118
INTRINSIC ABS, MAX, MIN
120
* .. Executable Statements ..
124
* Quick return if possible
128
IF( ILO.EQ.IHI ) THEN
129
WR( ILO ) = H( ILO, ILO )
137
* Set machine-dependent constants for the stopping criterion.
138
* If norm(H) <= sqrt(OVFL), overflow should not occur.
140
UNFL = SLAMCH( 'Safe minimum' )
142
CALL SLABAD( UNFL, OVFL )
143
ULP = SLAMCH( 'Precision' )
144
SMLNUM = UNFL*( NH / ULP )
146
* I1 and I2 are the indices of the first row and last column of H
147
* to which transformations must be applied. If eigenvalues only are
148
* being computed, I1 and I2 are set inside the main loop.
155
* ITN is the total number of QR iterations allowed.
159
* The main loop begins here. I is the loop index and decreases from
160
* IHI to ILO in steps of 1 or 2. Each iteration of the loop works
161
* with the active submatrix in rows and columns L to I.
162
* Eigenvalues I+1 to IHI have already converged. Either L = ILO or
163
* H(L,L-1) is negligible so that the matrix splits.
171
* Perform QR iterations on rows and columns ILO to I until a
172
* submatrix of order 1 or 2 splits off at the bottom because a
173
* subdiagonal element has become negligible.
177
* Look for a single small subdiagonal element.
179
DO 20 K = I, L + 1, -1
180
TST1 = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) )
182
$ TST1 = SLANHS( '1', I-L+1, H( L, L ), LDH, WORK )
183
IF( ABS( H( K, K-1 ) ).LE.MAX( ULP*TST1, SMLNUM ) )
190
* H(L,L-1) is negligible
195
* Exit from loop if a submatrix of order 1 or 2 has split off.
200
* Now the active submatrix is in rows and columns L to I. If
201
* eigenvalues only are being computed, only the active submatrix
202
* need be transformed.
204
IF( .NOT.WANTT ) THEN
209
IF( ITS.EQ.10 .OR. ITS.EQ.20 ) THEN
213
S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
219
* Prepare to use Wilkinson's double shift
223
H43H34 = H( I, I-1 )*H( I-1, I )
226
* Look for two consecutive small subdiagonal elements.
228
DO 40 M = I - 2, L, -1
230
* Determine the effect of starting the double-shift QR
231
* iteration at row M, and see if this would make H(M,M-1)
240
V1 = ( H33S*H44S-H43H34 ) / H21 + H12
241
V2 = H22 - H11 - H33S - H44S
243
S = ABS( V1 ) + ABS( V2 ) + ABS( V3 )
254
TST1 = ABS( V1 )*( ABS( H00 )+ABS( H11 )+ABS( H22 ) )
255
IF( ABS( H10 )*( ABS( V2 )+ABS( V3 ) ).LE.ULP*TST1 )
260
* Double-shift QR step
264
* The first iteration of this loop determines a reflection G
265
* from the vector V and applies it from left and right to H,
266
* thus creating a nonzero bulge below the subdiagonal.
268
* Each subsequent iteration determines a reflection G to
269
* restore the Hessenberg form in the (K-1)th column, and thus
270
* chases the bulge one step toward the bottom of the active
271
* submatrix. NR is the order of G.
275
$ CALL SCOPY( NR, H( K, K-1 ), 1, V, 1 )
276
CALL SLARFG( NR, V( 1 ), V( 2 ), 1, T1 )
281
$ H( K+2, K-1 ) = ZERO
282
ELSE IF( M.GT.L ) THEN
283
H( K, K-1 ) = -H( K, K-1 )
291
* Apply G from the left to transform the rows of the matrix
292
* in columns K to I2.
295
SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J )
296
H( K, J ) = H( K, J ) - SUM*T1
297
H( K+1, J ) = H( K+1, J ) - SUM*T2
298
H( K+2, J ) = H( K+2, J ) - SUM*T3
301
* Apply G from the right to transform the columns of the
302
* matrix in rows I1 to min(K+3,I).
304
DO 70 J = I1, MIN( K+3, I )
305
SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 )
306
H( J, K ) = H( J, K ) - SUM*T1
307
H( J, K+1 ) = H( J, K+1 ) - SUM*T2
308
H( J, K+2 ) = H( J, K+2 ) - SUM*T3
313
* Accumulate transformations in the matrix Z
316
SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 )
317
Z( J, K ) = Z( J, K ) - SUM*T1
318
Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
319
Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3
322
ELSE IF( NR.EQ.2 ) THEN
324
* Apply G from the left to transform the rows of the matrix
325
* in columns K to I2.
328
SUM = H( K, J ) + V2*H( K+1, J )
329
H( K, J ) = H( K, J ) - SUM*T1
330
H( K+1, J ) = H( K+1, J ) - SUM*T2
333
* Apply G from the right to transform the columns of the
334
* matrix in rows I1 to min(K+3,I).
337
SUM = H( J, K ) + V2*H( J, K+1 )
338
H( J, K ) = H( J, K ) - SUM*T1
339
H( J, K+1 ) = H( J, K+1 ) - SUM*T2
344
* Accumulate transformations in the matrix Z
346
DO 110 J = ILOZ, IHIZ
347
SUM = Z( J, K ) + V2*Z( J, K+1 )
348
Z( J, K ) = Z( J, K ) - SUM*T1
349
Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
357
* Failure to converge in remaining number of iterations
366
* H(I,I-1) is negligible: one eigenvalue has converged.
370
ELSE IF( L.EQ.I-1 ) THEN
372
* H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
374
* Transform the 2-by-2 submatrix to standard Schur form,
375
* and compute and store the eigenvalues.
377
CALL SLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ),
378
$ H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ),
383
* Apply the transformation to the rest of H.
386
$ CALL SROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH,
388
CALL SROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN )
392
* Apply the transformation to Z.
394
CALL SROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
398
* Decrement number of remaining iterations, and return to start of
399
* the main loop with new value of I.