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