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