1 SUBROUTINE DHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z,
2 $ LDZ, WORK, LWORK, INFO )
3 *
4 * -- LAPACK computational routine (version 3.2.2) --
5 * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
6 * June 2010
7 *
8 * .. Scalar Arguments ..
9 INTEGER IHI, ILO, INFO, LDH, LDZ, LWORK, N
10 CHARACTER COMPZ, JOB
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION H( LDH, * ), WI( * ), WORK( * ), WR( * ),
14 $ Z( LDZ, * )
15 * ..
16 * Purpose
17 * =======
18 *
19 * DHSEQR computes the eigenvalues of a Hessenberg matrix H
20 * and, optionally, the matrices T and Z from the Schur decomposition
21 * H = Z T Z**T, where T is an upper quasi-triangular matrix (the
22 * Schur form), and Z is the orthogonal matrix of Schur vectors.
23 *
24 * Optionally Z may be postmultiplied into an input orthogonal
25 * matrix Q so that this routine can give the Schur factorization
26 * of a matrix A which has been reduced to the Hessenberg form H
27 * by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
28 *
29 * Arguments
30 * =========
31 *
32 * JOB (input) CHARACTER*1
33 * = 'E': compute eigenvalues only;
34 * = 'S': compute eigenvalues and the Schur form T.
35 *
36 * COMPZ (input) CHARACTER*1
37 * = 'N': no Schur vectors are computed;
38 * = 'I': Z is initialized to the unit matrix and the matrix Z
39 * of Schur vectors of H is returned;
40 * = 'V': Z must contain an orthogonal matrix Q on entry, and
41 * the product Q*Z is returned.
42 *
43 * N (input) INTEGER
44 * The order of the matrix H. N .GE. 0.
45 *
46 * ILO (input) INTEGER
47 * IHI (input) INTEGER
48 * It is assumed that H is already upper triangular in rows
49 * and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
50 * set by a previous call to DGEBAL, and then passed to DGEHRD
51 * when the matrix output by DGEBAL is reduced to Hessenberg
52 * form. Otherwise ILO and IHI should be set to 1 and N
53 * respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
54 * If N = 0, then ILO = 1 and IHI = 0.
55 *
56 * H (input/output) DOUBLE PRECISION array, dimension (LDH,N)
57 * On entry, the upper Hessenberg matrix H.
58 * On exit, if INFO = 0 and JOB = 'S', then H contains the
59 * upper quasi-triangular matrix T from the Schur decomposition
60 * (the Schur form); 2-by-2 diagonal blocks (corresponding to
61 * complex conjugate pairs of eigenvalues) are returned in
62 * standard form, with H(i,i) = H(i+1,i+1) and
63 * H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
64 * contents of H are unspecified on exit. (The output value of
65 * H when INFO.GT.0 is given under the description of INFO
66 * below.)
67 *
68 * Unlike earlier versions of DHSEQR, this subroutine may
69 * explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
70 * or j = IHI+1, IHI+2, ... N.
71 *
72 * LDH (input) INTEGER
73 * The leading dimension of the array H. LDH .GE. max(1,N).
74 *
75 * WR (output) DOUBLE PRECISION array, dimension (N)
76 * WI (output) DOUBLE PRECISION array, dimension (N)
77 * The real and imaginary parts, respectively, of the computed
78 * eigenvalues. If two eigenvalues are computed as a complex
79 * conjugate pair, they are stored in consecutive elements of
80 * WR and WI, say the i-th and (i+1)th, with WI(i) .GT. 0 and
81 * WI(i+1) .LT. 0. If JOB = 'S', the eigenvalues are stored in
82 * the same order as on the diagonal of the Schur form returned
83 * in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2
84 * diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
85 * WI(i+1) = -WI(i).
86 *
87 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
88 * If COMPZ = 'N', Z is not referenced.
89 * If COMPZ = 'I', on entry Z need not be set and on exit,
90 * if INFO = 0, Z contains the orthogonal matrix Z of the Schur
91 * vectors of H. If COMPZ = 'V', on entry Z must contain an
92 * N-by-N matrix Q, which is assumed to be equal to the unit
93 * matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
94 * if INFO = 0, Z contains Q*Z.
95 * Normally Q is the orthogonal matrix generated by DORGHR
96 * after the call to DGEHRD which formed the Hessenberg matrix
97 * H. (The output value of Z when INFO.GT.0 is given under
98 * the description of INFO below.)
99 *
100 * LDZ (input) INTEGER
101 * The leading dimension of the array Z. if COMPZ = 'I' or
102 * COMPZ = 'V', then LDZ.GE.MAX(1,N). Otherwize, LDZ.GE.1.
103 *
104 * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
105 * On exit, if INFO = 0, WORK(1) returns an estimate of
106 * the optimal value for LWORK.
107 *
108 * LWORK (input) INTEGER
109 * The dimension of the array WORK. LWORK .GE. max(1,N)
110 * is sufficient and delivers very good and sometimes
111 * optimal performance. However, LWORK as large as 11*N
112 * may be required for optimal performance. A workspace
113 * query is recommended to determine the optimal workspace
114 * size.
115 *
116 * If LWORK = -1, then DHSEQR does a workspace query.
117 * In this case, DHSEQR checks the input parameters and
118 * estimates the optimal workspace size for the given
119 * values of N, ILO and IHI. The estimate is returned
120 * in WORK(1). No error message related to LWORK is
121 * issued by XERBLA. Neither H nor Z are accessed.
122 *
123 *
124 * INFO (output) INTEGER
125 * = 0: successful exit
126 * .LT. 0: if INFO = -i, the i-th argument had an illegal
127 * value
128 * .GT. 0: if INFO = i, DHSEQR failed to compute all of
129 * the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
130 * and WI contain those eigenvalues which have been
131 * successfully computed. (Failures are rare.)
132 *
133 * If INFO .GT. 0 and JOB = 'E', then on exit, the
134 * remaining unconverged eigenvalues are the eigen-
135 * values of the upper Hessenberg matrix rows and
136 * columns ILO through INFO of the final, output
137 * value of H.
138 *
139 * If INFO .GT. 0 and JOB = 'S', then on exit
140 *
141 * (*) (initial value of H)*U = U*(final value of H)
142 *
143 * where U is an orthogonal matrix. The final
144 * value of H is upper Hessenberg and quasi-triangular
145 * in rows and columns INFO+1 through IHI.
146 *
147 * If INFO .GT. 0 and COMPZ = 'V', then on exit
148 *
149 * (final value of Z) = (initial value of Z)*U
150 *
151 * where U is the orthogonal matrix in (*) (regard-
152 * less of the value of JOB.)
153 *
154 * If INFO .GT. 0 and COMPZ = 'I', then on exit
155 * (final value of Z) = U
156 * where U is the orthogonal matrix in (*) (regard-
157 * less of the value of JOB.)
158 *
159 * If INFO .GT. 0 and COMPZ = 'N', then Z is not
160 * accessed.
161 *
162 * ================================================================
163 * Default values supplied by
164 * ILAENV(ISPEC,'DHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
165 * It is suggested that these defaults be adjusted in order
166 * to attain best performance in each particular
167 * computational environment.
168 *
169 * ISPEC=12: The DLAHQR vs DLAQR0 crossover point.
170 * Default: 75. (Must be at least 11.)
171 *
172 * ISPEC=13: Recommended deflation window size.
173 * This depends on ILO, IHI and NS. NS is the
174 * number of simultaneous shifts returned
175 * by ILAENV(ISPEC=15). (See ISPEC=15 below.)
176 * The default for (IHI-ILO+1).LE.500 is NS.
177 * The default for (IHI-ILO+1).GT.500 is 3*NS/2.
178 *
179 * ISPEC=14: Nibble crossover point. (See IPARMQ for
180 * details.) Default: 14% of deflation window
181 * size.
182 *
183 * ISPEC=15: Number of simultaneous shifts in a multishift
184 * QR iteration.
185 *
186 * If IHI-ILO+1 is ...
187 *
188 * greater than ...but less ... the
189 * or equal to ... than default is
190 *
191 * 1 30 NS = 2(+)
192 * 30 60 NS = 4(+)
193 * 60 150 NS = 10(+)
194 * 150 590 NS = **
195 * 590 3000 NS = 64
196 * 3000 6000 NS = 128
197 * 6000 infinity NS = 256
198 *
199 * (+) By default some or all matrices of this order
200 * are passed to the implicit double shift routine
201 * DLAHQR and this parameter is ignored. See
202 * ISPEC=12 above and comments in IPARMQ for
203 * details.
204 *
205 * (**) The asterisks (**) indicate an ad-hoc
206 * function of N increasing from 10 to 64.
207 *
208 * ISPEC=16: Select structured matrix multiply.
209 * If the number of simultaneous shifts (specified
210 * by ISPEC=15) is less than 14, then the default
211 * for ISPEC=16 is 0. Otherwise the default for
212 * ISPEC=16 is 2.
213 *
214 * ================================================================
215 * Based on contributions by
216 * Karen Braman and Ralph Byers, Department of Mathematics,
217 * University of Kansas, USA
218 *
219 * ================================================================
220 * References:
221 * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
222 * Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
223 * Performance, SIAM Journal of Matrix Analysis, volume 23, pages
224 * 929--947, 2002.
225 *
226 * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
227 * Algorithm Part II: Aggressive Early Deflation, SIAM Journal
228 * of Matrix Analysis, volume 23, pages 948--973, 2002.
229 *
230 * ================================================================
231 * .. Parameters ..
232 *
233 * ==== Matrices of order NTINY or smaller must be processed by
234 * . DLAHQR because of insufficient subdiagonal scratch space.
235 * . (This is a hard limit.) ====
236 INTEGER NTINY
237 PARAMETER ( NTINY = 11 )
238 *
239 * ==== NL allocates some local workspace to help small matrices
240 * . through a rare DLAHQR failure. NL .GT. NTINY = 11 is
241 * . required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom-
242 * . mended. (The default value of NMIN is 75.) Using NL = 49
243 * . allows up to six simultaneous shifts and a 16-by-16
244 * . deflation window. ====
245 INTEGER NL
246 PARAMETER ( NL = 49 )
247 DOUBLE PRECISION ZERO, ONE
248 PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )
249 * ..
250 * .. Local Arrays ..
251 DOUBLE PRECISION HL( NL, NL ), WORKL( NL )
252 * ..
253 * .. Local Scalars ..
254 INTEGER I, KBOT, NMIN
255 LOGICAL INITZ, LQUERY, WANTT, WANTZ
256 * ..
257 * .. External Functions ..
258 INTEGER ILAENV
259 LOGICAL LSAME
260 EXTERNAL ILAENV, LSAME
261 * ..
262 * .. External Subroutines ..
263 EXTERNAL DLACPY, DLAHQR, DLAQR0, DLASET, XERBLA
264 * ..
265 * .. Intrinsic Functions ..
266 INTRINSIC DBLE, MAX, MIN
267 * ..
268 * .. Executable Statements ..
269 *
270 * ==== Decode and check the input parameters. ====
271 *
272 WANTT = LSAME( JOB, 'S' )
273 INITZ = LSAME( COMPZ, 'I' )
274 WANTZ = INITZ .OR. LSAME( COMPZ, 'V' )
275 WORK( 1 ) = DBLE( MAX( 1, N ) )
276 LQUERY = LWORK.EQ.-1
277 *
278 INFO = 0
279 IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN
280 INFO = -1
281 ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN
282 INFO = -2
283 ELSE IF( N.LT.0 ) THEN
284 INFO = -3
285 ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
286 INFO = -4
287 ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
288 INFO = -5
289 ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
290 INFO = -7
291 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
292 INFO = -11
293 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
294 INFO = -13
295 END IF
296 *
297 IF( INFO.NE.0 ) THEN
298 *
299 * ==== Quick return in case of invalid argument. ====
300 *
301 CALL XERBLA( 'DHSEQR', -INFO )
302 RETURN
303 *
304 ELSE IF( N.EQ.0 ) THEN
305 *
306 * ==== Quick return in case N = 0; nothing to do. ====
307 *
308 RETURN
309 *
310 ELSE IF( LQUERY ) THEN
311 *
312 * ==== Quick return in case of a workspace query ====
313 *
314 CALL DLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
315 $ IHI, Z, LDZ, WORK, LWORK, INFO )
316 * ==== Ensure reported workspace size is backward-compatible with
317 * . previous LAPACK versions. ====
318 WORK( 1 ) = MAX( DBLE( MAX( 1, N ) ), WORK( 1 ) )
319 RETURN
320 *
321 ELSE
322 *
323 * ==== copy eigenvalues isolated by DGEBAL ====
324 *
325 DO 10 I = 1, ILO - 1
326 WR( I ) = H( I, I )
327 WI( I ) = ZERO
328 10 CONTINUE
329 DO 20 I = IHI + 1, N
330 WR( I ) = H( I, I )
331 WI( I ) = ZERO
332 20 CONTINUE
333 *
334 * ==== Initialize Z, if requested ====
335 *
336 IF( INITZ )
337 $ CALL DLASET( 'A', N, N, ZERO, ONE, Z, LDZ )
338 *
339 * ==== Quick return if possible ====
340 *
341 IF( ILO.EQ.IHI ) THEN
342 WR( ILO ) = H( ILO, ILO )
343 WI( ILO ) = ZERO
344 RETURN
345 END IF
346 *
347 * ==== DLAHQR/DLAQR0 crossover point ====
348 *
349 NMIN = ILAENV( 12, 'DHSEQR', JOB( : 1 ) // COMPZ( : 1 ), N,
350 $ ILO, IHI, LWORK )
351 NMIN = MAX( NTINY, NMIN )
352 *
353 * ==== DLAQR0 for big matrices; DLAHQR for small ones ====
354 *
355 IF( N.GT.NMIN ) THEN
356 CALL DLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
357 $ IHI, Z, LDZ, WORK, LWORK, INFO )
358 ELSE
359 *
360 * ==== Small matrix ====
361 *
362 CALL DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
363 $ IHI, Z, LDZ, INFO )
364 *
365 IF( INFO.GT.0 ) THEN
366 *
367 * ==== A rare DLAHQR failure! DLAQR0 sometimes succeeds
368 * . when DLAHQR fails. ====
369 *
370 KBOT = INFO
371 *
372 IF( N.GE.NL ) THEN
373 *
374 * ==== Larger matrices have enough subdiagonal scratch
375 * . space to call DLAQR0 directly. ====
376 *
377 CALL DLAQR0( WANTT, WANTZ, N, ILO, KBOT, H, LDH, WR,
378 $ WI, ILO, IHI, Z, LDZ, WORK, LWORK, INFO )
379 *
380 ELSE
381 *
382 * ==== Tiny matrices don't have enough subdiagonal
383 * . scratch space to benefit from DLAQR0. Hence,
384 * . tiny matrices must be copied into a larger
385 * . array before calling DLAQR0. ====
386 *
387 CALL DLACPY( 'A', N, N, H, LDH, HL, NL )
388 HL( N+1, N ) = ZERO
389 CALL DLASET( 'A', NL, NL-N, ZERO, ZERO, HL( 1, N+1 ),
390 $ NL )
391 CALL DLAQR0( WANTT, WANTZ, NL, ILO, KBOT, HL, NL, WR,
392 $ WI, ILO, IHI, Z, LDZ, WORKL, NL, INFO )
393 IF( WANTT .OR. INFO.NE.0 )
394 $ CALL DLACPY( 'A', N, N, HL, NL, H, LDH )
395 END IF
396 END IF
397 END IF
398 *
399 * ==== Clear out the trash, if necessary. ====
400 *
401 IF( ( WANTT .OR. INFO.NE.0 ) .AND. N.GT.2 )
402 $ CALL DLASET( 'L', N-2, N-2, ZERO, ZERO, H( 3, 1 ), LDH )
403 *
404 * ==== Ensure reported workspace size is backward-compatible with
405 * . previous LAPACK versions. ====
406 *
407 WORK( 1 ) = MAX( DBLE( MAX( 1, N ) ), WORK( 1 ) )
408 END IF
409 *
410 * ==== End of DHSEQR ====
411 *
412 END
2 $ LDZ, WORK, LWORK, INFO )
3 *
4 * -- LAPACK computational routine (version 3.2.2) --
5 * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
6 * June 2010
7 *
8 * .. Scalar Arguments ..
9 INTEGER IHI, ILO, INFO, LDH, LDZ, LWORK, N
10 CHARACTER COMPZ, JOB
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION H( LDH, * ), WI( * ), WORK( * ), WR( * ),
14 $ Z( LDZ, * )
15 * ..
16 * Purpose
17 * =======
18 *
19 * DHSEQR computes the eigenvalues of a Hessenberg matrix H
20 * and, optionally, the matrices T and Z from the Schur decomposition
21 * H = Z T Z**T, where T is an upper quasi-triangular matrix (the
22 * Schur form), and Z is the orthogonal matrix of Schur vectors.
23 *
24 * Optionally Z may be postmultiplied into an input orthogonal
25 * matrix Q so that this routine can give the Schur factorization
26 * of a matrix A which has been reduced to the Hessenberg form H
27 * by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
28 *
29 * Arguments
30 * =========
31 *
32 * JOB (input) CHARACTER*1
33 * = 'E': compute eigenvalues only;
34 * = 'S': compute eigenvalues and the Schur form T.
35 *
36 * COMPZ (input) CHARACTER*1
37 * = 'N': no Schur vectors are computed;
38 * = 'I': Z is initialized to the unit matrix and the matrix Z
39 * of Schur vectors of H is returned;
40 * = 'V': Z must contain an orthogonal matrix Q on entry, and
41 * the product Q*Z is returned.
42 *
43 * N (input) INTEGER
44 * The order of the matrix H. N .GE. 0.
45 *
46 * ILO (input) INTEGER
47 * IHI (input) INTEGER
48 * It is assumed that H is already upper triangular in rows
49 * and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
50 * set by a previous call to DGEBAL, and then passed to DGEHRD
51 * when the matrix output by DGEBAL is reduced to Hessenberg
52 * form. Otherwise ILO and IHI should be set to 1 and N
53 * respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
54 * If N = 0, then ILO = 1 and IHI = 0.
55 *
56 * H (input/output) DOUBLE PRECISION array, dimension (LDH,N)
57 * On entry, the upper Hessenberg matrix H.
58 * On exit, if INFO = 0 and JOB = 'S', then H contains the
59 * upper quasi-triangular matrix T from the Schur decomposition
60 * (the Schur form); 2-by-2 diagonal blocks (corresponding to
61 * complex conjugate pairs of eigenvalues) are returned in
62 * standard form, with H(i,i) = H(i+1,i+1) and
63 * H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
64 * contents of H are unspecified on exit. (The output value of
65 * H when INFO.GT.0 is given under the description of INFO
66 * below.)
67 *
68 * Unlike earlier versions of DHSEQR, this subroutine may
69 * explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
70 * or j = IHI+1, IHI+2, ... N.
71 *
72 * LDH (input) INTEGER
73 * The leading dimension of the array H. LDH .GE. max(1,N).
74 *
75 * WR (output) DOUBLE PRECISION array, dimension (N)
76 * WI (output) DOUBLE PRECISION array, dimension (N)
77 * The real and imaginary parts, respectively, of the computed
78 * eigenvalues. If two eigenvalues are computed as a complex
79 * conjugate pair, they are stored in consecutive elements of
80 * WR and WI, say the i-th and (i+1)th, with WI(i) .GT. 0 and
81 * WI(i+1) .LT. 0. If JOB = 'S', the eigenvalues are stored in
82 * the same order as on the diagonal of the Schur form returned
83 * in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2
84 * diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
85 * WI(i+1) = -WI(i).
86 *
87 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
88 * If COMPZ = 'N', Z is not referenced.
89 * If COMPZ = 'I', on entry Z need not be set and on exit,
90 * if INFO = 0, Z contains the orthogonal matrix Z of the Schur
91 * vectors of H. If COMPZ = 'V', on entry Z must contain an
92 * N-by-N matrix Q, which is assumed to be equal to the unit
93 * matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
94 * if INFO = 0, Z contains Q*Z.
95 * Normally Q is the orthogonal matrix generated by DORGHR
96 * after the call to DGEHRD which formed the Hessenberg matrix
97 * H. (The output value of Z when INFO.GT.0 is given under
98 * the description of INFO below.)
99 *
100 * LDZ (input) INTEGER
101 * The leading dimension of the array Z. if COMPZ = 'I' or
102 * COMPZ = 'V', then LDZ.GE.MAX(1,N). Otherwize, LDZ.GE.1.
103 *
104 * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
105 * On exit, if INFO = 0, WORK(1) returns an estimate of
106 * the optimal value for LWORK.
107 *
108 * LWORK (input) INTEGER
109 * The dimension of the array WORK. LWORK .GE. max(1,N)
110 * is sufficient and delivers very good and sometimes
111 * optimal performance. However, LWORK as large as 11*N
112 * may be required for optimal performance. A workspace
113 * query is recommended to determine the optimal workspace
114 * size.
115 *
116 * If LWORK = -1, then DHSEQR does a workspace query.
117 * In this case, DHSEQR checks the input parameters and
118 * estimates the optimal workspace size for the given
119 * values of N, ILO and IHI. The estimate is returned
120 * in WORK(1). No error message related to LWORK is
121 * issued by XERBLA. Neither H nor Z are accessed.
122 *
123 *
124 * INFO (output) INTEGER
125 * = 0: successful exit
126 * .LT. 0: if INFO = -i, the i-th argument had an illegal
127 * value
128 * .GT. 0: if INFO = i, DHSEQR failed to compute all of
129 * the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
130 * and WI contain those eigenvalues which have been
131 * successfully computed. (Failures are rare.)
132 *
133 * If INFO .GT. 0 and JOB = 'E', then on exit, the
134 * remaining unconverged eigenvalues are the eigen-
135 * values of the upper Hessenberg matrix rows and
136 * columns ILO through INFO of the final, output
137 * value of H.
138 *
139 * If INFO .GT. 0 and JOB = 'S', then on exit
140 *
141 * (*) (initial value of H)*U = U*(final value of H)
142 *
143 * where U is an orthogonal matrix. The final
144 * value of H is upper Hessenberg and quasi-triangular
145 * in rows and columns INFO+1 through IHI.
146 *
147 * If INFO .GT. 0 and COMPZ = 'V', then on exit
148 *
149 * (final value of Z) = (initial value of Z)*U
150 *
151 * where U is the orthogonal matrix in (*) (regard-
152 * less of the value of JOB.)
153 *
154 * If INFO .GT. 0 and COMPZ = 'I', then on exit
155 * (final value of Z) = U
156 * where U is the orthogonal matrix in (*) (regard-
157 * less of the value of JOB.)
158 *
159 * If INFO .GT. 0 and COMPZ = 'N', then Z is not
160 * accessed.
161 *
162 * ================================================================
163 * Default values supplied by
164 * ILAENV(ISPEC,'DHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
165 * It is suggested that these defaults be adjusted in order
166 * to attain best performance in each particular
167 * computational environment.
168 *
169 * ISPEC=12: The DLAHQR vs DLAQR0 crossover point.
170 * Default: 75. (Must be at least 11.)
171 *
172 * ISPEC=13: Recommended deflation window size.
173 * This depends on ILO, IHI and NS. NS is the
174 * number of simultaneous shifts returned
175 * by ILAENV(ISPEC=15). (See ISPEC=15 below.)
176 * The default for (IHI-ILO+1).LE.500 is NS.
177 * The default for (IHI-ILO+1).GT.500 is 3*NS/2.
178 *
179 * ISPEC=14: Nibble crossover point. (See IPARMQ for
180 * details.) Default: 14% of deflation window
181 * size.
182 *
183 * ISPEC=15: Number of simultaneous shifts in a multishift
184 * QR iteration.
185 *
186 * If IHI-ILO+1 is ...
187 *
188 * greater than ...but less ... the
189 * or equal to ... than default is
190 *
191 * 1 30 NS = 2(+)
192 * 30 60 NS = 4(+)
193 * 60 150 NS = 10(+)
194 * 150 590 NS = **
195 * 590 3000 NS = 64
196 * 3000 6000 NS = 128
197 * 6000 infinity NS = 256
198 *
199 * (+) By default some or all matrices of this order
200 * are passed to the implicit double shift routine
201 * DLAHQR and this parameter is ignored. See
202 * ISPEC=12 above and comments in IPARMQ for
203 * details.
204 *
205 * (**) The asterisks (**) indicate an ad-hoc
206 * function of N increasing from 10 to 64.
207 *
208 * ISPEC=16: Select structured matrix multiply.
209 * If the number of simultaneous shifts (specified
210 * by ISPEC=15) is less than 14, then the default
211 * for ISPEC=16 is 0. Otherwise the default for
212 * ISPEC=16 is 2.
213 *
214 * ================================================================
215 * Based on contributions by
216 * Karen Braman and Ralph Byers, Department of Mathematics,
217 * University of Kansas, USA
218 *
219 * ================================================================
220 * References:
221 * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
222 * Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
223 * Performance, SIAM Journal of Matrix Analysis, volume 23, pages
224 * 929--947, 2002.
225 *
226 * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
227 * Algorithm Part II: Aggressive Early Deflation, SIAM Journal
228 * of Matrix Analysis, volume 23, pages 948--973, 2002.
229 *
230 * ================================================================
231 * .. Parameters ..
232 *
233 * ==== Matrices of order NTINY or smaller must be processed by
234 * . DLAHQR because of insufficient subdiagonal scratch space.
235 * . (This is a hard limit.) ====
236 INTEGER NTINY
237 PARAMETER ( NTINY = 11 )
238 *
239 * ==== NL allocates some local workspace to help small matrices
240 * . through a rare DLAHQR failure. NL .GT. NTINY = 11 is
241 * . required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom-
242 * . mended. (The default value of NMIN is 75.) Using NL = 49
243 * . allows up to six simultaneous shifts and a 16-by-16
244 * . deflation window. ====
245 INTEGER NL
246 PARAMETER ( NL = 49 )
247 DOUBLE PRECISION ZERO, ONE
248 PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )
249 * ..
250 * .. Local Arrays ..
251 DOUBLE PRECISION HL( NL, NL ), WORKL( NL )
252 * ..
253 * .. Local Scalars ..
254 INTEGER I, KBOT, NMIN
255 LOGICAL INITZ, LQUERY, WANTT, WANTZ
256 * ..
257 * .. External Functions ..
258 INTEGER ILAENV
259 LOGICAL LSAME
260 EXTERNAL ILAENV, LSAME
261 * ..
262 * .. External Subroutines ..
263 EXTERNAL DLACPY, DLAHQR, DLAQR0, DLASET, XERBLA
264 * ..
265 * .. Intrinsic Functions ..
266 INTRINSIC DBLE, MAX, MIN
267 * ..
268 * .. Executable Statements ..
269 *
270 * ==== Decode and check the input parameters. ====
271 *
272 WANTT = LSAME( JOB, 'S' )
273 INITZ = LSAME( COMPZ, 'I' )
274 WANTZ = INITZ .OR. LSAME( COMPZ, 'V' )
275 WORK( 1 ) = DBLE( MAX( 1, N ) )
276 LQUERY = LWORK.EQ.-1
277 *
278 INFO = 0
279 IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN
280 INFO = -1
281 ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN
282 INFO = -2
283 ELSE IF( N.LT.0 ) THEN
284 INFO = -3
285 ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
286 INFO = -4
287 ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
288 INFO = -5
289 ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
290 INFO = -7
291 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
292 INFO = -11
293 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
294 INFO = -13
295 END IF
296 *
297 IF( INFO.NE.0 ) THEN
298 *
299 * ==== Quick return in case of invalid argument. ====
300 *
301 CALL XERBLA( 'DHSEQR', -INFO )
302 RETURN
303 *
304 ELSE IF( N.EQ.0 ) THEN
305 *
306 * ==== Quick return in case N = 0; nothing to do. ====
307 *
308 RETURN
309 *
310 ELSE IF( LQUERY ) THEN
311 *
312 * ==== Quick return in case of a workspace query ====
313 *
314 CALL DLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
315 $ IHI, Z, LDZ, WORK, LWORK, INFO )
316 * ==== Ensure reported workspace size is backward-compatible with
317 * . previous LAPACK versions. ====
318 WORK( 1 ) = MAX( DBLE( MAX( 1, N ) ), WORK( 1 ) )
319 RETURN
320 *
321 ELSE
322 *
323 * ==== copy eigenvalues isolated by DGEBAL ====
324 *
325 DO 10 I = 1, ILO - 1
326 WR( I ) = H( I, I )
327 WI( I ) = ZERO
328 10 CONTINUE
329 DO 20 I = IHI + 1, N
330 WR( I ) = H( I, I )
331 WI( I ) = ZERO
332 20 CONTINUE
333 *
334 * ==== Initialize Z, if requested ====
335 *
336 IF( INITZ )
337 $ CALL DLASET( 'A', N, N, ZERO, ONE, Z, LDZ )
338 *
339 * ==== Quick return if possible ====
340 *
341 IF( ILO.EQ.IHI ) THEN
342 WR( ILO ) = H( ILO, ILO )
343 WI( ILO ) = ZERO
344 RETURN
345 END IF
346 *
347 * ==== DLAHQR/DLAQR0 crossover point ====
348 *
349 NMIN = ILAENV( 12, 'DHSEQR', JOB( : 1 ) // COMPZ( : 1 ), N,
350 $ ILO, IHI, LWORK )
351 NMIN = MAX( NTINY, NMIN )
352 *
353 * ==== DLAQR0 for big matrices; DLAHQR for small ones ====
354 *
355 IF( N.GT.NMIN ) THEN
356 CALL DLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
357 $ IHI, Z, LDZ, WORK, LWORK, INFO )
358 ELSE
359 *
360 * ==== Small matrix ====
361 *
362 CALL DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
363 $ IHI, Z, LDZ, INFO )
364 *
365 IF( INFO.GT.0 ) THEN
366 *
367 * ==== A rare DLAHQR failure! DLAQR0 sometimes succeeds
368 * . when DLAHQR fails. ====
369 *
370 KBOT = INFO
371 *
372 IF( N.GE.NL ) THEN
373 *
374 * ==== Larger matrices have enough subdiagonal scratch
375 * . space to call DLAQR0 directly. ====
376 *
377 CALL DLAQR0( WANTT, WANTZ, N, ILO, KBOT, H, LDH, WR,
378 $ WI, ILO, IHI, Z, LDZ, WORK, LWORK, INFO )
379 *
380 ELSE
381 *
382 * ==== Tiny matrices don't have enough subdiagonal
383 * . scratch space to benefit from DLAQR0. Hence,
384 * . tiny matrices must be copied into a larger
385 * . array before calling DLAQR0. ====
386 *
387 CALL DLACPY( 'A', N, N, H, LDH, HL, NL )
388 HL( N+1, N ) = ZERO
389 CALL DLASET( 'A', NL, NL-N, ZERO, ZERO, HL( 1, N+1 ),
390 $ NL )
391 CALL DLAQR0( WANTT, WANTZ, NL, ILO, KBOT, HL, NL, WR,
392 $ WI, ILO, IHI, Z, LDZ, WORKL, NL, INFO )
393 IF( WANTT .OR. INFO.NE.0 )
394 $ CALL DLACPY( 'A', N, N, HL, NL, H, LDH )
395 END IF
396 END IF
397 END IF
398 *
399 * ==== Clear out the trash, if necessary. ====
400 *
401 IF( ( WANTT .OR. INFO.NE.0 ) .AND. N.GT.2 )
402 $ CALL DLASET( 'L', N-2, N-2, ZERO, ZERO, H( 3, 1 ), LDH )
403 *
404 * ==== Ensure reported workspace size is backward-compatible with
405 * . previous LAPACK versions. ====
406 *
407 WORK( 1 ) = MAX( DBLE( MAX( 1, N ) ), WORK( 1 ) )
408 END IF
409 *
410 * ==== End of DHSEQR ====
411 *
412 END