LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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