LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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)*H*(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 .GE. 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.GT.0, then 1.LE.ILO.LE.IHI.LE.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.GT.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.GT.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 .GE. 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.GT.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.GE.MAX(1,N). Otherwize, LDZ.GE.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 .GE. 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 *> .LT. 0: if INFO = -i, the i-th argument had an illegal
175 *> value
176 *> .GT. 0: if INFO = i, CHSEQR failed to compute all of
177 *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
178 *> and WI contain those eigenvalues which have been
179 *> successfully computed. (Failures are rare.)
180 *>
181 *> If INFO .GT. 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 .GT. 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 .GT. 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 .GT. 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 .GT. 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 *> \date November 2011
220 *
221 *> \ingroup complexOTHERcomputational
222 *
223 *> \par Contributors:
224 * ==================
225 *>
226 *> Karen Braman and Ralph Byers, Department of Mathematics,
227 *> University of Kansas, USA
228 *
229 *> \par Further Details:
230 * =====================
231 *>
232 *> \verbatim
233 *>
234 *> Default values supplied by
235 *> ILAENV(ISPEC,'CHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
236 *> It is suggested that these defaults be adjusted in order
237 *> to attain best performance in each particular
238 *> computational environment.
239 *>
240 *> ISPEC=12: The CLAHQR vs CLAQR0 crossover point.
241 *> Default: 75. (Must be at least 11.)
242 *>
243 *> ISPEC=13: Recommended deflation window size.
244 *> This depends on ILO, IHI and NS. NS is the
245 *> number of simultaneous shifts returned
246 *> by ILAENV(ISPEC=15). (See ISPEC=15 below.)
247 *> The default for (IHI-ILO+1).LE.500 is NS.
248 *> The default for (IHI-ILO+1).GT.500 is 3*NS/2.
249 *>
250 *> ISPEC=14: Nibble crossover point. (See IPARMQ for
251 *> details.) Default: 14% of deflation window
252 *> size.
253 *>
254 *> ISPEC=15: Number of simultaneous shifts in a multishift
255 *> QR iteration.
256 *>
257 *> If IHI-ILO+1 is ...
258 *>
259 *> greater than ...but less ... the
260 *> or equal to ... than default is
261 *>
262 *> 1 30 NS = 2(+)
263 *> 30 60 NS = 4(+)
264 *> 60 150 NS = 10(+)
265 *> 150 590 NS = **
266 *> 590 3000 NS = 64
267 *> 3000 6000 NS = 128
268 *> 6000 infinity NS = 256
269 *>
270 *> (+) By default some or all matrices of this order
271 *> are passed to the implicit double shift routine
272 *> CLAHQR and this parameter is ignored. See
273 *> ISPEC=12 above and comments in IPARMQ for
274 *> details.
275 *>
276 *> (**) The asterisks (**) indicate an ad-hoc
277 *> function of N increasing from 10 to 64.
278 *>
279 *> ISPEC=16: Select structured matrix multiply.
280 *> If the number of simultaneous shifts (specified
281 *> by ISPEC=15) is less than 14, then the default
282 *> for ISPEC=16 is 0. Otherwise the default for
283 *> ISPEC=16 is 2.
284 *> \endverbatim
285 *
286 *> \par References:
287 * ================
288 *>
289 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
290 *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
291 *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
292 *> 929--947, 2002.
293 *> \n
294 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
295 *> Algorithm Part II: Aggressive Early Deflation, SIAM Journal
296 *> of Matrix Analysis, volume 23, pages 948--973, 2002.
297 *
298 * =====================================================================
299  SUBROUTINE chseqr( JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ,
300  $ work, lwork, info )
301 *
302 * -- LAPACK computational routine (version 3.4.0) --
303 * -- LAPACK is a software package provided by Univ. of Tennessee, --
304 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
305 * November 2011
306 *
307 * .. Scalar Arguments ..
308  INTEGER ihi, ilo, info, ldh, ldz, lwork, n
309  CHARACTER compz, job
310 * ..
311 * .. Array Arguments ..
312  COMPLEX h( ldh, * ), w( * ), work( * ), z( ldz, * )
313 * ..
314 *
315 * =====================================================================
316 *
317 * .. Parameters ..
318 *
319 * ==== Matrices of order NTINY or smaller must be processed by
320 * . CLAHQR because of insufficient subdiagonal scratch space.
321 * . (This is a hard limit.) ====
322  INTEGER ntiny
323  parameter( ntiny = 11 )
324 *
325 * ==== NL allocates some local workspace to help small matrices
326 * . through a rare CLAHQR failure. NL .GT. NTINY = 11 is
327 * . required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom-
328 * . mended. (The default value of NMIN is 75.) Using NL = 49
329 * . allows up to six simultaneous shifts and a 16-by-16
330 * . deflation window. ====
331  INTEGER nl
332  parameter( nl = 49 )
333  COMPLEX zero, one
334  parameter( zero = ( 0.0e0, 0.0e0 ),
335  $ one = ( 1.0e0, 0.0e0 ) )
336  REAL rzero
337  parameter( rzero = 0.0e0 )
338 * ..
339 * .. Local Arrays ..
340  COMPLEX hl( nl, nl ), workl( nl )
341 * ..
342 * .. Local Scalars ..
343  INTEGER kbot, nmin
344  LOGICAL initz, lquery, wantt, wantz
345 * ..
346 * .. External Functions ..
347  INTEGER ilaenv
348  LOGICAL lsame
349  EXTERNAL ilaenv, lsame
350 * ..
351 * .. External Subroutines ..
352  EXTERNAL ccopy, clacpy, clahqr, claqr0, claset, xerbla
353 * ..
354 * .. Intrinsic Functions ..
355  INTRINSIC cmplx, max, min, real
356 * ..
357 * .. Executable Statements ..
358 *
359 * ==== Decode and check the input parameters. ====
360 *
361  wantt = lsame( job, 'S' )
362  initz = lsame( compz, 'I' )
363  wantz = initz .OR. lsame( compz, 'V' )
364  work( 1 ) = cmplx( REAL( MAX( 1, N ) ), rzero )
365  lquery = lwork.EQ.-1
366 *
367  info = 0
368  IF( .NOT.lsame( job, 'E' ) .AND. .NOT.wantt ) THEN
369  info = -1
370  ELSE IF( .NOT.lsame( compz, 'N' ) .AND. .NOT.wantz ) THEN
371  info = -2
372  ELSE IF( n.LT.0 ) THEN
373  info = -3
374  ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
375  info = -4
376  ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
377  info = -5
378  ELSE IF( ldh.LT.max( 1, n ) ) THEN
379  info = -7
380  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.max( 1, n ) ) ) THEN
381  info = -10
382  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
383  info = -12
384  END IF
385 *
386  IF( info.NE.0 ) THEN
387 *
388 * ==== Quick return in case of invalid argument. ====
389 *
390  CALL xerbla( 'CHSEQR', -info )
391  return
392 *
393  ELSE IF( n.EQ.0 ) THEN
394 *
395 * ==== Quick return in case N = 0; nothing to do. ====
396 *
397  return
398 *
399  ELSE IF( lquery ) THEN
400 *
401 * ==== Quick return in case of a workspace query ====
402 *
403  CALL claqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi, z,
404  $ ldz, work, lwork, info )
405 * ==== Ensure reported workspace size is backward-compatible with
406 * . previous LAPACK versions. ====
407  work( 1 ) = cmplx( max( REAL( WORK( 1 ) ), REAL( MAX( 1, $ N ) ) ), rzero )
408  return
409 *
410  ELSE
411 *
412 * ==== copy eigenvalues isolated by CGEBAL ====
413 *
414  IF( ilo.GT.1 )
415  $ CALL ccopy( ilo-1, h, ldh+1, w, 1 )
416  IF( ihi.LT.n )
417  $ CALL ccopy( n-ihi, h( ihi+1, ihi+1 ), ldh+1, w( ihi+1 ), 1 )
418 *
419 * ==== Initialize Z, if requested ====
420 *
421  IF( initz )
422  $ CALL claset( 'A', n, n, zero, one, z, ldz )
423 *
424 * ==== Quick return if possible ====
425 *
426  IF( ilo.EQ.ihi ) THEN
427  w( ilo ) = h( ilo, ilo )
428  return
429  END IF
430 *
431 * ==== CLAHQR/CLAQR0 crossover point ====
432 *
433  nmin = ilaenv( 12, 'CHSEQR', job( : 1 ) // compz( : 1 ), n,
434  $ ilo, ihi, lwork )
435  nmin = max( ntiny, nmin )
436 *
437 * ==== CLAQR0 for big matrices; CLAHQR for small ones ====
438 *
439  IF( n.GT.nmin ) THEN
440  CALL claqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,
441  $ z, ldz, work, lwork, info )
442  ELSE
443 *
444 * ==== Small matrix ====
445 *
446  CALL clahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,
447  $ z, ldz, info )
448 *
449  IF( info.GT.0 ) THEN
450 *
451 * ==== A rare CLAHQR failure! CLAQR0 sometimes succeeds
452 * . when CLAHQR fails. ====
453 *
454  kbot = info
455 *
456  IF( n.GE.nl ) THEN
457 *
458 * ==== Larger matrices have enough subdiagonal scratch
459 * . space to call CLAQR0 directly. ====
460 *
461  CALL claqr0( wantt, wantz, n, ilo, kbot, h, ldh, w,
462  $ ilo, ihi, z, ldz, work, lwork, info )
463 *
464  ELSE
465 *
466 * ==== Tiny matrices don't have enough subdiagonal
467 * . scratch space to benefit from CLAQR0. Hence,
468 * . tiny matrices must be copied into a larger
469 * . array before calling CLAQR0. ====
470 *
471  CALL clacpy( 'A', n, n, h, ldh, hl, nl )
472  hl( n+1, n ) = zero
473  CALL claset( 'A', nl, nl-n, zero, zero, hl( 1, n+1 ),
474  $ nl )
475  CALL claqr0( wantt, wantz, nl, ilo, kbot, hl, nl, w,
476  $ ilo, ihi, z, ldz, workl, nl, info )
477  IF( wantt .OR. info.NE.0 )
478  $ CALL clacpy( 'A', n, n, hl, nl, h, ldh )
479  END IF
480  END IF
481  END IF
482 *
483 * ==== Clear out the trash, if necessary. ====
484 *
485  IF( ( wantt .OR. info.NE.0 ) .AND. n.GT.2 )
486  $ CALL claset( 'L', n-2, n-2, zero, zero, h( 3, 1 ), ldh )
487 *
488 * ==== Ensure reported workspace size is backward-compatible with
489 * . previous LAPACK versions. ====
490 *
491  work( 1 ) = cmplx( max( REAL( MAX( 1, N ) ),
492  $ REAL( WORK( 1 ) ) ), rzero )
493  END IF
494 *
495 * ==== End of CHSEQR ====
496 *
497  END
498