LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
chseqr.f
Go to the documentation of this file.
1*> \brief \b CHSEQR
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CHSEQR + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chseqr.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chseqr.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chseqr.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CHSEQR( 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 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> CHSEQR 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 CGEBAL, and then passed to ZGEHRD
87*> when the matrix output by CGEBAL 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 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 CHSEQR, 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 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 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 CUNGHR
133*> after the call to CGEHRD 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 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 CHSEQR does a workspace query.
163*> In this case, CHSEQR 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, CHSEQR 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,'CHSEQR',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 CLAHQR vs CLAQR0 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*> CLAHQR 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 chseqr( 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 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* . CLAHQR 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 CLAHQR 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 ZERO, ONE
331 parameter( zero = ( 0.0e0, 0.0e0 ),
332 $ one = ( 1.0e0, 0.0e0 ) )
333 REAL RZERO
334 parameter( rzero = 0.0e0 )
335* ..
336* .. Local Arrays ..
337 COMPLEX 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 REAL SROUNDUP_LWORK
347 EXTERNAL ilaenv, lsame, sroundup_lwork
348* ..
349* .. External Subroutines ..
350 EXTERNAL ccopy, clacpy, clahqr, claqr0, claset, xerbla
351* ..
352* .. Intrinsic Functions ..
353 INTRINSIC cmplx, max, min, real
354* ..
355* .. Executable Statements ..
356*
357* ==== Decode and check the input parameters. ====
358*
359 wantt = lsame( job, 'S' )
360 initz = lsame( compz, 'I' )
361 wantz = initz .OR. lsame( compz, 'V' )
362 work( 1 ) = cmplx( real( max( 1, n ) ), rzero )
363 lquery = lwork.EQ.-1
364*
365 info = 0
366 IF( .NOT.lsame( job, 'E' ) .AND. .NOT.wantt ) THEN
367 info = -1
368 ELSE IF( .NOT.lsame( compz, 'N' ) .AND. .NOT.wantz ) THEN
369 info = -2
370 ELSE IF( n.LT.0 ) THEN
371 info = -3
372 ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
373 info = -4
374 ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
375 info = -5
376 ELSE IF( ldh.LT.max( 1, n ) ) THEN
377 info = -7
378 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.max( 1, n ) ) ) THEN
379 info = -10
380 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
381 info = -12
382 END IF
383*
384 IF( info.NE.0 ) THEN
385*
386* ==== Quick return in case of invalid argument. ====
387*
388 CALL xerbla( 'CHSEQR', -info )
389 RETURN
390*
391 ELSE IF( n.EQ.0 ) THEN
392*
393* ==== Quick return in case N = 0; nothing to do. ====
394*
395 RETURN
396*
397 ELSE IF( lquery ) THEN
398*
399* ==== Quick return in case of a workspace query ====
400*
401 CALL claqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi, z,
402 $ ldz, work, lwork, info )
403* ==== Ensure reported workspace size is backward-compatible with
404* . previous LAPACK versions. ====
405 work( 1 ) = cmplx( max( real( work( 1 ) ), real( max( 1,
406 $ n ) ) ), rzero )
407 RETURN
408*
409 ELSE
410*
411* ==== copy eigenvalues isolated by CGEBAL ====
412*
413 IF( ilo.GT.1 )
414 $ CALL ccopy( ilo-1, h, ldh+1, w, 1 )
415 IF( ihi.LT.n )
416 $ CALL ccopy( n-ihi, h( ihi+1, ihi+1 ), ldh+1, w( ihi+1 ), 1 )
417*
418* ==== Initialize Z, if requested ====
419*
420 IF( initz )
421 $ CALL claset( 'A', n, n, zero, one, z, ldz )
422*
423* ==== Quick return if possible ====
424*
425 IF( ilo.EQ.ihi ) THEN
426 w( ilo ) = h( ilo, ilo )
427 RETURN
428 END IF
429*
430* ==== CLAHQR/CLAQR0 crossover point ====
431*
432 nmin = ilaenv( 12, 'CHSEQR', job( : 1 ) // compz( : 1 ), n,
433 $ ilo, ihi, lwork )
434 nmin = max( ntiny, nmin )
435*
436* ==== CLAQR0 for big matrices; CLAHQR for small ones ====
437*
438 IF( n.GT.nmin ) THEN
439 CALL claqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,
440 $ z, ldz, work, lwork, info )
441 ELSE
442*
443* ==== Small matrix ====
444*
445 CALL clahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,
446 $ z, ldz, info )
447*
448 IF( info.GT.0 ) THEN
449*
450* ==== A rare CLAHQR failure! CLAQR0 sometimes succeeds
451* . when CLAHQR fails. ====
452*
453 kbot = info
454*
455 IF( n.GE.nl ) THEN
456*
457* ==== Larger matrices have enough subdiagonal scratch
458* . space to call CLAQR0 directly. ====
459*
460 CALL claqr0( wantt, wantz, n, ilo, kbot, h, ldh, w,
461 $ ilo, ihi, z, ldz, work, lwork, info )
462*
463 ELSE
464*
465* ==== Tiny matrices don't have enough subdiagonal
466* . scratch space to benefit from CLAQR0. Hence,
467* . tiny matrices must be copied into a larger
468* . array before calling CLAQR0. ====
469*
470 CALL clacpy( 'A', n, n, h, ldh, hl, nl )
471 hl( n+1, n ) = zero
472 CALL claset( 'A', nl, nl-n, zero, zero, hl( 1, n+1 ),
473 $ nl )
474 CALL claqr0( wantt, wantz, nl, ilo, kbot, hl, nl, w,
475 $ ilo, ihi, z, ldz, workl, nl, info )
476 IF( wantt .OR. info.NE.0 )
477 $ CALL clacpy( 'A', n, n, hl, nl, h, ldh )
478 END IF
479 END IF
480 END IF
481*
482* ==== Clear out the trash, if necessary. ====
483*
484 IF( ( wantt .OR. info.NE.0 ) .AND. n.GT.2 )
485 $ CALL claset( 'L', n-2, n-2, zero, zero, h( 3, 1 ), ldh )
486*
487* ==== Ensure reported workspace size is backward-compatible with
488* . previous LAPACK versions. ====
489*
490 work( 1 ) = cmplx( max( real( max( 1, n ) ),
491 $ real( work( 1 ) ) ), rzero )
492 END IF
493*
494* ==== End of CHSEQR ====
495*
496 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine chseqr(job, compz, n, ilo, ihi, h, ldh, w, z, ldz, work, lwork, info)
CHSEQR
Definition chseqr.f:299
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
subroutine clahqr(wantt, wantz, n, ilo, ihi, h, ldh, w, iloz, ihiz, z, ldz, info)
CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition clahqr.f:195
subroutine claqr0(wantt, wantz, n, ilo, ihi, h, ldh, w, iloz, ihiz, z, ldz, work, lwork, info)
CLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition claqr0.f:240
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106