Mercurial > octave-nkf
comparison libcruft/lapack/dlaqr4.f @ 7034:68db500cb558
[project @ 2007-10-16 18:54:19 by jwe]
author | jwe |
---|---|
date | Tue, 16 Oct 2007 18:54:23 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
7033:f0142f2afdc6 | 7034:68db500cb558 |
---|---|
1 SUBROUTINE DLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, | |
2 $ ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO ) | |
3 * | |
4 * -- LAPACK auxiliary routine (version 3.1) -- | |
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. | |
6 * November 2006 | |
7 * | |
8 * .. Scalar Arguments .. | |
9 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N | |
10 LOGICAL WANTT, WANTZ | |
11 * .. | |
12 * .. Array Arguments .. | |
13 DOUBLE PRECISION H( LDH, * ), WI( * ), WORK( * ), WR( * ), | |
14 $ Z( LDZ, * ) | |
15 * .. | |
16 * | |
17 * This subroutine implements one level of recursion for DLAQR0. | |
18 * It is a complete implementation of the small bulge multi-shift | |
19 * QR algorithm. It may be called by DLAQR0 and, for large enough | |
20 * deflation window size, it may be called by DLAQR3. This | |
21 * subroutine is identical to DLAQR0 except that it calls DLAQR2 | |
22 * instead of DLAQR3. | |
23 * | |
24 * Purpose | |
25 * ======= | |
26 * | |
27 * DLAQR4 computes the eigenvalues of a Hessenberg matrix H | |
28 * and, optionally, the matrices T and Z from the Schur decomposition | |
29 * H = Z T Z**T, where T is an upper quasi-triangular matrix (the | |
30 * Schur form), and Z is the orthogonal matrix of Schur vectors. | |
31 * | |
32 * Optionally Z may be postmultiplied into an input orthogonal | |
33 * matrix Q so that this routine can give the Schur factorization | |
34 * of a matrix A which has been reduced to the Hessenberg form H | |
35 * by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. | |
36 * | |
37 * Arguments | |
38 * ========= | |
39 * | |
40 * WANTT (input) LOGICAL | |
41 * = .TRUE. : the full Schur form T is required; | |
42 * = .FALSE.: only eigenvalues are required. | |
43 * | |
44 * WANTZ (input) LOGICAL | |
45 * = .TRUE. : the matrix of Schur vectors Z is required; | |
46 * = .FALSE.: Schur vectors are not required. | |
47 * | |
48 * N (input) INTEGER | |
49 * The order of the matrix H. N .GE. 0. | |
50 * | |
51 * ILO (input) INTEGER | |
52 * IHI (input) INTEGER | |
53 * It is assumed that H is already upper triangular in rows | |
54 * and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1, | |
55 * H(ILO,ILO-1) is zero. ILO and IHI are normally set by a | |
56 * previous call to DGEBAL, and then passed to DGEHRD when the | |
57 * matrix output by DGEBAL is reduced to Hessenberg form. | |
58 * Otherwise, ILO and IHI should be set to 1 and N, | |
59 * respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. | |
60 * If N = 0, then ILO = 1 and IHI = 0. | |
61 * | |
62 * H (input/output) DOUBLE PRECISION array, dimension (LDH,N) | |
63 * On entry, the upper Hessenberg matrix H. | |
64 * On exit, if INFO = 0 and WANTT is .TRUE., then H contains | |
65 * the upper quasi-triangular matrix T from the Schur | |
66 * decomposition (the Schur form); 2-by-2 diagonal blocks | |
67 * (corresponding to complex conjugate pairs of eigenvalues) | |
68 * are returned in standard form, with H(i,i) = H(i+1,i+1) | |
69 * and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is | |
70 * .FALSE., then the contents of H are unspecified on exit. | |
71 * (The output value of H when INFO.GT.0 is given under the | |
72 * description of INFO below.) | |
73 * | |
74 * This subroutine may explicitly set H(i,j) = 0 for i.GT.j and | |
75 * j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. | |
76 * | |
77 * LDH (input) INTEGER | |
78 * The leading dimension of the array H. LDH .GE. max(1,N). | |
79 * | |
80 * WR (output) DOUBLE PRECISION array, dimension (IHI) | |
81 * WI (output) DOUBLE PRECISION array, dimension (IHI) | |
82 * The real and imaginary parts, respectively, of the computed | |
83 * eigenvalues of H(ILO:IHI,ILO:IHI) are stored WR(ILO:IHI) | |
84 * and WI(ILO:IHI). If two eigenvalues are computed as a | |
85 * complex conjugate pair, they are stored in consecutive | |
86 * elements of WR and WI, say the i-th and (i+1)th, with | |
87 * WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then | |
88 * the eigenvalues are stored in the same order as on the | |
89 * diagonal of the Schur form returned in H, with | |
90 * WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal | |
91 * block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and | |
92 * WI(i+1) = -WI(i). | |
93 * | |
94 * ILOZ (input) INTEGER | |
95 * IHIZ (input) INTEGER | |
96 * Specify the rows of Z to which transformations must be | |
97 * applied if WANTZ is .TRUE.. | |
98 * 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N. | |
99 * | |
100 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,IHI) | |
101 * If WANTZ is .FALSE., then Z is not referenced. | |
102 * If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is | |
103 * replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the | |
104 * orthogonal Schur factor of H(ILO:IHI,ILO:IHI). | |
105 * (The output value of Z when INFO.GT.0 is given under | |
106 * the description of INFO below.) | |
107 * | |
108 * LDZ (input) INTEGER | |
109 * The leading dimension of the array Z. if WANTZ is .TRUE. | |
110 * then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1. | |
111 * | |
112 * WORK (workspace/output) DOUBLE PRECISION array, dimension LWORK | |
113 * On exit, if LWORK = -1, WORK(1) returns an estimate of | |
114 * the optimal value for LWORK. | |
115 * | |
116 * LWORK (input) INTEGER | |
117 * The dimension of the array WORK. LWORK .GE. max(1,N) | |
118 * is sufficient, but LWORK typically as large as 6*N may | |
119 * be required for optimal performance. A workspace query | |
120 * to determine the optimal workspace size is recommended. | |
121 * | |
122 * If LWORK = -1, then DLAQR4 does a workspace query. | |
123 * In this case, DLAQR4 checks the input parameters and | |
124 * estimates the optimal workspace size for the given | |
125 * values of N, ILO and IHI. The estimate is returned | |
126 * in WORK(1). No error message related to LWORK is | |
127 * issued by XERBLA. Neither H nor Z are accessed. | |
128 * | |
129 * | |
130 * INFO (output) INTEGER | |
131 * = 0: successful exit | |
132 * .GT. 0: if INFO = i, DLAQR4 failed to compute all of | |
133 * the eigenvalues. Elements 1:ilo-1 and i+1:n of WR | |
134 * and WI contain those eigenvalues which have been | |
135 * successfully computed. (Failures are rare.) | |
136 * | |
137 * If INFO .GT. 0 and WANT is .FALSE., then on exit, | |
138 * the remaining unconverged eigenvalues are the eigen- | |
139 * values of the upper Hessenberg matrix rows and | |
140 * columns ILO through INFO of the final, output | |
141 * value of H. | |
142 * | |
143 * If INFO .GT. 0 and WANTT is .TRUE., then on exit | |
144 * | |
145 * (*) (initial value of H)*U = U*(final value of H) | |
146 * | |
147 * where U is an orthogonal matrix. The final | |
148 * value of H is upper Hessenberg and quasi-triangular | |
149 * in rows and columns INFO+1 through IHI. | |
150 * | |
151 * If INFO .GT. 0 and WANTZ is .TRUE., then on exit | |
152 * | |
153 * (final value of Z(ILO:IHI,ILOZ:IHIZ) | |
154 * = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U | |
155 * | |
156 * where U is the orthogonal matrix in (*) (regard- | |
157 * less of the value of WANTT.) | |
158 * | |
159 * If INFO .GT. 0 and WANTZ is .FALSE., then Z is not | |
160 * accessed. | |
161 * | |
162 * ================================================================ | |
163 * Based on contributions by | |
164 * Karen Braman and Ralph Byers, Department of Mathematics, | |
165 * University of Kansas, USA | |
166 * | |
167 * ================================================================ | |
168 * References: | |
169 * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR | |
170 * Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 | |
171 * Performance, SIAM Journal of Matrix Analysis, volume 23, pages | |
172 * 929--947, 2002. | |
173 * | |
174 * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR | |
175 * Algorithm Part II: Aggressive Early Deflation, SIAM Journal | |
176 * of Matrix Analysis, volume 23, pages 948--973, 2002. | |
177 * | |
178 * ================================================================ | |
179 * .. Parameters .. | |
180 * | |
181 * ==== Matrices of order NTINY or smaller must be processed by | |
182 * . DLAHQR because of insufficient subdiagonal scratch space. | |
183 * . (This is a hard limit.) ==== | |
184 * | |
185 * ==== Exceptional deflation windows: try to cure rare | |
186 * . slow convergence by increasing the size of the | |
187 * . deflation window after KEXNW iterations. ===== | |
188 * | |
189 * ==== Exceptional shifts: try to cure rare slow convergence | |
190 * . with ad-hoc exceptional shifts every KEXSH iterations. | |
191 * . The constants WILK1 and WILK2 are used to form the | |
192 * . exceptional shifts. ==== | |
193 * | |
194 INTEGER NTINY | |
195 PARAMETER ( NTINY = 11 ) | |
196 INTEGER KEXNW, KEXSH | |
197 PARAMETER ( KEXNW = 5, KEXSH = 6 ) | |
198 DOUBLE PRECISION WILK1, WILK2 | |
199 PARAMETER ( WILK1 = 0.75d0, WILK2 = -0.4375d0 ) | |
200 DOUBLE PRECISION ZERO, ONE | |
201 PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 ) | |
202 * .. | |
203 * .. Local Scalars .. | |
204 DOUBLE PRECISION AA, BB, CC, CS, DD, SN, SS, SWAP | |
205 INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS, | |
206 $ KT, KTOP, KU, KV, KWH, KWTOP, KWV, LD, LS, | |
207 $ LWKOPT, NDFL, NH, NHO, NIBBLE, NMIN, NS, NSMAX, | |
208 $ NSR, NVE, NW, NWMAX, NWR | |
209 LOGICAL NWINC, SORTED | |
210 CHARACTER JBCMPZ*2 | |
211 * .. | |
212 * .. External Functions .. | |
213 INTEGER ILAENV | |
214 EXTERNAL ILAENV | |
215 * .. | |
216 * .. Local Arrays .. | |
217 DOUBLE PRECISION ZDUM( 1, 1 ) | |
218 * .. | |
219 * .. External Subroutines .. | |
220 EXTERNAL DLACPY, DLAHQR, DLANV2, DLAQR2, DLAQR5 | |
221 * .. | |
222 * .. Intrinsic Functions .. | |
223 INTRINSIC ABS, DBLE, INT, MAX, MIN, MOD | |
224 * .. | |
225 * .. Executable Statements .. | |
226 INFO = 0 | |
227 * | |
228 * ==== Quick return for N = 0: nothing to do. ==== | |
229 * | |
230 IF( N.EQ.0 ) THEN | |
231 WORK( 1 ) = ONE | |
232 RETURN | |
233 END IF | |
234 * | |
235 * ==== Set up job flags for ILAENV. ==== | |
236 * | |
237 IF( WANTT ) THEN | |
238 JBCMPZ( 1: 1 ) = 'S' | |
239 ELSE | |
240 JBCMPZ( 1: 1 ) = 'E' | |
241 END IF | |
242 IF( WANTZ ) THEN | |
243 JBCMPZ( 2: 2 ) = 'V' | |
244 ELSE | |
245 JBCMPZ( 2: 2 ) = 'N' | |
246 END IF | |
247 * | |
248 * ==== Tiny matrices must use DLAHQR. ==== | |
249 * | |
250 IF( N.LE.NTINY ) THEN | |
251 * | |
252 * ==== Estimate optimal workspace. ==== | |
253 * | |
254 LWKOPT = 1 | |
255 IF( LWORK.NE.-1 ) | |
256 $ CALL DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, | |
257 $ ILOZ, IHIZ, Z, LDZ, INFO ) | |
258 ELSE | |
259 * | |
260 * ==== Use small bulge multi-shift QR with aggressive early | |
261 * . deflation on larger-than-tiny matrices. ==== | |
262 * | |
263 * ==== Hope for the best. ==== | |
264 * | |
265 INFO = 0 | |
266 * | |
267 * ==== NWR = recommended deflation window size. At this | |
268 * . point, N .GT. NTINY = 11, so there is enough | |
269 * . subdiagonal workspace for NWR.GE.2 as required. | |
270 * . (In fact, there is enough subdiagonal space for | |
271 * . NWR.GE.3.) ==== | |
272 * | |
273 NWR = ILAENV( 13, 'DLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) | |
274 NWR = MAX( 2, NWR ) | |
275 NWR = MIN( IHI-ILO+1, ( N-1 ) / 3, NWR ) | |
276 NW = NWR | |
277 * | |
278 * ==== NSR = recommended number of simultaneous shifts. | |
279 * . At this point N .GT. NTINY = 11, so there is at | |
280 * . enough subdiagonal workspace for NSR to be even | |
281 * . and greater than or equal to two as required. ==== | |
282 * | |
283 NSR = ILAENV( 15, 'DLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) | |
284 NSR = MIN( NSR, ( N+6 ) / 9, IHI-ILO ) | |
285 NSR = MAX( 2, NSR-MOD( NSR, 2 ) ) | |
286 * | |
287 * ==== Estimate optimal workspace ==== | |
288 * | |
289 * ==== Workspace query call to DLAQR2 ==== | |
290 * | |
291 CALL DLAQR2( WANTT, WANTZ, N, ILO, IHI, NWR+1, H, LDH, ILOZ, | |
292 $ IHIZ, Z, LDZ, LS, LD, WR, WI, H, LDH, N, H, LDH, | |
293 $ N, H, LDH, WORK, -1 ) | |
294 * | |
295 * ==== Optimal workspace = MAX(DLAQR5, DLAQR2) ==== | |
296 * | |
297 LWKOPT = MAX( 3*NSR / 2, INT( WORK( 1 ) ) ) | |
298 * | |
299 * ==== Quick return in case of workspace query. ==== | |
300 * | |
301 IF( LWORK.EQ.-1 ) THEN | |
302 WORK( 1 ) = DBLE( LWKOPT ) | |
303 RETURN | |
304 END IF | |
305 * | |
306 * ==== DLAHQR/DLAQR0 crossover point ==== | |
307 * | |
308 NMIN = ILAENV( 12, 'DLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) | |
309 NMIN = MAX( NTINY, NMIN ) | |
310 * | |
311 * ==== Nibble crossover point ==== | |
312 * | |
313 NIBBLE = ILAENV( 14, 'DLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) | |
314 NIBBLE = MAX( 0, NIBBLE ) | |
315 * | |
316 * ==== Accumulate reflections during ttswp? Use block | |
317 * . 2-by-2 structure during matrix-matrix multiply? ==== | |
318 * | |
319 KACC22 = ILAENV( 16, 'DLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) | |
320 KACC22 = MAX( 0, KACC22 ) | |
321 KACC22 = MIN( 2, KACC22 ) | |
322 * | |
323 * ==== NWMAX = the largest possible deflation window for | |
324 * . which there is sufficient workspace. ==== | |
325 * | |
326 NWMAX = MIN( ( N-1 ) / 3, LWORK / 2 ) | |
327 * | |
328 * ==== NSMAX = the Largest number of simultaneous shifts | |
329 * . for which there is sufficient workspace. ==== | |
330 * | |
331 NSMAX = MIN( ( N+6 ) / 9, 2*LWORK / 3 ) | |
332 NSMAX = NSMAX - MOD( NSMAX, 2 ) | |
333 * | |
334 * ==== NDFL: an iteration count restarted at deflation. ==== | |
335 * | |
336 NDFL = 1 | |
337 * | |
338 * ==== ITMAX = iteration limit ==== | |
339 * | |
340 ITMAX = MAX( 30, 2*KEXSH )*MAX( 10, ( IHI-ILO+1 ) ) | |
341 * | |
342 * ==== Last row and column in the active block ==== | |
343 * | |
344 KBOT = IHI | |
345 * | |
346 * ==== Main Loop ==== | |
347 * | |
348 DO 80 IT = 1, ITMAX | |
349 * | |
350 * ==== Done when KBOT falls below ILO ==== | |
351 * | |
352 IF( KBOT.LT.ILO ) | |
353 $ GO TO 90 | |
354 * | |
355 * ==== Locate active block ==== | |
356 * | |
357 DO 10 K = KBOT, ILO + 1, -1 | |
358 IF( H( K, K-1 ).EQ.ZERO ) | |
359 $ GO TO 20 | |
360 10 CONTINUE | |
361 K = ILO | |
362 20 CONTINUE | |
363 KTOP = K | |
364 * | |
365 * ==== Select deflation window size ==== | |
366 * | |
367 NH = KBOT - KTOP + 1 | |
368 IF( NDFL.LT.KEXNW .OR. NH.LT.NW ) THEN | |
369 * | |
370 * ==== Typical deflation window. If possible and | |
371 * . advisable, nibble the entire active block. | |
372 * . If not, use size NWR or NWR+1 depending upon | |
373 * . which has the smaller corresponding subdiagonal | |
374 * . entry (a heuristic). ==== | |
375 * | |
376 NWINC = .TRUE. | |
377 IF( NH.LE.MIN( NMIN, NWMAX ) ) THEN | |
378 NW = NH | |
379 ELSE | |
380 NW = MIN( NWR, NH, NWMAX ) | |
381 IF( NW.LT.NWMAX ) THEN | |
382 IF( NW.GE.NH-1 ) THEN | |
383 NW = NH | |
384 ELSE | |
385 KWTOP = KBOT - NW + 1 | |
386 IF( ABS( H( KWTOP, KWTOP-1 ) ).GT. | |
387 $ ABS( H( KWTOP-1, KWTOP-2 ) ) )NW = NW + 1 | |
388 END IF | |
389 END IF | |
390 END IF | |
391 ELSE | |
392 * | |
393 * ==== Exceptional deflation window. If there have | |
394 * . been no deflations in KEXNW or more iterations, | |
395 * . then vary the deflation window size. At first, | |
396 * . because, larger windows are, in general, more | |
397 * . powerful than smaller ones, rapidly increase the | |
398 * . window up to the maximum reasonable and possible. | |
399 * . Then maybe try a slightly smaller window. ==== | |
400 * | |
401 IF( NWINC .AND. NW.LT.MIN( NWMAX, NH ) ) THEN | |
402 NW = MIN( NWMAX, NH, 2*NW ) | |
403 ELSE | |
404 NWINC = .FALSE. | |
405 IF( NW.EQ.NH .AND. NH.GT.2 ) | |
406 $ NW = NH - 1 | |
407 END IF | |
408 END IF | |
409 * | |
410 * ==== Aggressive early deflation: | |
411 * . split workspace under the subdiagonal into | |
412 * . - an nw-by-nw work array V in the lower | |
413 * . left-hand-corner, | |
414 * . - an NW-by-at-least-NW-but-more-is-better | |
415 * . (NW-by-NHO) horizontal work array along | |
416 * . the bottom edge, | |
417 * . - an at-least-NW-but-more-is-better (NHV-by-NW) | |
418 * . vertical work array along the left-hand-edge. | |
419 * . ==== | |
420 * | |
421 KV = N - NW + 1 | |
422 KT = NW + 1 | |
423 NHO = ( N-NW-1 ) - KT + 1 | |
424 KWV = NW + 2 | |
425 NVE = ( N-NW ) - KWV + 1 | |
426 * | |
427 * ==== Aggressive early deflation ==== | |
428 * | |
429 CALL DLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, | |
430 $ IHIZ, Z, LDZ, LS, LD, WR, WI, H( KV, 1 ), LDH, | |
431 $ NHO, H( KV, KT ), LDH, NVE, H( KWV, 1 ), LDH, | |
432 $ WORK, LWORK ) | |
433 * | |
434 * ==== Adjust KBOT accounting for new deflations. ==== | |
435 * | |
436 KBOT = KBOT - LD | |
437 * | |
438 * ==== KS points to the shifts. ==== | |
439 * | |
440 KS = KBOT - LS + 1 | |
441 * | |
442 * ==== Skip an expensive QR sweep if there is a (partly | |
443 * . heuristic) reason to expect that many eigenvalues | |
444 * . will deflate without it. Here, the QR sweep is | |
445 * . skipped if many eigenvalues have just been deflated | |
446 * . or if the remaining active block is small. | |
447 * | |
448 IF( ( LD.EQ.0 ) .OR. ( ( 100*LD.LE.NW*NIBBLE ) .AND. ( KBOT- | |
449 $ KTOP+1.GT.MIN( NMIN, NWMAX ) ) ) ) THEN | |
450 * | |
451 * ==== NS = nominal number of simultaneous shifts. | |
452 * . This may be lowered (slightly) if DLAQR2 | |
453 * . did not provide that many shifts. ==== | |
454 * | |
455 NS = MIN( NSMAX, NSR, MAX( 2, KBOT-KTOP ) ) | |
456 NS = NS - MOD( NS, 2 ) | |
457 * | |
458 * ==== If there have been no deflations | |
459 * . in a multiple of KEXSH iterations, | |
460 * . then try exceptional shifts. | |
461 * . Otherwise use shifts provided by | |
462 * . DLAQR2 above or from the eigenvalues | |
463 * . of a trailing principal submatrix. ==== | |
464 * | |
465 IF( MOD( NDFL, KEXSH ).EQ.0 ) THEN | |
466 KS = KBOT - NS + 1 | |
467 DO 30 I = KBOT, MAX( KS+1, KTOP+2 ), -2 | |
468 SS = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) ) | |
469 AA = WILK1*SS + H( I, I ) | |
470 BB = SS | |
471 CC = WILK2*SS | |
472 DD = AA | |
473 CALL DLANV2( AA, BB, CC, DD, WR( I-1 ), WI( I-1 ), | |
474 $ WR( I ), WI( I ), CS, SN ) | |
475 30 CONTINUE | |
476 IF( KS.EQ.KTOP ) THEN | |
477 WR( KS+1 ) = H( KS+1, KS+1 ) | |
478 WI( KS+1 ) = ZERO | |
479 WR( KS ) = WR( KS+1 ) | |
480 WI( KS ) = WI( KS+1 ) | |
481 END IF | |
482 ELSE | |
483 * | |
484 * ==== Got NS/2 or fewer shifts? Use DLAHQR | |
485 * . on a trailing principal submatrix to | |
486 * . get more. (Since NS.LE.NSMAX.LE.(N+6)/9, | |
487 * . there is enough space below the subdiagonal | |
488 * . to fit an NS-by-NS scratch array.) ==== | |
489 * | |
490 IF( KBOT-KS+1.LE.NS / 2 ) THEN | |
491 KS = KBOT - NS + 1 | |
492 KT = N - NS + 1 | |
493 CALL DLACPY( 'A', NS, NS, H( KS, KS ), LDH, | |
494 $ H( KT, 1 ), LDH ) | |
495 CALL DLAHQR( .false., .false., NS, 1, NS, | |
496 $ H( KT, 1 ), LDH, WR( KS ), WI( KS ), | |
497 $ 1, 1, ZDUM, 1, INF ) | |
498 KS = KS + INF | |
499 * | |
500 * ==== In case of a rare QR failure use | |
501 * . eigenvalues of the trailing 2-by-2 | |
502 * . principal submatrix. ==== | |
503 * | |
504 IF( KS.GE.KBOT ) THEN | |
505 AA = H( KBOT-1, KBOT-1 ) | |
506 CC = H( KBOT, KBOT-1 ) | |
507 BB = H( KBOT-1, KBOT ) | |
508 DD = H( KBOT, KBOT ) | |
509 CALL DLANV2( AA, BB, CC, DD, WR( KBOT-1 ), | |
510 $ WI( KBOT-1 ), WR( KBOT ), | |
511 $ WI( KBOT ), CS, SN ) | |
512 KS = KBOT - 1 | |
513 END IF | |
514 END IF | |
515 * | |
516 IF( KBOT-KS+1.GT.NS ) THEN | |
517 * | |
518 * ==== Sort the shifts (Helps a little) | |
519 * . Bubble sort keeps complex conjugate | |
520 * . pairs together. ==== | |
521 * | |
522 SORTED = .false. | |
523 DO 50 K = KBOT, KS + 1, -1 | |
524 IF( SORTED ) | |
525 $ GO TO 60 | |
526 SORTED = .true. | |
527 DO 40 I = KS, K - 1 | |
528 IF( ABS( WR( I ) )+ABS( WI( I ) ).LT. | |
529 $ ABS( WR( I+1 ) )+ABS( WI( I+1 ) ) ) THEN | |
530 SORTED = .false. | |
531 * | |
532 SWAP = WR( I ) | |
533 WR( I ) = WR( I+1 ) | |
534 WR( I+1 ) = SWAP | |
535 * | |
536 SWAP = WI( I ) | |
537 WI( I ) = WI( I+1 ) | |
538 WI( I+1 ) = SWAP | |
539 END IF | |
540 40 CONTINUE | |
541 50 CONTINUE | |
542 60 CONTINUE | |
543 END IF | |
544 * | |
545 * ==== Shuffle shifts into pairs of real shifts | |
546 * . and pairs of complex conjugate shifts | |
547 * . assuming complex conjugate shifts are | |
548 * . already adjacent to one another. (Yes, | |
549 * . they are.) ==== | |
550 * | |
551 DO 70 I = KBOT, KS + 2, -2 | |
552 IF( WI( I ).NE.-WI( I-1 ) ) THEN | |
553 * | |
554 SWAP = WR( I ) | |
555 WR( I ) = WR( I-1 ) | |
556 WR( I-1 ) = WR( I-2 ) | |
557 WR( I-2 ) = SWAP | |
558 * | |
559 SWAP = WI( I ) | |
560 WI( I ) = WI( I-1 ) | |
561 WI( I-1 ) = WI( I-2 ) | |
562 WI( I-2 ) = SWAP | |
563 END IF | |
564 70 CONTINUE | |
565 END IF | |
566 * | |
567 * ==== If there are only two shifts and both are | |
568 * . real, then use only one. ==== | |
569 * | |
570 IF( KBOT-KS+1.EQ.2 ) THEN | |
571 IF( WI( KBOT ).EQ.ZERO ) THEN | |
572 IF( ABS( WR( KBOT )-H( KBOT, KBOT ) ).LT. | |
573 $ ABS( WR( KBOT-1 )-H( KBOT, KBOT ) ) ) THEN | |
574 WR( KBOT-1 ) = WR( KBOT ) | |
575 ELSE | |
576 WR( KBOT ) = WR( KBOT-1 ) | |
577 END IF | |
578 END IF | |
579 END IF | |
580 * | |
581 * ==== Use up to NS of the the smallest magnatiude | |
582 * . shifts. If there aren't NS shifts available, | |
583 * . then use them all, possibly dropping one to | |
584 * . make the number of shifts even. ==== | |
585 * | |
586 NS = MIN( NS, KBOT-KS+1 ) | |
587 NS = NS - MOD( NS, 2 ) | |
588 KS = KBOT - NS + 1 | |
589 * | |
590 * ==== Small-bulge multi-shift QR sweep: | |
591 * . split workspace under the subdiagonal into | |
592 * . - a KDU-by-KDU work array U in the lower | |
593 * . left-hand-corner, | |
594 * . - a KDU-by-at-least-KDU-but-more-is-better | |
595 * . (KDU-by-NHo) horizontal work array WH along | |
596 * . the bottom edge, | |
597 * . - and an at-least-KDU-but-more-is-better-by-KDU | |
598 * . (NVE-by-KDU) vertical work WV arrow along | |
599 * . the left-hand-edge. ==== | |
600 * | |
601 KDU = 3*NS - 3 | |
602 KU = N - KDU + 1 | |
603 KWH = KDU + 1 | |
604 NHO = ( N-KDU+1-4 ) - ( KDU+1 ) + 1 | |
605 KWV = KDU + 4 | |
606 NVE = N - KDU - KWV + 1 | |
607 * | |
608 * ==== Small-bulge multi-shift QR sweep ==== | |
609 * | |
610 CALL DLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NS, | |
611 $ WR( KS ), WI( KS ), H, LDH, ILOZ, IHIZ, Z, | |
612 $ LDZ, WORK, 3, H( KU, 1 ), LDH, NVE, | |
613 $ H( KWV, 1 ), LDH, NHO, H( KU, KWH ), LDH ) | |
614 END IF | |
615 * | |
616 * ==== Note progress (or the lack of it). ==== | |
617 * | |
618 IF( LD.GT.0 ) THEN | |
619 NDFL = 1 | |
620 ELSE | |
621 NDFL = NDFL + 1 | |
622 END IF | |
623 * | |
624 * ==== End of main loop ==== | |
625 80 CONTINUE | |
626 * | |
627 * ==== Iteration limit exceeded. Set INFO to show where | |
628 * . the problem occurred and exit. ==== | |
629 * | |
630 INFO = KBOT | |
631 90 CONTINUE | |
632 END IF | |
633 * | |
634 * ==== Return the optimal value of LWORK. ==== | |
635 * | |
636 WORK( 1 ) = DBLE( LWKOPT ) | |
637 * | |
638 * ==== End of DLAQR4 ==== | |
639 * | |
640 END |