LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
dhseqr.f
Go to the documentation of this file.
1 *> \brief \b DHSEQR
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DHSEQR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dhseqr.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dhseqr.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dhseqr.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DHSEQR( 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 * DOUBLE PRECISION H( LDH, * ), WI( * ), WORK( * ), WR( * ),
30 * $ Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DHSEQR 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 DGEBAL, and then passed to ZGEHRD
88 *> when the matrix output by DGEBAL 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 DOUBLE PRECISION 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 DHSEQR, 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 DOUBLE PRECISION array, dimension (N)
122 *> \endverbatim
123 *>
124 *> \param[out] WI
125 *> \verbatim
126 *> WI is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DORGHR
150 *> after the call to DGEHRD 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 DOUBLE PRECISION 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 DHSEQR does a workspace query.
180 *> In this case, DHSEQR 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, DHSEQR 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 doubleOTHERcomputational
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,'DHSEQR',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 DLAHQR vs DLAQR0 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 *> DLAHQR 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 dhseqr( 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  DOUBLE PRECISION 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 * . DLAHQR 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 DLAHQR 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  DOUBLE PRECISION ZERO, ONE
352  parameter ( zero = 0.0d0, one = 1.0d0 )
353 * ..
354 * .. Local Arrays ..
355  DOUBLE PRECISION 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 dlacpy, dlahqr, dlaqr0, dlaset, xerbla
368 * ..
369 * .. Intrinsic Functions ..
370  INTRINSIC dble, max, min
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 ) = dble( 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( 'DHSEQR', -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 dlaqr0( 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( dble( max( 1, n ) ), work( 1 ) )
423  RETURN
424 *
425  ELSE
426 *
427 * ==== copy eigenvalues isolated by DGEBAL ====
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 dlaset( '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 * ==== DLAHQR/DLAQR0 crossover point ====
452 *
453  nmin = ilaenv( 12, 'DHSEQR', job( : 1 ) // compz( : 1 ), n,
454  $ ilo, ihi, lwork )
455  nmin = max( ntiny, nmin )
456 *
457 * ==== DLAQR0 for big matrices; DLAHQR for small ones ====
458 *
459  IF( n.GT.nmin ) THEN
460  CALL dlaqr0( 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 dlahqr( 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 DLAHQR failure! DLAQR0 sometimes succeeds
472 * . when DLAHQR fails. ====
473 *
474  kbot = info
475 *
476  IF( n.GE.nl ) THEN
477 *
478 * ==== Larger matrices have enough subdiagonal scratch
479 * . space to call DLAQR0 directly. ====
480 *
481  CALL dlaqr0( 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 DLAQR0. Hence,
488 * . tiny matrices must be copied into a larger
489 * . array before calling DLAQR0. ====
490 *
491  CALL dlacpy( 'A', n, n, h, ldh, hl, nl )
492  hl( n+1, n ) = zero
493  CALL dlaset( 'A', nl, nl-n, zero, zero, hl( 1, n+1 ),
494  $ nl )
495  CALL dlaqr0( 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 dlacpy( '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 dlaset( '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( dble( max( 1, n ) ), work( 1 ) )
512  END IF
513 *
514 * ==== End of DHSEQR ====
515 *
516  END
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)
DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
Definition: dlahqr.f:209
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlaqr0(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
DLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition: dlaqr0.f:258
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
Definition: dhseqr.f:318