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