001:       SUBROUTINE ZGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
002:      $                  LDVS, WORK, LWORK, RWORK, BWORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          JOBVS, SORT
011:       INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
012: *     ..
013: *     .. Array Arguments ..
014:       LOGICAL            BWORK( * )
015:       DOUBLE PRECISION   RWORK( * )
016:       COMPLEX*16         A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
017: *     ..
018: *     .. Function Arguments ..
019:       LOGICAL            SELECT
020:       EXTERNAL           SELECT
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: *
026: *  ZGEES computes for an N-by-N complex nonsymmetric matrix A, the
027: *  eigenvalues, the Schur form T, and, optionally, the matrix of Schur
028: *  vectors Z.  This gives the Schur factorization A = Z*T*(Z**H).
029: *
030: *  Optionally, it also orders the eigenvalues on the diagonal of the
031: *  Schur form so that selected eigenvalues are at the top left.
032: *  The leading columns of Z then form an orthonormal basis for the
033: *  invariant subspace corresponding to the selected eigenvalues.
034: *
035: *  A complex matrix is in Schur form if it is upper triangular.
036: *
037: *  Arguments
038: *  =========
039: *
040: *  JOBVS   (input) CHARACTER*1
041: *          = 'N': Schur vectors are not computed;
042: *          = 'V': Schur vectors are computed.
043: *
044: *  SORT    (input) CHARACTER*1
045: *          Specifies whether or not to order the eigenvalues on the
046: *          diagonal of the Schur form.
047: *          = 'N': Eigenvalues are not ordered:
048: *          = 'S': Eigenvalues are ordered (see SELECT).
049: *
050: *  SELECT  (external procedure) LOGICAL FUNCTION of one COMPLEX*16 argument
051: *          SELECT must be declared EXTERNAL in the calling subroutine.
052: *          If SORT = 'S', SELECT is used to select eigenvalues to order
053: *          to the top left of the Schur form.
054: *          IF SORT = 'N', SELECT is not referenced.
055: *          The eigenvalue W(j) is selected if SELECT(W(j)) is true.
056: *
057: *  N       (input) INTEGER
058: *          The order of the matrix A. N >= 0.
059: *
060: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
061: *          On entry, the N-by-N matrix A.
062: *          On exit, A has been overwritten by its Schur form T.
063: *
064: *  LDA     (input) INTEGER
065: *          The leading dimension of the array A.  LDA >= max(1,N).
066: *
067: *  SDIM    (output) INTEGER
068: *          If SORT = 'N', SDIM = 0.
069: *          If SORT = 'S', SDIM = number of eigenvalues for which
070: *                         SELECT is true.
071: *
072: *  W       (output) COMPLEX*16 array, dimension (N)
073: *          W contains the computed eigenvalues, in the same order that
074: *          they appear on the diagonal of the output Schur form T.
075: *
076: *  VS      (output) COMPLEX*16 array, dimension (LDVS,N)
077: *          If JOBVS = 'V', VS contains the unitary matrix Z of Schur
078: *          vectors.
079: *          If JOBVS = 'N', VS is not referenced.
080: *
081: *  LDVS    (input) INTEGER
082: *          The leading dimension of the array VS.  LDVS >= 1; if
083: *          JOBVS = 'V', LDVS >= N.
084: *
085: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
086: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
087: *
088: *  LWORK   (input) INTEGER
089: *          The dimension of the array WORK.  LWORK >= max(1,2*N).
090: *          For good performance, LWORK must generally be larger.
091: *
092: *          If LWORK = -1, then a workspace query is assumed; the routine
093: *          only calculates the optimal size of the WORK array, returns
094: *          this value as the first entry of the WORK array, and no error
095: *          message related to LWORK is issued by XERBLA.
096: *
097: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
098: *
099: *  BWORK   (workspace) LOGICAL array, dimension (N)
100: *          Not referenced if SORT = 'N'.
101: *
102: *  INFO    (output) INTEGER
103: *          = 0: successful exit
104: *          < 0: if INFO = -i, the i-th argument had an illegal value.
105: *          > 0: if INFO = i, and i is
106: *               <= N:  the QR algorithm failed to compute all the
107: *                      eigenvalues; elements 1:ILO-1 and i+1:N of W
108: *                      contain those eigenvalues which have converged;
109: *                      if JOBVS = 'V', VS contains the matrix which
110: *                      reduces A to its partially converged Schur form.
111: *               = N+1: the eigenvalues could not be reordered because
112: *                      some eigenvalues were too close to separate (the
113: *                      problem is very ill-conditioned);
114: *               = N+2: after reordering, roundoff changed values of
115: *                      some complex eigenvalues so that leading
116: *                      eigenvalues in the Schur form no longer satisfy
117: *                      SELECT = .TRUE..  This could also be caused by
118: *                      underflow due to scaling.
119: *
120: *  =====================================================================
121: *
122: *     .. Parameters ..
123:       DOUBLE PRECISION   ZERO, ONE
124:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
125: *     ..
126: *     .. Local Scalars ..
127:       LOGICAL            LQUERY, SCALEA, WANTST, WANTVS
128:       INTEGER            HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
129:      $                   ITAU, IWRK, MAXWRK, MINWRK
130:       DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
131: *     ..
132: *     .. Local Arrays ..
133:       DOUBLE PRECISION   DUM( 1 )
134: *     ..
135: *     .. External Subroutines ..
136:       EXTERNAL           DLABAD, XERBLA, ZCOPY, ZGEBAK, ZGEBAL, ZGEHRD,
137:      $                   ZHSEQR, ZLACPY, ZLASCL, ZTRSEN, ZUNGHR
138: *     ..
139: *     .. External Functions ..
140:       LOGICAL            LSAME
141:       INTEGER            ILAENV
142:       DOUBLE PRECISION   DLAMCH, ZLANGE
143:       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
144: *     ..
145: *     .. Intrinsic Functions ..
146:       INTRINSIC          MAX, SQRT
147: *     ..
148: *     .. Executable Statements ..
149: *
150: *     Test the input arguments
151: *
152:       INFO = 0
153:       LQUERY = ( LWORK.EQ.-1 )
154:       WANTVS = LSAME( JOBVS, 'V' )
155:       WANTST = LSAME( SORT, 'S' )
156:       IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
157:          INFO = -1
158:       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
159:          INFO = -2
160:       ELSE IF( N.LT.0 ) THEN
161:          INFO = -4
162:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
163:          INFO = -6
164:       ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
165:          INFO = -10
166:       END IF
167: *
168: *     Compute workspace
169: *      (Note: Comments in the code beginning "Workspace:" describe the
170: *       minimal amount of workspace needed at that point in the code,
171: *       as well as the preferred amount for good performance.
172: *       CWorkspace refers to complex workspace, and RWorkspace to real
173: *       workspace. NB refers to the optimal block size for the
174: *       immediately following subroutine, as returned by ILAENV.
175: *       HSWORK refers to the workspace preferred by ZHSEQR, as
176: *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
177: *       the worst case.)
178: *
179:       IF( INFO.EQ.0 ) THEN
180:          IF( N.EQ.0 ) THEN
181:             MINWRK = 1
182:             MAXWRK = 1
183:          ELSE
184:             MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
185:             MINWRK = 2*N
186: *
187:             CALL ZHSEQR( 'S', JOBVS, N, 1, N, A, LDA, W, VS, LDVS,
188:      $             WORK, -1, IEVAL )
189:             HSWORK = WORK( 1 )
190: *
191:             IF( .NOT.WANTVS ) THEN
192:                MAXWRK = MAX( MAXWRK, HSWORK )
193:             ELSE
194:                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
195:      $                       ' ', N, 1, N, -1 ) )
196:                MAXWRK = MAX( MAXWRK, HSWORK )
197:             END IF
198:          END IF
199:          WORK( 1 ) = MAXWRK
200: *
201:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
202:             INFO = -12
203:          END IF
204:       END IF
205: *
206:       IF( INFO.NE.0 ) THEN
207:          CALL XERBLA( 'ZGEES ', -INFO )
208:          RETURN
209:       ELSE IF( LQUERY ) THEN
210:          RETURN
211:       END IF
212: *
213: *     Quick return if possible
214: *
215:       IF( N.EQ.0 ) THEN
216:          SDIM = 0
217:          RETURN
218:       END IF
219: *
220: *     Get machine constants
221: *
222:       EPS = DLAMCH( 'P' )
223:       SMLNUM = DLAMCH( 'S' )
224:       BIGNUM = ONE / SMLNUM
225:       CALL DLABAD( SMLNUM, BIGNUM )
226:       SMLNUM = SQRT( SMLNUM ) / EPS
227:       BIGNUM = ONE / SMLNUM
228: *
229: *     Scale A if max element outside range [SMLNUM,BIGNUM]
230: *
231:       ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
232:       SCALEA = .FALSE.
233:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
234:          SCALEA = .TRUE.
235:          CSCALE = SMLNUM
236:       ELSE IF( ANRM.GT.BIGNUM ) THEN
237:          SCALEA = .TRUE.
238:          CSCALE = BIGNUM
239:       END IF
240:       IF( SCALEA )
241:      $   CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
242: *
243: *     Permute the matrix to make it more nearly triangular
244: *     (CWorkspace: none)
245: *     (RWorkspace: need N)
246: *
247:       IBAL = 1
248:       CALL ZGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
249: *
250: *     Reduce to upper Hessenberg form
251: *     (CWorkspace: need 2*N, prefer N+N*NB)
252: *     (RWorkspace: none)
253: *
254:       ITAU = 1
255:       IWRK = N + ITAU
256:       CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
257:      $             LWORK-IWRK+1, IERR )
258: *
259:       IF( WANTVS ) THEN
260: *
261: *        Copy Householder vectors to VS
262: *
263:          CALL ZLACPY( 'L', N, N, A, LDA, VS, LDVS )
264: *
265: *        Generate unitary matrix in VS
266: *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
267: *        (RWorkspace: none)
268: *
269:          CALL ZUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
270:      $                LWORK-IWRK+1, IERR )
271:       END IF
272: *
273:       SDIM = 0
274: *
275: *     Perform QR iteration, accumulating Schur vectors in VS if desired
276: *     (CWorkspace: need 1, prefer HSWORK (see comments) )
277: *     (RWorkspace: none)
278: *
279:       IWRK = ITAU
280:       CALL ZHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS,
281:      $             WORK( IWRK ), LWORK-IWRK+1, IEVAL )
282:       IF( IEVAL.GT.0 )
283:      $   INFO = IEVAL
284: *
285: *     Sort eigenvalues if desired
286: *
287:       IF( WANTST .AND. INFO.EQ.0 ) THEN
288:          IF( SCALEA )
289:      $      CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR )
290:          DO 10 I = 1, N
291:             BWORK( I ) = SELECT( W( I ) )
292:    10    CONTINUE
293: *
294: *        Reorder eigenvalues and transform Schur vectors
295: *        (CWorkspace: none)
296: *        (RWorkspace: none)
297: *
298:          CALL ZTRSEN( 'N', JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM,
299:      $                S, SEP, WORK( IWRK ), LWORK-IWRK+1, ICOND )
300:       END IF
301: *
302:       IF( WANTVS ) THEN
303: *
304: *        Undo balancing
305: *        (CWorkspace: none)
306: *        (RWorkspace: need N)
307: *
308:          CALL ZGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS,
309:      $                IERR )
310:       END IF
311: *
312:       IF( SCALEA ) THEN
313: *
314: *        Undo scaling for the Schur form of A
315: *
316:          CALL ZLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
317:          CALL ZCOPY( N, A, LDA+1, W, 1 )
318:       END IF
319: *
320:       WORK( 1 ) = MAXWRK
321:       RETURN
322: *
323: *     End of ZGEES
324: *
325:       END
326: