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