LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
zhseqr.f
Go to the documentation of this file.
1 *> \brief \b ZHSEQR
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZHSEQR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhseqr.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhseqr.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhseqr.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ,
22 * WORK, LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER IHI, ILO, INFO, LDH, LDZ, LWORK, N
26 * CHARACTER COMPZ, JOB
27 * ..
28 * .. Array Arguments ..
29 * COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> ZHSEQR computes the eigenvalues of a Hessenberg matrix H
39 *> and, optionally, the matrices T and Z from the Schur decomposition
40 *> H = Z T Z**H, where T is an upper triangular matrix (the
41 *> Schur form), and Z is the unitary matrix of Schur vectors.
42 *>
43 *> Optionally Z may be postmultiplied into an input unitary
44 *> matrix Q so that this routine can give the Schur factorization
45 *> of a matrix A which has been reduced to the Hessenberg form H
46 *> by the unitary matrix Q: A = Q*H*Q**H = (QZ)*T*(QZ)**H.
47 *> \endverbatim
48 *
49 * Arguments:
50 * ==========
51 *
52 *> \param[in] JOB
53 *> \verbatim
54 *> JOB is CHARACTER*1
55 *> = 'E': compute eigenvalues only;
56 *> = 'S': compute eigenvalues and the Schur form T.
57 *> \endverbatim
58 *>
59 *> \param[in] COMPZ
60 *> \verbatim
61 *> COMPZ is CHARACTER*1
62 *> = 'N': no Schur vectors are computed;
63 *> = 'I': Z is initialized to the unit matrix and the matrix Z
64 *> of Schur vectors of H is returned;
65 *> = 'V': Z must contain an unitary matrix Q on entry, and
66 *> the product Q*Z is returned.
67 *> \endverbatim
68 *>
69 *> \param[in] N
70 *> \verbatim
71 *> N is INTEGER
72 *> The order of the matrix H. N >= 0.
73 *> \endverbatim
74 *>
75 *> \param[in] ILO
76 *> \verbatim
77 *> ILO is INTEGER
78 *> \endverbatim
79 *>
80 *> \param[in] IHI
81 *> \verbatim
82 *> IHI is INTEGER
83 *>
84 *> It is assumed that H is already upper triangular in rows
85 *> and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
86 *> set by a previous call to ZGEBAL, and then passed to ZGEHRD
87 *> when the matrix output by ZGEBAL is reduced to Hessenberg
88 *> form. Otherwise ILO and IHI should be set to 1 and N
89 *> respectively. If N > 0, then 1 <= ILO <= IHI <= N.
90 *> If N = 0, then ILO = 1 and IHI = 0.
91 *> \endverbatim
92 *>
93 *> \param[in,out] H
94 *> \verbatim
95 *> H is COMPLEX*16 array, dimension (LDH,N)
96 *> On entry, the upper Hessenberg matrix H.
97 *> On exit, if INFO = 0 and JOB = 'S', H contains the upper
98 *> triangular matrix T from the Schur decomposition (the
99 *> Schur form). If INFO = 0 and JOB = 'E', the contents of
100 *> H are unspecified on exit. (The output value of H when
101 *> INFO > 0 is given under the description of INFO below.)
102 *>
103 *> Unlike earlier versions of ZHSEQR, this subroutine may
104 *> explicitly H(i,j) = 0 for i > j and j = 1, 2, ... ILO-1
105 *> or j = IHI+1, IHI+2, ... N.
106 *> \endverbatim
107 *>
108 *> \param[in] LDH
109 *> \verbatim
110 *> LDH is INTEGER
111 *> The leading dimension of the array H. LDH >= max(1,N).
112 *> \endverbatim
113 *>
114 *> \param[out] W
115 *> \verbatim
116 *> W is COMPLEX*16 array, dimension (N)
117 *> The computed eigenvalues. If JOB = 'S', the eigenvalues are
118 *> stored in the same order as on the diagonal of the Schur
119 *> form returned in H, with W(i) = H(i,i).
120 *> \endverbatim
121 *>
122 *> \param[in,out] Z
123 *> \verbatim
124 *> Z is COMPLEX*16 array, dimension (LDZ,N)
125 *> If COMPZ = 'N', Z is not referenced.
126 *> If COMPZ = 'I', on entry Z need not be set and on exit,
127 *> if INFO = 0, Z contains the unitary matrix Z of the Schur
128 *> vectors of H. If COMPZ = 'V', on entry Z must contain an
129 *> N-by-N matrix Q, which is assumed to be equal to the unit
130 *> matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
131 *> if INFO = 0, Z contains Q*Z.
132 *> Normally Q is the unitary matrix generated by ZUNGHR
133 *> after the call to ZGEHRD which formed the Hessenberg matrix
134 *> H. (The output value of Z when INFO > 0 is given under
135 *> the description of INFO below.)
136 *> \endverbatim
137 *>
138 *> \param[in] LDZ
139 *> \verbatim
140 *> LDZ is INTEGER
141 *> The leading dimension of the array Z. if COMPZ = 'I' or
142 *> COMPZ = 'V', then LDZ >= MAX(1,N). Otherwise, LDZ >= 1.
143 *> \endverbatim
144 *>
145 *> \param[out] WORK
146 *> \verbatim
147 *> WORK is COMPLEX*16 array, dimension (LWORK)
148 *> On exit, if INFO = 0, WORK(1) returns an estimate of
149 *> the optimal value for LWORK.
150 *> \endverbatim
151 *>
152 *> \param[in] LWORK
153 *> \verbatim
154 *> LWORK is INTEGER
155 *> The dimension of the array WORK. LWORK >= max(1,N)
156 *> is sufficient and delivers very good and sometimes
157 *> optimal performance. However, LWORK as large as 11*N
158 *> may be required for optimal performance. A workspace
159 *> query is recommended to determine the optimal workspace
160 *> size.
161 *>
162 *> If LWORK = -1, then ZHSEQR does a workspace query.
163 *> In this case, ZHSEQR checks the input parameters and
164 *> estimates the optimal workspace size for the given
165 *> values of N, ILO and IHI. The estimate is returned
166 *> in WORK(1). No error message related to LWORK is
167 *> issued by XERBLA. Neither H nor Z are accessed.
168 *> \endverbatim
169 *>
170 *> \param[out] INFO
171 *> \verbatim
172 *> INFO is INTEGER
173 *> = 0: successful exit
174 *> < 0: if INFO = -i, the i-th argument had an illegal
175 *> value
176 *> > 0: if INFO = i, ZHSEQR failed to compute all of
177 *> the eigenvalues. Elements 1:ilo-1 and i+1:n of W
178 *> contain those eigenvalues which have been
179 *> successfully computed. (Failures are rare.)
180 *>
181 *> If INFO > 0 and JOB = 'E', then on exit, the
182 *> remaining unconverged eigenvalues are the eigen-
183 *> values of the upper Hessenberg matrix rows and
184 *> columns ILO through INFO of the final, output
185 *> value of H.
186 *>
187 *> If INFO > 0 and JOB = 'S', then on exit
188 *>
189 *> (*) (initial value of H)*U = U*(final value of H)
190 *>
191 *> where U is a unitary matrix. The final
192 *> value of H is upper Hessenberg and triangular in
193 *> rows and columns INFO+1 through IHI.
194 *>
195 *> If INFO > 0 and COMPZ = 'V', then on exit
196 *>
197 *> (final value of Z) = (initial value of Z)*U
198 *>
199 *> where U is the unitary matrix in (*) (regard-
200 *> less of the value of JOB.)
201 *>
202 *> If INFO > 0 and COMPZ = 'I', then on exit
203 *> (final value of Z) = U
204 *> where U is the unitary matrix in (*) (regard-
205 *> less of the value of JOB.)
206 *>
207 *> If INFO > 0 and COMPZ = 'N', then Z is not
208 *> accessed.
209 *> \endverbatim
210 *
211 * Authors:
212 * ========
213 *
214 *> \author Univ. of Tennessee
215 *> \author Univ. of California Berkeley
216 *> \author Univ. of Colorado Denver
217 *> \author NAG Ltd.
218 *
219 *> \ingroup complex16OTHERcomputational
220 *
221 *> \par Contributors:
222 * ==================
223 *>
224 *> Karen Braman and Ralph Byers, Department of Mathematics,
225 *> University of Kansas, USA
226 *
227 *> \par Further Details:
228 * =====================
229 *>
230 *> \verbatim
231 *>
232 *> Default values supplied by
233 *> ILAENV(ISPEC,'ZHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
234 *> It is suggested that these defaults be adjusted in order
235 *> to attain best performance in each particular
236 *> computational environment.
237 *>
238 *> ISPEC=12: The ZLAHQR vs ZLAQR0 crossover point.
239 *> Default: 75. (Must be at least 11.)
240 *>
241 *> ISPEC=13: Recommended deflation window size.
242 *> This depends on ILO, IHI and NS. NS is the
243 *> number of simultaneous shifts returned
244 *> by ILAENV(ISPEC=15). (See ISPEC=15 below.)
245 *> The default for (IHI-ILO+1) <= 500 is NS.
246 *> The default for (IHI-ILO+1) > 500 is 3*NS/2.
247 *>
248 *> ISPEC=14: Nibble crossover point. (See IPARMQ for
249 *> details.) Default: 14% of deflation window
250 *> size.
251 *>
252 *> ISPEC=15: Number of simultaneous shifts in a multishift
253 *> QR iteration.
254 *>
255 *> If IHI-ILO+1 is ...
256 *>
257 *> greater than ...but less ... the
258 *> or equal to ... than default is
259 *>
260 *> 1 30 NS = 2(+)
261 *> 30 60 NS = 4(+)
262 *> 60 150 NS = 10(+)
263 *> 150 590 NS = **
264 *> 590 3000 NS = 64
265 *> 3000 6000 NS = 128
266 *> 6000 infinity NS = 256
267 *>
268 *> (+) By default some or all matrices of this order
269 *> are passed to the implicit double shift routine
270 *> ZLAHQR and this parameter is ignored. See
271 *> ISPEC=12 above and comments in IPARMQ for
272 *> details.
273 *>
274 *> (**) The asterisks (**) indicate an ad-hoc
275 *> function of N increasing from 10 to 64.
276 *>
277 *> ISPEC=16: Select structured matrix multiply.
278 *> If the number of simultaneous shifts (specified
279 *> by ISPEC=15) is less than 14, then the default
280 *> for ISPEC=16 is 0. Otherwise the default for
281 *> ISPEC=16 is 2.
282 *> \endverbatim
283 *
284 *> \par References:
285 * ================
286 *>
287 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
288 *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
289 *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
290 *> 929--947, 2002.
291 *> \n
292 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
293 *> Algorithm Part II: Aggressive Early Deflation, SIAM Journal
294 *> of Matrix Analysis, volume 23, pages 948--973, 2002.
295 *
296 * =====================================================================
297  SUBROUTINE zhseqr( JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ,
298  $ WORK, LWORK, INFO )
299 *
300 * -- LAPACK computational routine --
301 * -- LAPACK is a software package provided by Univ. of Tennessee, --
302 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
303 *
304 * .. Scalar Arguments ..
305  INTEGER IHI, ILO, INFO, LDH, LDZ, LWORK, N
306  CHARACTER COMPZ, JOB
307 * ..
308 * .. Array Arguments ..
309  COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
310 * ..
311 *
312 * =====================================================================
313 *
314 * .. Parameters ..
315 *
316 * ==== Matrices of order NTINY or smaller must be processed by
317 * . ZLAHQR because of insufficient subdiagonal scratch space.
318 * . (This is a hard limit.) ====
319  INTEGER NTINY
320  parameter( ntiny = 15 )
321 *
322 * ==== NL allocates some local workspace to help small matrices
323 * . through a rare ZLAHQR failure. NL > NTINY = 15 is
324 * . required and NL <= NMIN = ILAENV(ISPEC=12,...) is recom-
325 * . mended. (The default value of NMIN is 75.) Using NL = 49
326 * . allows up to six simultaneous shifts and a 16-by-16
327 * . deflation window. ====
328  INTEGER NL
329  parameter( nl = 49 )
330  COMPLEX*16 ZERO, ONE
331  parameter( zero = ( 0.0d0, 0.0d0 ),
332  $ one = ( 1.0d0, 0.0d0 ) )
333  DOUBLE PRECISION RZERO
334  parameter( rzero = 0.0d0 )
335 * ..
336 * .. Local Arrays ..
337  COMPLEX*16 HL( NL, NL ), WORKL( NL )
338 * ..
339 * .. Local Scalars ..
340  INTEGER KBOT, NMIN
341  LOGICAL INITZ, LQUERY, WANTT, WANTZ
342 * ..
343 * .. External Functions ..
344  INTEGER ILAENV
345  LOGICAL LSAME
346  EXTERNAL ilaenv, lsame
347 * ..
348 * .. External Subroutines ..
349  EXTERNAL xerbla, zcopy, zlacpy, zlahqr, zlaqr0, zlaset
350 * ..
351 * .. Intrinsic Functions ..
352  INTRINSIC dble, dcmplx, max, min
353 * ..
354 * .. Executable Statements ..
355 *
356 * ==== Decode and check the input parameters. ====
357 *
358  wantt = lsame( job, 'S' )
359  initz = lsame( compz, 'I' )
360  wantz = initz .OR. lsame( compz, 'V' )
361  work( 1 ) = dcmplx( dble( max( 1, n ) ), rzero )
362  lquery = lwork.EQ.-1
363 *
364  info = 0
365  IF( .NOT.lsame( job, 'E' ) .AND. .NOT.wantt ) THEN
366  info = -1
367  ELSE IF( .NOT.lsame( compz, 'N' ) .AND. .NOT.wantz ) THEN
368  info = -2
369  ELSE IF( n.LT.0 ) THEN
370  info = -3
371  ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
372  info = -4
373  ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
374  info = -5
375  ELSE IF( ldh.LT.max( 1, n ) ) THEN
376  info = -7
377  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.max( 1, n ) ) ) THEN
378  info = -10
379  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
380  info = -12
381  END IF
382 *
383  IF( info.NE.0 ) THEN
384 *
385 * ==== Quick return in case of invalid argument. ====
386 *
387  CALL xerbla( 'ZHSEQR', -info )
388  RETURN
389 *
390  ELSE IF( n.EQ.0 ) THEN
391 *
392 * ==== Quick return in case N = 0; nothing to do. ====
393 *
394  RETURN
395 *
396  ELSE IF( lquery ) THEN
397 *
398 * ==== Quick return in case of a workspace query ====
399 *
400  CALL zlaqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi, z,
401  $ ldz, work, lwork, info )
402 * ==== Ensure reported workspace size is backward-compatible with
403 * . previous LAPACK versions. ====
404  work( 1 ) = dcmplx( max( dble( work( 1 ) ), dble( max( 1,
405  $ n ) ) ), rzero )
406  RETURN
407 *
408  ELSE
409 *
410 * ==== copy eigenvalues isolated by ZGEBAL ====
411 *
412  IF( ilo.GT.1 )
413  $ CALL zcopy( ilo-1, h, ldh+1, w, 1 )
414  IF( ihi.LT.n )
415  $ CALL zcopy( n-ihi, h( ihi+1, ihi+1 ), ldh+1, w( ihi+1 ), 1 )
416 *
417 * ==== Initialize Z, if requested ====
418 *
419  IF( initz )
420  $ CALL zlaset( 'A', n, n, zero, one, z, ldz )
421 *
422 * ==== Quick return if possible ====
423 *
424  IF( ilo.EQ.ihi ) THEN
425  w( ilo ) = h( ilo, ilo )
426  RETURN
427  END IF
428 *
429 * ==== ZLAHQR/ZLAQR0 crossover point ====
430 *
431  nmin = ilaenv( 12, 'ZHSEQR', job( : 1 ) // compz( : 1 ), n,
432  $ ilo, ihi, lwork )
433  nmin = max( ntiny, nmin )
434 *
435 * ==== ZLAQR0 for big matrices; ZLAHQR for small ones ====
436 *
437  IF( n.GT.nmin ) THEN
438  CALL zlaqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,
439  $ z, ldz, work, lwork, info )
440  ELSE
441 *
442 * ==== Small matrix ====
443 *
444  CALL zlahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,
445  $ z, ldz, info )
446 *
447  IF( info.GT.0 ) THEN
448 *
449 * ==== A rare ZLAHQR failure! ZLAQR0 sometimes succeeds
450 * . when ZLAHQR fails. ====
451 *
452  kbot = info
453 *
454  IF( n.GE.nl ) THEN
455 *
456 * ==== Larger matrices have enough subdiagonal scratch
457 * . space to call ZLAQR0 directly. ====
458 *
459  CALL zlaqr0( wantt, wantz, n, ilo, kbot, h, ldh, w,
460  $ ilo, ihi, z, ldz, work, lwork, info )
461 *
462  ELSE
463 *
464 * ==== Tiny matrices don't have enough subdiagonal
465 * . scratch space to benefit from ZLAQR0. Hence,
466 * . tiny matrices must be copied into a larger
467 * . array before calling ZLAQR0. ====
468 *
469  CALL zlacpy( 'A', n, n, h, ldh, hl, nl )
470  hl( n+1, n ) = zero
471  CALL zlaset( 'A', nl, nl-n, zero, zero, hl( 1, n+1 ),
472  $ nl )
473  CALL zlaqr0( wantt, wantz, nl, ilo, kbot, hl, nl, w,
474  $ ilo, ihi, z, ldz, workl, nl, info )
475  IF( wantt .OR. info.NE.0 )
476  $ CALL zlacpy( 'A', n, n, hl, nl, h, ldh )
477  END IF
478  END IF
479  END IF
480 *
481 * ==== Clear out the trash, if necessary. ====
482 *
483  IF( ( wantt .OR. info.NE.0 ) .AND. n.GT.2 )
484  $ CALL zlaset( 'L', n-2, n-2, zero, zero, h( 3, 1 ), ldh )
485 *
486 * ==== Ensure reported workspace size is backward-compatible with
487 * . previous LAPACK versions. ====
488 *
489  work( 1 ) = dcmplx( max( dble( max( 1, n ) ),
490  $ dble( work( 1 ) ) ), rzero )
491  END IF
492 *
493 * ==== End of ZHSEQR ====
494 *
495  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition: zlahqr.f:195
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:106
subroutine zlaqr0(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
ZLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition: zlaqr0.f:241
subroutine zhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
ZHSEQR
Definition: zhseqr.f:299