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