LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
zgees.f
Go to the documentation of this file.
1 *> \brief <b> ZGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZGEES + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgees.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgees.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgees.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
22 * LDVS, WORK, LWORK, RWORK, BWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBVS, SORT
26 * INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
27 * ..
28 * .. Array Arguments ..
29 * LOGICAL BWORK( * )
30 * DOUBLE PRECISION RWORK( * )
31 * COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
32 * ..
33 * .. Function Arguments ..
34 * LOGICAL SELECT
35 * EXTERNAL SELECT
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> ZGEES computes for an N-by-N complex nonsymmetric matrix A, the
45 *> eigenvalues, the Schur form T, and, optionally, the matrix of Schur
46 *> vectors Z. This gives the Schur factorization A = Z*T*(Z**H).
47 *>
48 *> Optionally, it also orders the eigenvalues on the diagonal of the
49 *> Schur form so that selected eigenvalues are at the top left.
50 *> The leading columns of Z then form an orthonormal basis for the
51 *> invariant subspace corresponding to the selected eigenvalues.
52 *>
53 *> A complex matrix is in Schur form if it is upper triangular.
54 *> \endverbatim
55 *
56 * Arguments:
57 * ==========
58 *
59 *> \param[in] JOBVS
60 *> \verbatim
61 *> JOBVS is CHARACTER*1
62 *> = 'N': Schur vectors are not computed;
63 *> = 'V': Schur vectors are computed.
64 *> \endverbatim
65 *>
66 *> \param[in] SORT
67 *> \verbatim
68 *> SORT is CHARACTER*1
69 *> Specifies whether or not to order the eigenvalues on the
70 *> diagonal of the Schur form.
71 *> = 'N': Eigenvalues are not ordered:
72 *> = 'S': Eigenvalues are ordered (see SELECT).
73 *> \endverbatim
74 *>
75 *> \param[in] SELECT
76 *> \verbatim
77 *> SELECT is a LOGICAL FUNCTION of one COMPLEX*16 argument
78 *> SELECT must be declared EXTERNAL in the calling subroutine.
79 *> If SORT = 'S', SELECT is used to select eigenvalues to order
80 *> to the top left of the Schur form.
81 *> IF SORT = 'N', SELECT is not referenced.
82 *> The eigenvalue W(j) is selected if SELECT(W(j)) is true.
83 *> \endverbatim
84 *>
85 *> \param[in] N
86 *> \verbatim
87 *> N is INTEGER
88 *> The order of the matrix A. N >= 0.
89 *> \endverbatim
90 *>
91 *> \param[in,out] A
92 *> \verbatim
93 *> A is COMPLEX*16 array, dimension (LDA,N)
94 *> On entry, the N-by-N matrix A.
95 *> On exit, A has been overwritten by its Schur form T.
96 *> \endverbatim
97 *>
98 *> \param[in] LDA
99 *> \verbatim
100 *> LDA is INTEGER
101 *> The leading dimension of the array A. LDA >= max(1,N).
102 *> \endverbatim
103 *>
104 *> \param[out] SDIM
105 *> \verbatim
106 *> SDIM is INTEGER
107 *> If SORT = 'N', SDIM = 0.
108 *> If SORT = 'S', SDIM = number of eigenvalues for which
109 *> SELECT is true.
110 *> \endverbatim
111 *>
112 *> \param[out] W
113 *> \verbatim
114 *> W is COMPLEX*16 array, dimension (N)
115 *> W contains the computed eigenvalues, in the same order that
116 *> they appear on the diagonal of the output Schur form T.
117 *> \endverbatim
118 *>
119 *> \param[out] VS
120 *> \verbatim
121 *> VS is COMPLEX*16 array, dimension (LDVS,N)
122 *> If JOBVS = 'V', VS contains the unitary matrix Z of Schur
123 *> vectors.
124 *> If JOBVS = 'N', VS is not referenced.
125 *> \endverbatim
126 *>
127 *> \param[in] LDVS
128 *> \verbatim
129 *> LDVS is INTEGER
130 *> The leading dimension of the array VS. LDVS >= 1; if
131 *> JOBVS = 'V', LDVS >= N.
132 *> \endverbatim
133 *>
134 *> \param[out] WORK
135 *> \verbatim
136 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
137 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
138 *> \endverbatim
139 *>
140 *> \param[in] LWORK
141 *> \verbatim
142 *> LWORK is INTEGER
143 *> The dimension of the array WORK. LWORK >= max(1,2*N).
144 *> For good performance, LWORK must generally be larger.
145 *>
146 *> If LWORK = -1, then a workspace query is assumed; the routine
147 *> only calculates the optimal size of the WORK array, returns
148 *> this value as the first entry of the WORK array, and no error
149 *> message related to LWORK is issued by XERBLA.
150 *> \endverbatim
151 *>
152 *> \param[out] RWORK
153 *> \verbatim
154 *> RWORK is DOUBLE PRECISION array, dimension (N)
155 *> \endverbatim
156 *>
157 *> \param[out] BWORK
158 *> \verbatim
159 *> BWORK is LOGICAL array, dimension (N)
160 *> Not referenced if SORT = 'N'.
161 *> \endverbatim
162 *>
163 *> \param[out] INFO
164 *> \verbatim
165 *> INFO is INTEGER
166 *> = 0: successful exit
167 *> < 0: if INFO = -i, the i-th argument had an illegal value.
168 *> > 0: if INFO = i, and i is
169 *> <= N: the QR algorithm failed to compute all the
170 *> eigenvalues; elements 1:ILO-1 and i+1:N of W
171 *> contain those eigenvalues which have converged;
172 *> if JOBVS = 'V', VS contains the matrix which
173 *> reduces A to its partially converged Schur form.
174 *> = N+1: the eigenvalues could not be reordered because
175 *> some eigenvalues were too close to separate (the
176 *> problem is very ill-conditioned);
177 *> = N+2: after reordering, roundoff changed values of
178 *> some complex eigenvalues so that leading
179 *> eigenvalues in the Schur form no longer satisfy
180 *> SELECT = .TRUE.. This could also be caused by
181 *> underflow due to scaling.
182 *> \endverbatim
183 *
184 * Authors:
185 * ========
186 *
187 *> \author Univ. of Tennessee
188 *> \author Univ. of California Berkeley
189 *> \author Univ. of Colorado Denver
190 *> \author NAG Ltd.
191 *
192 *> \date November 2011
193 *
194 *> \ingroup complex16GEeigen
195 *
196 * =====================================================================
197  SUBROUTINE zgees( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
198  $ ldvs, work, lwork, rwork, bwork, info )
199 *
200 * -- LAPACK driver routine (version 3.4.0) --
201 * -- LAPACK is a software package provided by Univ. of Tennessee, --
202 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
203 * November 2011
204 *
205 * .. Scalar Arguments ..
206  CHARACTER JOBVS, SORT
207  INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
208 * ..
209 * .. Array Arguments ..
210  LOGICAL BWORK( * )
211  DOUBLE PRECISION RWORK( * )
212  COMPLEX*16 A( lda, * ), VS( ldvs, * ), W( * ), WORK( * )
213 * ..
214 * .. Function Arguments ..
215  LOGICAL SELECT
216  EXTERNAL SELECT
217 * ..
218 *
219 * =====================================================================
220 *
221 * .. Parameters ..
222  DOUBLE PRECISION ZERO, ONE
223  parameter ( zero = 0.0d0, one = 1.0d0 )
224 * ..
225 * .. Local Scalars ..
226  LOGICAL LQUERY, SCALEA, WANTST, WANTVS
227  INTEGER HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
228  $ itau, iwrk, maxwrk, minwrk
229  DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
230 * ..
231 * .. Local Arrays ..
232  DOUBLE PRECISION DUM( 1 )
233 * ..
234 * .. External Subroutines ..
235  EXTERNAL dlabad, xerbla, zcopy, zgebak, zgebal, zgehrd,
237 * ..
238 * .. External Functions ..
239  LOGICAL LSAME
240  INTEGER ILAENV
241  DOUBLE PRECISION DLAMCH, ZLANGE
242  EXTERNAL lsame, ilaenv, dlamch, zlange
243 * ..
244 * .. Intrinsic Functions ..
245  INTRINSIC max, sqrt
246 * ..
247 * .. Executable Statements ..
248 *
249 * Test the input arguments
250 *
251  info = 0
252  lquery = ( lwork.EQ.-1 )
253  wantvs = lsame( jobvs, 'V' )
254  wantst = lsame( sort, 'S' )
255  IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
256  info = -1
257  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
258  info = -2
259  ELSE IF( n.LT.0 ) THEN
260  info = -4
261  ELSE IF( lda.LT.max( 1, n ) ) THEN
262  info = -6
263  ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
264  info = -10
265  END IF
266 *
267 * Compute workspace
268 * (Note: Comments in the code beginning "Workspace:" describe the
269 * minimal amount of workspace needed at that point in the code,
270 * as well as the preferred amount for good performance.
271 * CWorkspace refers to complex workspace, and RWorkspace to real
272 * workspace. NB refers to the optimal block size for the
273 * immediately following subroutine, as returned by ILAENV.
274 * HSWORK refers to the workspace preferred by ZHSEQR, as
275 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
276 * the worst case.)
277 *
278  IF( info.EQ.0 ) THEN
279  IF( n.EQ.0 ) THEN
280  minwrk = 1
281  maxwrk = 1
282  ELSE
283  maxwrk = n + n*ilaenv( 1, 'ZGEHRD', ' ', n, 1, n, 0 )
284  minwrk = 2*n
285 *
286  CALL zhseqr( 'S', jobvs, n, 1, n, a, lda, w, vs, ldvs,
287  $ work, -1, ieval )
288  hswork = work( 1 )
289 *
290  IF( .NOT.wantvs ) THEN
291  maxwrk = max( maxwrk, hswork )
292  ELSE
293  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'ZUNGHR',
294  $ ' ', n, 1, n, -1 ) )
295  maxwrk = max( maxwrk, hswork )
296  END IF
297  END IF
298  work( 1 ) = maxwrk
299 *
300  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
301  info = -12
302  END IF
303  END IF
304 *
305  IF( info.NE.0 ) THEN
306  CALL xerbla( 'ZGEES ', -info )
307  RETURN
308  ELSE IF( lquery ) THEN
309  RETURN
310  END IF
311 *
312 * Quick return if possible
313 *
314  IF( n.EQ.0 ) THEN
315  sdim = 0
316  RETURN
317  END IF
318 *
319 * Get machine constants
320 *
321  eps = dlamch( 'P' )
322  smlnum = dlamch( 'S' )
323  bignum = one / smlnum
324  CALL dlabad( smlnum, bignum )
325  smlnum = sqrt( smlnum ) / eps
326  bignum = one / smlnum
327 *
328 * Scale A if max element outside range [SMLNUM,BIGNUM]
329 *
330  anrm = zlange( 'M', n, n, a, lda, dum )
331  scalea = .false.
332  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
333  scalea = .true.
334  cscale = smlnum
335  ELSE IF( anrm.GT.bignum ) THEN
336  scalea = .true.
337  cscale = bignum
338  END IF
339  IF( scalea )
340  $ CALL zlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
341 *
342 * Permute the matrix to make it more nearly triangular
343 * (CWorkspace: none)
344 * (RWorkspace: need N)
345 *
346  ibal = 1
347  CALL zgebal( 'P', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
348 *
349 * Reduce to upper Hessenberg form
350 * (CWorkspace: need 2*N, prefer N+N*NB)
351 * (RWorkspace: none)
352 *
353  itau = 1
354  iwrk = n + itau
355  CALL zgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
356  $ lwork-iwrk+1, ierr )
357 *
358  IF( wantvs ) THEN
359 *
360 * Copy Householder vectors to VS
361 *
362  CALL zlacpy( 'L', n, n, a, lda, vs, ldvs )
363 *
364 * Generate unitary matrix in VS
365 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
366 * (RWorkspace: none)
367 *
368  CALL zunghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
369  $ lwork-iwrk+1, ierr )
370  END IF
371 *
372  sdim = 0
373 *
374 * Perform QR iteration, accumulating Schur vectors in VS if desired
375 * (CWorkspace: need 1, prefer HSWORK (see comments) )
376 * (RWorkspace: none)
377 *
378  iwrk = itau
379  CALL zhseqr( 'S', jobvs, n, ilo, ihi, a, lda, w, vs, ldvs,
380  $ work( iwrk ), lwork-iwrk+1, ieval )
381  IF( ieval.GT.0 )
382  $ info = ieval
383 *
384 * Sort eigenvalues if desired
385 *
386  IF( wantst .AND. info.EQ.0 ) THEN
387  IF( scalea )
388  $ CALL zlascl( 'G', 0, 0, cscale, anrm, n, 1, w, n, ierr )
389  DO 10 i = 1, n
390  bwork( i ) = SELECT( w( i ) )
391  10 CONTINUE
392 *
393 * Reorder eigenvalues and transform Schur vectors
394 * (CWorkspace: none)
395 * (RWorkspace: none)
396 *
397  CALL ztrsen( 'N', jobvs, bwork, n, a, lda, vs, ldvs, w, sdim,
398  $ s, sep, work( iwrk ), lwork-iwrk+1, icond )
399  END IF
400 *
401  IF( wantvs ) THEN
402 *
403 * Undo balancing
404 * (CWorkspace: none)
405 * (RWorkspace: need N)
406 *
407  CALL zgebak( 'P', 'R', n, ilo, ihi, rwork( ibal ), n, vs, ldvs,
408  $ ierr )
409  END IF
410 *
411  IF( scalea ) THEN
412 *
413 * Undo scaling for the Schur form of A
414 *
415  CALL zlascl( 'U', 0, 0, cscale, anrm, n, n, a, lda, ierr )
416  CALL zcopy( n, a, lda+1, w, 1 )
417  END IF
418 *
419  work( 1 ) = maxwrk
420  RETURN
421 *
422 * End of ZGEES
423 *
424  END
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
ZGEBAK
Definition: zgebak.f:133
subroutine zgees(JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS, LDVS, WORK, LWORK, RWORK, BWORK, INFO)
ZGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition: zgees.f:199
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
Definition: zgehrd.f:169
subroutine ztrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, SEP, WORK, LWORK, INFO)
ZTRSEN
Definition: ztrsen.f:266
subroutine zunghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGHR
Definition: zunghr.f:128
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine zgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
ZGEBAL
Definition: zgebal.f:162
subroutine zhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
ZHSEQR
Definition: zhseqr.f:301
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:145