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