|
2329
|
1 SUBROUTINE DHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, |
|
|
2 $ LDZ, WORK, LWORK, INFO ) |
|
|
3 * |
|
3333
|
4 * -- LAPACK routine (version 3.0) -- |
|
2329
|
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
6 * Courant Institute, Argonne National Lab, and Rice University |
|
3333
|
7 * June 30, 1999 |
|
2329
|
8 * |
|
|
9 * .. Scalar Arguments .. |
|
|
10 CHARACTER COMPZ, JOB |
|
|
11 INTEGER IHI, ILO, INFO, LDH, LDZ, LWORK, N |
|
|
12 * .. |
|
|
13 * .. Array Arguments .. |
|
|
14 DOUBLE PRECISION H( LDH, * ), WI( * ), WORK( * ), WR( * ), |
|
|
15 $ Z( LDZ, * ) |
|
|
16 * .. |
|
|
17 * |
|
|
18 * Purpose |
|
|
19 * ======= |
|
|
20 * |
|
|
21 * DHSEQR computes the eigenvalues of a real upper Hessenberg matrix H |
|
|
22 * and, optionally, the matrices T and Z from the Schur decomposition |
|
|
23 * H = Z T Z**T, where T is an upper quasi-triangular matrix (the Schur |
|
|
24 * form), and Z is the orthogonal matrix of Schur vectors. |
|
|
25 * |
|
|
26 * Optionally Z may be postmultiplied into an input orthogonal matrix Q, |
|
|
27 * so that this routine can give the Schur factorization of a matrix A |
|
|
28 * which has been reduced to the Hessenberg form H by the orthogonal |
|
|
29 * matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. |
|
|
30 * |
|
|
31 * Arguments |
|
|
32 * ========= |
|
|
33 * |
|
|
34 * JOB (input) CHARACTER*1 |
|
|
35 * = 'E': compute eigenvalues only; |
|
|
36 * = 'S': compute eigenvalues and the Schur form T. |
|
|
37 * |
|
|
38 * COMPZ (input) CHARACTER*1 |
|
|
39 * = 'N': no Schur vectors are computed; |
|
|
40 * = 'I': Z is initialized to the unit matrix and the matrix Z |
|
|
41 * of Schur vectors of H is returned; |
|
|
42 * = 'V': Z must contain an orthogonal matrix Q on entry, and |
|
|
43 * the product Q*Z is returned. |
|
|
44 * |
|
|
45 * N (input) INTEGER |
|
|
46 * The order of the matrix H. N >= 0. |
|
|
47 * |
|
|
48 * ILO (input) INTEGER |
|
|
49 * IHI (input) INTEGER |
|
|
50 * It is assumed that H is already upper triangular in rows |
|
|
51 * and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally |
|
|
52 * set by a previous call to DGEBAL, and then passed to SGEHRD |
|
|
53 * when the matrix output by DGEBAL is reduced to Hessenberg |
|
|
54 * form. Otherwise ILO and IHI should be set to 1 and N |
|
|
55 * respectively. |
|
|
56 * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. |
|
|
57 * |
|
|
58 * H (input/output) DOUBLE PRECISION array, dimension (LDH,N) |
|
|
59 * On entry, the upper Hessenberg matrix H. |
|
|
60 * On exit, if JOB = 'S', H contains the upper quasi-triangular |
|
|
61 * matrix T from the Schur decomposition (the Schur form); |
|
|
62 * 2-by-2 diagonal blocks (corresponding to complex conjugate |
|
|
63 * pairs of eigenvalues) are returned in standard form, with |
|
|
64 * H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0. If JOB = 'E', |
|
|
65 * the contents of H are unspecified on exit. |
|
|
66 * |
|
|
67 * LDH (input) INTEGER |
|
|
68 * The leading dimension of the array H. LDH >= max(1,N). |
|
|
69 * |
|
|
70 * WR (output) DOUBLE PRECISION array, dimension (N) |
|
|
71 * WI (output) DOUBLE PRECISION array, dimension (N) |
|
|
72 * The real and imaginary parts, respectively, of the computed |
|
|
73 * eigenvalues. If two eigenvalues are computed as a complex |
|
|
74 * conjugate pair, they are stored in consecutive elements of |
|
|
75 * WR and WI, say the i-th and (i+1)th, with WI(i) > 0 and |
|
|
76 * WI(i+1) < 0. If JOB = 'S', the eigenvalues are stored in the |
|
|
77 * same order as on the diagonal of the Schur form returned in |
|
|
78 * H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 |
|
|
79 * diagonal block, WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and |
|
|
80 * WI(i+1) = -WI(i). |
|
|
81 * |
|
|
82 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) |
|
|
83 * If COMPZ = 'N': Z is not referenced. |
|
|
84 * If COMPZ = 'I': on entry, Z need not be set, and on exit, Z |
|
|
85 * contains the orthogonal matrix Z of the Schur vectors of H. |
|
|
86 * If COMPZ = 'V': on entry Z must contain an N-by-N matrix Q, |
|
|
87 * which is assumed to be equal to the unit matrix except for |
|
|
88 * the submatrix Z(ILO:IHI,ILO:IHI); on exit Z contains Q*Z. |
|
|
89 * Normally Q is the orthogonal matrix generated by DORGHR after |
|
|
90 * the call to DGEHRD which formed the Hessenberg matrix H. |
|
|
91 * |
|
|
92 * LDZ (input) INTEGER |
|
|
93 * The leading dimension of the array Z. |
|
|
94 * LDZ >= max(1,N) if COMPZ = 'I' or 'V'; LDZ >= 1 otherwise. |
|
|
95 * |
|
3333
|
96 * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) |
|
|
97 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
|
2329
|
98 * |
|
|
99 * LWORK (input) INTEGER |
|
3333
|
100 * The dimension of the array WORK. LWORK >= max(1,N). |
|
|
101 * |
|
|
102 * If LWORK = -1, then a workspace query is assumed; the routine |
|
|
103 * only calculates the optimal size of the WORK array, returns |
|
|
104 * this value as the first entry of the WORK array, and no error |
|
|
105 * message related to LWORK is issued by XERBLA. |
|
2329
|
106 * |
|
|
107 * INFO (output) INTEGER |
|
|
108 * = 0: successful exit |
|
|
109 * < 0: if INFO = -i, the i-th argument had an illegal value |
|
|
110 * > 0: if INFO = i, DHSEQR failed to compute all of the |
|
|
111 * eigenvalues in a total of 30*(IHI-ILO+1) iterations; |
|
|
112 * elements 1:ilo-1 and i+1:n of WR and WI contain those |
|
|
113 * eigenvalues which have been successfully computed. |
|
|
114 * |
|
|
115 * ===================================================================== |
|
|
116 * |
|
|
117 * .. Parameters .. |
|
|
118 DOUBLE PRECISION ZERO, ONE, TWO |
|
|
119 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) |
|
|
120 DOUBLE PRECISION CONST |
|
|
121 PARAMETER ( CONST = 1.5D+0 ) |
|
|
122 INTEGER NSMAX, LDS |
|
|
123 PARAMETER ( NSMAX = 15, LDS = NSMAX ) |
|
|
124 * .. |
|
|
125 * .. Local Scalars .. |
|
3333
|
126 LOGICAL INITZ, LQUERY, WANTT, WANTZ |
|
2329
|
127 INTEGER I, I1, I2, IERR, II, ITEMP, ITN, ITS, J, K, L, |
|
|
128 $ MAXB, NH, NR, NS, NV |
|
|
129 DOUBLE PRECISION ABSW, OVFL, SMLNUM, TAU, TEMP, TST1, ULP, UNFL |
|
|
130 * .. |
|
|
131 * .. Local Arrays .. |
|
|
132 DOUBLE PRECISION S( LDS, NSMAX ), V( NSMAX+1 ), VV( NSMAX+1 ) |
|
|
133 * .. |
|
|
134 * .. External Functions .. |
|
|
135 LOGICAL LSAME |
|
|
136 INTEGER IDAMAX, ILAENV |
|
|
137 DOUBLE PRECISION DLAMCH, DLANHS, DLAPY2 |
|
|
138 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DLANHS, DLAPY2 |
|
|
139 * .. |
|
|
140 * .. External Subroutines .. |
|
3333
|
141 EXTERNAL DCOPY, DGEMV, DLACPY, DLAHQR, DLARFG, DLARFX, |
|
|
142 $ DLASET, DSCAL, XERBLA |
|
2329
|
143 * .. |
|
|
144 * .. Intrinsic Functions .. |
|
|
145 INTRINSIC ABS, MAX, MIN |
|
|
146 * .. |
|
|
147 * .. Executable Statements .. |
|
|
148 * |
|
|
149 * Decode and test the input parameters |
|
|
150 * |
|
|
151 WANTT = LSAME( JOB, 'S' ) |
|
|
152 INITZ = LSAME( COMPZ, 'I' ) |
|
|
153 WANTZ = INITZ .OR. LSAME( COMPZ, 'V' ) |
|
|
154 * |
|
|
155 INFO = 0 |
|
3333
|
156 WORK( 1 ) = MAX( 1, N ) |
|
|
157 LQUERY = ( LWORK.EQ.-1 ) |
|
2329
|
158 IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN |
|
|
159 INFO = -1 |
|
|
160 ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN |
|
|
161 INFO = -2 |
|
|
162 ELSE IF( N.LT.0 ) THEN |
|
|
163 INFO = -3 |
|
|
164 ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN |
|
|
165 INFO = -4 |
|
|
166 ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN |
|
|
167 INFO = -5 |
|
|
168 ELSE IF( LDH.LT.MAX( 1, N ) ) THEN |
|
|
169 INFO = -7 |
|
|
170 ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. LDZ.LT.MAX( 1, N ) ) THEN |
|
|
171 INFO = -11 |
|
3333
|
172 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN |
|
|
173 INFO = -13 |
|
2329
|
174 END IF |
|
|
175 IF( INFO.NE.0 ) THEN |
|
|
176 CALL XERBLA( 'DHSEQR', -INFO ) |
|
|
177 RETURN |
|
3333
|
178 ELSE IF( LQUERY ) THEN |
|
|
179 RETURN |
|
2329
|
180 END IF |
|
|
181 * |
|
|
182 * Initialize Z, if necessary |
|
|
183 * |
|
|
184 IF( INITZ ) |
|
|
185 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) |
|
|
186 * |
|
|
187 * Store the eigenvalues isolated by DGEBAL. |
|
|
188 * |
|
|
189 DO 10 I = 1, ILO - 1 |
|
|
190 WR( I ) = H( I, I ) |
|
|
191 WI( I ) = ZERO |
|
|
192 10 CONTINUE |
|
|
193 DO 20 I = IHI + 1, N |
|
|
194 WR( I ) = H( I, I ) |
|
|
195 WI( I ) = ZERO |
|
|
196 20 CONTINUE |
|
|
197 * |
|
|
198 * Quick return if possible. |
|
|
199 * |
|
|
200 IF( N.EQ.0 ) |
|
|
201 $ RETURN |
|
|
202 IF( ILO.EQ.IHI ) THEN |
|
|
203 WR( ILO ) = H( ILO, ILO ) |
|
|
204 WI( ILO ) = ZERO |
|
|
205 RETURN |
|
|
206 END IF |
|
|
207 * |
|
|
208 * Set rows and columns ILO to IHI to zero below the first |
|
|
209 * subdiagonal. |
|
|
210 * |
|
|
211 DO 40 J = ILO, IHI - 2 |
|
|
212 DO 30 I = J + 2, N |
|
|
213 H( I, J ) = ZERO |
|
|
214 30 CONTINUE |
|
|
215 40 CONTINUE |
|
|
216 NH = IHI - ILO + 1 |
|
|
217 * |
|
|
218 * Determine the order of the multi-shift QR algorithm to be used. |
|
|
219 * |
|
|
220 NS = ILAENV( 4, 'DHSEQR', JOB // COMPZ, N, ILO, IHI, -1 ) |
|
|
221 MAXB = ILAENV( 8, 'DHSEQR', JOB // COMPZ, N, ILO, IHI, -1 ) |
|
|
222 IF( NS.LE.2 .OR. NS.GT.NH .OR. MAXB.GE.NH ) THEN |
|
|
223 * |
|
|
224 * Use the standard double-shift algorithm |
|
|
225 * |
|
|
226 CALL DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO, |
|
|
227 $ IHI, Z, LDZ, INFO ) |
|
|
228 RETURN |
|
|
229 END IF |
|
|
230 MAXB = MAX( 3, MAXB ) |
|
|
231 NS = MIN( NS, MAXB, NSMAX ) |
|
|
232 * |
|
|
233 * Now 2 < NS <= MAXB < NH. |
|
|
234 * |
|
|
235 * Set machine-dependent constants for the stopping criterion. |
|
|
236 * If norm(H) <= sqrt(OVFL), overflow should not occur. |
|
|
237 * |
|
|
238 UNFL = DLAMCH( 'Safe minimum' ) |
|
|
239 OVFL = ONE / UNFL |
|
|
240 CALL DLABAD( UNFL, OVFL ) |
|
|
241 ULP = DLAMCH( 'Precision' ) |
|
|
242 SMLNUM = UNFL*( NH / ULP ) |
|
|
243 * |
|
|
244 * I1 and I2 are the indices of the first row and last column of H |
|
|
245 * to which transformations must be applied. If eigenvalues only are |
|
|
246 * being computed, I1 and I2 are set inside the main loop. |
|
|
247 * |
|
|
248 IF( WANTT ) THEN |
|
|
249 I1 = 1 |
|
|
250 I2 = N |
|
|
251 END IF |
|
|
252 * |
|
|
253 * ITN is the total number of multiple-shift QR iterations allowed. |
|
|
254 * |
|
|
255 ITN = 30*NH |
|
|
256 * |
|
|
257 * The main loop begins here. I is the loop index and decreases from |
|
|
258 * IHI to ILO in steps of at most MAXB. Each iteration of the loop |
|
|
259 * works with the active submatrix in rows and columns L to I. |
|
|
260 * Eigenvalues I+1 to IHI have already converged. Either L = ILO or |
|
|
261 * H(L,L-1) is negligible so that the matrix splits. |
|
|
262 * |
|
|
263 I = IHI |
|
|
264 50 CONTINUE |
|
|
265 L = ILO |
|
|
266 IF( I.LT.ILO ) |
|
|
267 $ GO TO 170 |
|
|
268 * |
|
|
269 * Perform multiple-shift QR iterations on rows and columns ILO to I |
|
|
270 * until a submatrix of order at most MAXB splits off at the bottom |
|
|
271 * because a subdiagonal element has become negligible. |
|
|
272 * |
|
|
273 DO 150 ITS = 0, ITN |
|
|
274 * |
|
|
275 * Look for a single small subdiagonal element. |
|
|
276 * |
|
|
277 DO 60 K = I, L + 1, -1 |
|
|
278 TST1 = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) ) |
|
|
279 IF( TST1.EQ.ZERO ) |
|
|
280 $ TST1 = DLANHS( '1', I-L+1, H( L, L ), LDH, WORK ) |
|
|
281 IF( ABS( H( K, K-1 ) ).LE.MAX( ULP*TST1, SMLNUM ) ) |
|
|
282 $ GO TO 70 |
|
|
283 60 CONTINUE |
|
|
284 70 CONTINUE |
|
|
285 L = K |
|
|
286 IF( L.GT.ILO ) THEN |
|
|
287 * |
|
|
288 * H(L,L-1) is negligible. |
|
|
289 * |
|
|
290 H( L, L-1 ) = ZERO |
|
|
291 END IF |
|
|
292 * |
|
|
293 * Exit from loop if a submatrix of order <= MAXB has split off. |
|
|
294 * |
|
|
295 IF( L.GE.I-MAXB+1 ) |
|
|
296 $ GO TO 160 |
|
|
297 * |
|
|
298 * Now the active submatrix is in rows and columns L to I. If |
|
|
299 * eigenvalues only are being computed, only the active submatrix |
|
|
300 * need be transformed. |
|
|
301 * |
|
|
302 IF( .NOT.WANTT ) THEN |
|
|
303 I1 = L |
|
|
304 I2 = I |
|
|
305 END IF |
|
|
306 * |
|
|
307 IF( ITS.EQ.20 .OR. ITS.EQ.30 ) THEN |
|
|
308 * |
|
|
309 * Exceptional shifts. |
|
|
310 * |
|
|
311 DO 80 II = I - NS + 1, I |
|
|
312 WR( II ) = CONST*( ABS( H( II, II-1 ) )+ |
|
|
313 $ ABS( H( II, II ) ) ) |
|
|
314 WI( II ) = ZERO |
|
|
315 80 CONTINUE |
|
|
316 ELSE |
|
|
317 * |
|
|
318 * Use eigenvalues of trailing submatrix of order NS as shifts. |
|
|
319 * |
|
|
320 CALL DLACPY( 'Full', NS, NS, H( I-NS+1, I-NS+1 ), LDH, S, |
|
|
321 $ LDS ) |
|
|
322 CALL DLAHQR( .FALSE., .FALSE., NS, 1, NS, S, LDS, |
|
|
323 $ WR( I-NS+1 ), WI( I-NS+1 ), 1, NS, Z, LDZ, |
|
|
324 $ IERR ) |
|
|
325 IF( IERR.GT.0 ) THEN |
|
|
326 * |
|
|
327 * If DLAHQR failed to compute all NS eigenvalues, use the |
|
|
328 * unconverged diagonal elements as the remaining shifts. |
|
|
329 * |
|
|
330 DO 90 II = 1, IERR |
|
|
331 WR( I-NS+II ) = S( II, II ) |
|
|
332 WI( I-NS+II ) = ZERO |
|
|
333 90 CONTINUE |
|
|
334 END IF |
|
|
335 END IF |
|
|
336 * |
|
|
337 * Form the first column of (G-w(1)) (G-w(2)) . . . (G-w(ns)) |
|
|
338 * where G is the Hessenberg submatrix H(L:I,L:I) and w is |
|
|
339 * the vector of shifts (stored in WR and WI). The result is |
|
|
340 * stored in the local array V. |
|
|
341 * |
|
|
342 V( 1 ) = ONE |
|
|
343 DO 100 II = 2, NS + 1 |
|
|
344 V( II ) = ZERO |
|
|
345 100 CONTINUE |
|
|
346 NV = 1 |
|
|
347 DO 120 J = I - NS + 1, I |
|
|
348 IF( WI( J ).GE.ZERO ) THEN |
|
|
349 IF( WI( J ).EQ.ZERO ) THEN |
|
|
350 * |
|
|
351 * real shift |
|
|
352 * |
|
|
353 CALL DCOPY( NV+1, V, 1, VV, 1 ) |
|
|
354 CALL DGEMV( 'No transpose', NV+1, NV, ONE, H( L, L ), |
|
|
355 $ LDH, VV, 1, -WR( J ), V, 1 ) |
|
|
356 NV = NV + 1 |
|
|
357 ELSE IF( WI( J ).GT.ZERO ) THEN |
|
|
358 * |
|
|
359 * complex conjugate pair of shifts |
|
|
360 * |
|
|
361 CALL DCOPY( NV+1, V, 1, VV, 1 ) |
|
|
362 CALL DGEMV( 'No transpose', NV+1, NV, ONE, H( L, L ), |
|
|
363 $ LDH, V, 1, -TWO*WR( J ), VV, 1 ) |
|
|
364 ITEMP = IDAMAX( NV+1, VV, 1 ) |
|
|
365 TEMP = ONE / MAX( ABS( VV( ITEMP ) ), SMLNUM ) |
|
|
366 CALL DSCAL( NV+1, TEMP, VV, 1 ) |
|
|
367 ABSW = DLAPY2( WR( J ), WI( J ) ) |
|
|
368 TEMP = ( TEMP*ABSW )*ABSW |
|
|
369 CALL DGEMV( 'No transpose', NV+2, NV+1, ONE, |
|
|
370 $ H( L, L ), LDH, VV, 1, TEMP, V, 1 ) |
|
|
371 NV = NV + 2 |
|
|
372 END IF |
|
|
373 * |
|
|
374 * Scale V(1:NV) so that max(abs(V(i))) = 1. If V is zero, |
|
|
375 * reset it to the unit vector. |
|
|
376 * |
|
|
377 ITEMP = IDAMAX( NV, V, 1 ) |
|
|
378 TEMP = ABS( V( ITEMP ) ) |
|
|
379 IF( TEMP.EQ.ZERO ) THEN |
|
|
380 V( 1 ) = ONE |
|
|
381 DO 110 II = 2, NV |
|
|
382 V( II ) = ZERO |
|
|
383 110 CONTINUE |
|
|
384 ELSE |
|
|
385 TEMP = MAX( TEMP, SMLNUM ) |
|
|
386 CALL DSCAL( NV, ONE / TEMP, V, 1 ) |
|
|
387 END IF |
|
|
388 END IF |
|
|
389 120 CONTINUE |
|
|
390 * |
|
|
391 * Multiple-shift QR step |
|
|
392 * |
|
|
393 DO 140 K = L, I - 1 |
|
|
394 * |
|
|
395 * The first iteration of this loop determines a reflection G |
|
|
396 * from the vector V and applies it from left and right to H, |
|
|
397 * thus creating a nonzero bulge below the subdiagonal. |
|
|
398 * |
|
|
399 * Each subsequent iteration determines a reflection G to |
|
|
400 * restore the Hessenberg form in the (K-1)th column, and thus |
|
|
401 * chases the bulge one step toward the bottom of the active |
|
|
402 * submatrix. NR is the order of G. |
|
|
403 * |
|
|
404 NR = MIN( NS+1, I-K+1 ) |
|
|
405 IF( K.GT.L ) |
|
|
406 $ CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 ) |
|
|
407 CALL DLARFG( NR, V( 1 ), V( 2 ), 1, TAU ) |
|
|
408 IF( K.GT.L ) THEN |
|
|
409 H( K, K-1 ) = V( 1 ) |
|
|
410 DO 130 II = K + 1, I |
|
|
411 H( II, K-1 ) = ZERO |
|
|
412 130 CONTINUE |
|
|
413 END IF |
|
|
414 V( 1 ) = ONE |
|
|
415 * |
|
|
416 * Apply G from the left to transform the rows of the matrix in |
|
|
417 * columns K to I2. |
|
|
418 * |
|
|
419 CALL DLARFX( 'Left', NR, I2-K+1, V, TAU, H( K, K ), LDH, |
|
|
420 $ WORK ) |
|
|
421 * |
|
|
422 * Apply G from the right to transform the columns of the |
|
|
423 * matrix in rows I1 to min(K+NR,I). |
|
|
424 * |
|
|
425 CALL DLARFX( 'Right', MIN( K+NR, I )-I1+1, NR, V, TAU, |
|
|
426 $ H( I1, K ), LDH, WORK ) |
|
|
427 * |
|
|
428 IF( WANTZ ) THEN |
|
|
429 * |
|
|
430 * Accumulate transformations in the matrix Z |
|
|
431 * |
|
|
432 CALL DLARFX( 'Right', NH, NR, V, TAU, Z( ILO, K ), LDZ, |
|
|
433 $ WORK ) |
|
|
434 END IF |
|
|
435 140 CONTINUE |
|
|
436 * |
|
|
437 150 CONTINUE |
|
|
438 * |
|
|
439 * Failure to converge in remaining number of iterations |
|
|
440 * |
|
|
441 INFO = I |
|
|
442 RETURN |
|
|
443 * |
|
|
444 160 CONTINUE |
|
|
445 * |
|
|
446 * A submatrix of order <= MAXB in rows and columns L to I has split |
|
|
447 * off. Use the double-shift QR algorithm to handle it. |
|
|
448 * |
|
|
449 CALL DLAHQR( WANTT, WANTZ, N, L, I, H, LDH, WR, WI, ILO, IHI, Z, |
|
|
450 $ LDZ, INFO ) |
|
|
451 IF( INFO.GT.0 ) |
|
|
452 $ RETURN |
|
|
453 * |
|
|
454 * Decrement number of remaining iterations, and return to start of |
|
|
455 * the main loop with a new value of I. |
|
|
456 * |
|
|
457 ITN = ITN - ITS |
|
|
458 I = L - 1 |
|
|
459 GO TO 50 |
|
|
460 * |
|
|
461 170 CONTINUE |
|
3333
|
462 WORK( 1 ) = MAX( 1, N ) |
|
2329
|
463 RETURN |
|
|
464 * |
|
|
465 * End of DHSEQR |
|
|
466 * |
|
|
467 END |
