LAPACK  3.10.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 *> \ingroup complex16GEeigen
193 *
194 * =====================================================================
195  SUBROUTINE zgees( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
196  $ LDVS, WORK, LWORK, RWORK, BWORK, INFO )
197 *
198 * -- LAPACK driver routine --
199 * -- LAPACK is a software package provided by Univ. of Tennessee, --
200 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
201 *
202 * .. Scalar Arguments ..
203  CHARACTER JOBVS, SORT
204  INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
205 * ..
206 * .. Array Arguments ..
207  LOGICAL BWORK( * )
208  DOUBLE PRECISION RWORK( * )
209  COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
210 * ..
211 * .. Function Arguments ..
212  LOGICAL SELECT
213  EXTERNAL SELECT
214 * ..
215 *
216 * =====================================================================
217 *
218 * .. Parameters ..
219  DOUBLE PRECISION ZERO, ONE
220  parameter( zero = 0.0d0, one = 1.0d0 )
221 * ..
222 * .. Local Scalars ..
223  LOGICAL LQUERY, SCALEA, WANTST, WANTVS
224  INTEGER HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
225  $ itau, iwrk, maxwrk, minwrk
226  DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
227 * ..
228 * .. Local Arrays ..
229  DOUBLE PRECISION DUM( 1 )
230 * ..
231 * .. External Subroutines ..
232  EXTERNAL dlabad, xerbla, zcopy, zgebak, zgebal, zgehrd,
234 * ..
235 * .. External Functions ..
236  LOGICAL LSAME
237  INTEGER ILAENV
238  DOUBLE PRECISION DLAMCH, ZLANGE
239  EXTERNAL lsame, ilaenv, dlamch, zlange
240 * ..
241 * .. Intrinsic Functions ..
242  INTRINSIC max, sqrt
243 * ..
244 * .. Executable Statements ..
245 *
246 * Test the input arguments
247 *
248  info = 0
249  lquery = ( lwork.EQ.-1 )
250  wantvs = lsame( jobvs, 'V' )
251  wantst = lsame( sort, 'S' )
252  IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
253  info = -1
254  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
255  info = -2
256  ELSE IF( n.LT.0 ) THEN
257  info = -4
258  ELSE IF( lda.LT.max( 1, n ) ) THEN
259  info = -6
260  ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
261  info = -10
262  END IF
263 *
264 * Compute workspace
265 * (Note: Comments in the code beginning "Workspace:" describe the
266 * minimal amount of workspace needed at that point in the code,
267 * as well as the preferred amount for good performance.
268 * CWorkspace refers to complex workspace, and RWorkspace to real
269 * workspace. NB refers to the optimal block size for the
270 * immediately following subroutine, as returned by ILAENV.
271 * HSWORK refers to the workspace preferred by ZHSEQR, as
272 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
273 * the worst case.)
274 *
275  IF( info.EQ.0 ) THEN
276  IF( n.EQ.0 ) THEN
277  minwrk = 1
278  maxwrk = 1
279  ELSE
280  maxwrk = n + n*ilaenv( 1, 'ZGEHRD', ' ', n, 1, n, 0 )
281  minwrk = 2*n
282 *
283  CALL zhseqr( 'S', jobvs, n, 1, n, a, lda, w, vs, ldvs,
284  $ work, -1, ieval )
285  hswork = dble( work( 1 ) )
286 *
287  IF( .NOT.wantvs ) THEN
288  maxwrk = max( maxwrk, hswork )
289  ELSE
290  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'ZUNGHR',
291  $ ' ', n, 1, n, -1 ) )
292  maxwrk = max( maxwrk, hswork )
293  END IF
294  END IF
295  work( 1 ) = maxwrk
296 *
297  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
298  info = -12
299  END IF
300  END IF
301 *
302  IF( info.NE.0 ) THEN
303  CALL xerbla( 'ZGEES ', -info )
304  RETURN
305  ELSE IF( lquery ) THEN
306  RETURN
307  END IF
308 *
309 * Quick return if possible
310 *
311  IF( n.EQ.0 ) THEN
312  sdim = 0
313  RETURN
314  END IF
315 *
316 * Get machine constants
317 *
318  eps = dlamch( 'P' )
319  smlnum = dlamch( 'S' )
320  bignum = one / smlnum
321  CALL dlabad( smlnum, bignum )
322  smlnum = sqrt( smlnum ) / eps
323  bignum = one / smlnum
324 *
325 * Scale A if max element outside range [SMLNUM,BIGNUM]
326 *
327  anrm = zlange( 'M', n, n, a, lda, dum )
328  scalea = .false.
329  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
330  scalea = .true.
331  cscale = smlnum
332  ELSE IF( anrm.GT.bignum ) THEN
333  scalea = .true.
334  cscale = bignum
335  END IF
336  IF( scalea )
337  $ CALL zlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
338 *
339 * Permute the matrix to make it more nearly triangular
340 * (CWorkspace: none)
341 * (RWorkspace: need N)
342 *
343  ibal = 1
344  CALL zgebal( 'P', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
345 *
346 * Reduce to upper Hessenberg form
347 * (CWorkspace: need 2*N, prefer N+N*NB)
348 * (RWorkspace: none)
349 *
350  itau = 1
351  iwrk = n + itau
352  CALL zgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
353  $ lwork-iwrk+1, ierr )
354 *
355  IF( wantvs ) THEN
356 *
357 * Copy Householder vectors to VS
358 *
359  CALL zlacpy( 'L', n, n, a, lda, vs, ldvs )
360 *
361 * Generate unitary matrix in VS
362 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
363 * (RWorkspace: none)
364 *
365  CALL zunghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
366  $ lwork-iwrk+1, ierr )
367  END IF
368 *
369  sdim = 0
370 *
371 * Perform QR iteration, accumulating Schur vectors in VS if desired
372 * (CWorkspace: need 1, prefer HSWORK (see comments) )
373 * (RWorkspace: none)
374 *
375  iwrk = itau
376  CALL zhseqr( 'S', jobvs, n, ilo, ihi, a, lda, w, vs, ldvs,
377  $ work( iwrk ), lwork-iwrk+1, ieval )
378  IF( ieval.GT.0 )
379  $ info = ieval
380 *
381 * Sort eigenvalues if desired
382 *
383  IF( wantst .AND. info.EQ.0 ) THEN
384  IF( scalea )
385  $ CALL zlascl( 'G', 0, 0, cscale, anrm, n, 1, w, n, ierr )
386  DO 10 i = 1, n
387  bwork( i ) = SELECT( w( i ) )
388  10 CONTINUE
389 *
390 * Reorder eigenvalues and transform Schur vectors
391 * (CWorkspace: none)
392 * (RWorkspace: none)
393 *
394  CALL ztrsen( 'N', jobvs, bwork, n, a, lda, vs, ldvs, w, sdim,
395  $ s, sep, work( iwrk ), lwork-iwrk+1, icond )
396  END IF
397 *
398  IF( wantvs ) THEN
399 *
400 * Undo balancing
401 * (CWorkspace: none)
402 * (RWorkspace: need N)
403 *
404  CALL zgebak( 'P', 'R', n, ilo, ihi, rwork( ibal ), n, vs, ldvs,
405  $ ierr )
406  END IF
407 *
408  IF( scalea ) THEN
409 *
410 * Undo scaling for the Schur form of A
411 *
412  CALL zlascl( 'U', 0, 0, cscale, anrm, n, n, a, lda, ierr )
413  CALL zcopy( n, a, lda+1, w, 1 )
414  END IF
415 *
416  work( 1 ) = maxwrk
417  RETURN
418 *
419 * End of ZGEES
420 *
421  END
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
ZGEBAL
Definition: zgebal.f:162
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
Definition: zgehrd.f:167
subroutine zgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
ZGEBAK
Definition: zgebak.f:131
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 m...
Definition: zgees.f:197
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:143
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
ZHSEQR
Definition: zhseqr.f:299
subroutine ztrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, SEP, WORK, LWORK, INFO)
ZTRSEN
Definition: ztrsen.f:264
subroutine zunghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGHR
Definition: zunghr.f:126