LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 hseqr
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)
Definition cblat2.f:3285
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
subroutine zhseqr(job, compz, n, ilo, ihi, h, ldh, w, z, ldz, work, lwork, info)
ZHSEQR
Definition zhseqr.f:299
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 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 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 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