001:       SUBROUTINE SGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM,
002:      $                   WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK,
003:      $                   IWORK, LIWORK, 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, LIWORK, LWORK, N, SDIM
013:       REAL               RCONDE, RCONDV
014: *     ..
015: *     .. Array Arguments ..
016:       LOGICAL            BWORK( * )
017:       INTEGER            IWORK( * )
018:       REAL               A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
019:      $                   WR( * )
020: *     ..
021: *     .. Function Arguments ..
022:       LOGICAL            SELECT
023:       EXTERNAL           SELECT
024: *     ..
025: *
026: *  Purpose
027: *  =======
028: *
029: *  SGEESX computes for an N-by-N real nonsymmetric matrix A, the
030: *  eigenvalues, the real Schur form T, and, optionally, the matrix of
031: *  Schur vectors Z.  This gives the Schur factorization A = Z*T*(Z**T).
032: *
033: *  Optionally, it also orders the eigenvalues on the diagonal of the
034: *  real Schur form so that selected eigenvalues are at the top left;
035: *  computes a reciprocal condition number for the average of the
036: *  selected eigenvalues (RCONDE); and computes a reciprocal condition
037: *  number for the right invariant subspace corresponding to the
038: *  selected eigenvalues (RCONDV).  The leading columns of Z form an
039: *  orthonormal basis for this invariant subspace.
040: *
041: *  For further explanation of the reciprocal condition numbers RCONDE
042: *  and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
043: *  these quantities are called s and sep respectively).
044: *
045: *  A real matrix is in real Schur form if it is upper quasi-triangular
046: *  with 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in
047: *  the form
048: *            [  a  b  ]
049: *            [  c  a  ]
050: *
051: *  where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
052: *
053: *  Arguments
054: *  =========
055: *
056: *  JOBVS   (input) CHARACTER*1
057: *          = 'N': Schur vectors are not computed;
058: *          = 'V': Schur vectors are computed.
059: *
060: *  SORT    (input) CHARACTER*1
061: *          Specifies whether or not to order the eigenvalues on the
062: *          diagonal of the Schur form.
063: *          = 'N': Eigenvalues are not ordered;
064: *          = 'S': Eigenvalues are ordered (see SELECT).
065: *
066: *  SELECT  (external procedure) LOGICAL FUNCTION of two REAL arguments
067: *          SELECT must be declared EXTERNAL in the calling subroutine.
068: *          If SORT = 'S', SELECT is used to select eigenvalues to sort
069: *          to the top left of the Schur form.
070: *          If SORT = 'N', SELECT is not referenced.
071: *          An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
072: *          SELECT(WR(j),WI(j)) is true; i.e., if either one of a
073: *          complex conjugate pair of eigenvalues is selected, then both
074: *          are.  Note that a selected complex eigenvalue may no longer
075: *          satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
076: *          ordering may change the value of complex eigenvalues
077: *          (especially if the eigenvalue is ill-conditioned); in this
078: *          case INFO may be set to N+3 (see INFO below).
079: *
080: *  SENSE   (input) CHARACTER*1
081: *          Determines which reciprocal condition numbers are computed.
082: *          = 'N': None are computed;
083: *          = 'E': Computed for average of selected eigenvalues only;
084: *          = 'V': Computed for selected right invariant subspace only;
085: *          = 'B': Computed for both.
086: *          If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
087: *
088: *  N       (input) INTEGER
089: *          The order of the matrix A. N >= 0.
090: *
091: *  A       (input/output) REAL array, dimension (LDA, N)
092: *          On entry, the N-by-N matrix A.
093: *          On exit, A is overwritten by its real Schur form T.
094: *
095: *  LDA     (input) INTEGER
096: *          The leading dimension of the array A.  LDA >= max(1,N).
097: *
098: *  SDIM    (output) INTEGER
099: *          If SORT = 'N', SDIM = 0.
100: *          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
101: *                         for which SELECT is true. (Complex conjugate
102: *                         pairs for which SELECT is true for either
103: *                         eigenvalue count as 2.)
104: *
105: *  WR      (output) REAL array, dimension (N)
106: *  WI      (output) REAL array, dimension (N)
107: *          WR and WI contain the real and imaginary parts, respectively,
108: *          of the computed eigenvalues, in the same order that they
109: *          appear on the diagonal of the output Schur form T.  Complex
110: *          conjugate pairs of eigenvalues appear consecutively with the
111: *          eigenvalue having the positive imaginary part first.
112: *
113: *  VS      (output) REAL array, dimension (LDVS,N)
114: *          If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
115: *          vectors.
116: *          If JOBVS = 'N', VS is not referenced.
117: *
118: *  LDVS    (input) INTEGER
119: *          The leading dimension of the array VS.  LDVS >= 1, and if
120: *          JOBVS = 'V', LDVS >= N.
121: *
122: *  RCONDE  (output) REAL
123: *          If SENSE = 'E' or 'B', RCONDE contains the reciprocal
124: *          condition number for the average of the selected eigenvalues.
125: *          Not referenced if SENSE = 'N' or 'V'.
126: *
127: *  RCONDV  (output) REAL
128: *          If SENSE = 'V' or 'B', RCONDV contains the reciprocal
129: *          condition number for the selected right invariant subspace.
130: *          Not referenced if SENSE = 'N' or 'E'.
131: *
132: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
133: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
134: *
135: *  LWORK   (input) INTEGER
136: *          The dimension of the array WORK.  LWORK >= max(1,3*N).
137: *          Also, if SENSE = 'E' or 'V' or 'B',
138: *          LWORK >= N+2*SDIM*(N-SDIM), where SDIM is the number of
139: *          selected eigenvalues computed by this routine.  Note that
140: *          N+2*SDIM*(N-SDIM) <= N+N*N/2. Note also that an error is only
141: *          returned if LWORK < max(1,3*N), but if SENSE = 'E' or 'V' or
142: *          'B' this may not be large enough.
143: *          For good performance, LWORK must generally be larger.
144: *
145: *          If LWORK = -1, then a workspace query is assumed; the routine
146: *          only calculates upper bounds on the optimal sizes of the
147: *          arrays WORK and IWORK, returns these values as the first
148: *          entries of the WORK and IWORK arrays, and no error messages
149: *          related to LWORK or LIWORK are issued by XERBLA.
150: *
151: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
152: *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
153: *
154: *  LIWORK  (input) INTEGER
155: *          The dimension of the array IWORK.
156: *          LIWORK >= 1; if SENSE = 'V' or 'B', LIWORK >= SDIM*(N-SDIM).
157: *          Note that SDIM*(N-SDIM) <= N*N/4. Note also that an error is
158: *          only returned if LIWORK < 1, but if SENSE = 'V' or 'B' this
159: *          may not be large enough.
160: *
161: *          If LIWORK = -1, then a workspace query is assumed; the
162: *          routine only calculates upper bounds on the optimal sizes of
163: *          the arrays WORK and IWORK, returns these values as the first
164: *          entries of the WORK and IWORK arrays, and no error messages
165: *          related to LWORK or LIWORK are issued by XERBLA.
166: *
167: *  BWORK   (workspace) LOGICAL array, dimension (N)
168: *          Not referenced if SORT = 'N'.
169: *
170: *  INFO    (output) INTEGER
171: *          = 0: successful exit
172: *          < 0: if INFO = -i, the i-th argument had an illegal value.
173: *          > 0: if INFO = i, and i is
174: *             <= N: the QR algorithm failed to compute all the
175: *                   eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
176: *                   contain those eigenvalues which have converged; if
177: *                   JOBVS = 'V', VS contains the transformation which
178: *                   reduces A to its partially converged Schur form.
179: *             = N+1: the eigenvalues could not be reordered because some
180: *                   eigenvalues were too close to separate (the problem
181: *                   is very ill-conditioned);
182: *             = N+2: after reordering, roundoff changed values of some
183: *                   complex eigenvalues so that leading eigenvalues in
184: *                   the Schur form no longer satisfy SELECT=.TRUE.  This
185: *                   could also be caused by underflow due to scaling.
186: *
187: *  =====================================================================
188: *
189: *     .. Parameters ..
190:       REAL               ZERO, ONE
191:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
192: *     ..
193: *     .. Local Scalars ..
194:       LOGICAL            CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTSB,
195:      $                   WANTSE, WANTSN, WANTST, WANTSV, WANTVS
196:       INTEGER            HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
197:      $                   IHI, ILO, INXT, IP, ITAU, IWRK, LWRK, LIWRK,
198:      $                   MAXWRK, MINWRK
199:       REAL               ANRM, BIGNUM, CSCALE, EPS, SMLNUM
200: *     ..
201: *     .. Local Arrays ..
202:       REAL               DUM( 1 )
203: *     ..
204: *     .. External Subroutines ..
205:       EXTERNAL           SCOPY, SGEBAK, SGEBAL, SGEHRD, SHSEQR, SLABAD,
206:      $                   SLACPY, SLASCL, SORGHR, SSWAP, STRSEN, XERBLA
207: *     ..
208: *     .. External Functions ..
209:       LOGICAL            LSAME
210:       INTEGER            ILAENV
211:       REAL               SLAMCH, SLANGE
212:       EXTERNAL           LSAME, ILAENV, SLAMCH, SLANGE
213: *     ..
214: *     .. Intrinsic Functions ..
215:       INTRINSIC          MAX, SQRT
216: *     ..
217: *     .. Executable Statements ..
218: *
219: *     Test the input arguments
220: *
221:       INFO = 0
222:       WANTVS = LSAME( JOBVS, 'V' )
223:       WANTST = LSAME( SORT, 'S' )
224:       WANTSN = LSAME( SENSE, 'N' )
225:       WANTSE = LSAME( SENSE, 'E' )
226:       WANTSV = LSAME( SENSE, 'V' )
227:       WANTSB = LSAME( SENSE, 'B' )
228:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
229:       IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
230:          INFO = -1
231:       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
232:          INFO = -2
233:       ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
234:      $         ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
235:          INFO = -4
236:       ELSE IF( N.LT.0 ) THEN
237:          INFO = -5
238:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
239:          INFO = -7
240:       ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
241:          INFO = -12
242:       END IF
243: *
244: *     Compute workspace
245: *      (Note: Comments in the code beginning "RWorkspace:" describe the
246: *       minimal amount of real workspace needed at that point in the
247: *       code, as well as the preferred amount for good performance.
248: *       IWorkspace refers to integer workspace.
249: *       NB refers to the optimal block size for the immediately
250: *       following subroutine, as returned by ILAENV.
251: *       HSWORK refers to the workspace preferred by SHSEQR, as
252: *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
253: *       the worst case.
254: *       If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
255: *       depends on SDIM, which is computed by the routine STRSEN later
256: *       in the code.)
257: *
258:       IF( INFO.EQ.0 ) THEN
259:          LIWRK = 1
260:          IF( N.EQ.0 ) THEN
261:             MINWRK = 1
262:             LWRK = 1
263:          ELSE
264:             MAXWRK = 2*N + N*ILAENV( 1, 'SGEHRD', ' ', N, 1, N, 0 )
265:             MINWRK = 3*N
266: *
267:             CALL SHSEQR( 'S', JOBVS, N, 1, N, A, LDA, WR, WI, VS, LDVS,
268:      $             WORK, -1, IEVAL )
269:             HSWORK = WORK( 1 )
270: *
271:             IF( .NOT.WANTVS ) THEN
272:                MAXWRK = MAX( MAXWRK, N + HSWORK )
273:             ELSE
274:                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
275:      $                       'SORGHR', ' ', N, 1, N, -1 ) )
276:                MAXWRK = MAX( MAXWRK, N + HSWORK )
277:             END IF
278:             LWRK = MAXWRK
279:             IF( .NOT.WANTSN )
280:      $         LWRK = MAX( LWRK, N + ( N*N )/2 )
281:             IF( WANTSV .OR. WANTSB )
282:      $         LIWRK = ( N*N )/4
283:          END IF
284:          IWORK( 1 ) = LIWRK
285:          WORK( 1 ) = LWRK
286: *
287:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
288:             INFO = -16
289:          ELSE IF( LIWORK.LT.1 .AND. .NOT.LQUERY ) THEN
290:             INFO = -18
291:          END IF
292:       END IF
293: *
294:       IF( INFO.NE.0 ) THEN
295:          CALL XERBLA( 'SGEESX', -INFO )
296:          RETURN
297:       END IF
298: *
299: *     Quick return if possible
300: *
301:       IF( N.EQ.0 ) THEN
302:          SDIM = 0
303:          RETURN
304:       END IF
305: *
306: *     Get machine constants
307: *
308:       EPS = SLAMCH( 'P' )
309:       SMLNUM = SLAMCH( 'S' )
310:       BIGNUM = ONE / SMLNUM
311:       CALL SLABAD( SMLNUM, BIGNUM )
312:       SMLNUM = SQRT( SMLNUM ) / EPS
313:       BIGNUM = ONE / SMLNUM
314: *
315: *     Scale A if max element outside range [SMLNUM,BIGNUM]
316: *
317:       ANRM = SLANGE( 'M', N, N, A, LDA, DUM )
318:       SCALEA = .FALSE.
319:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
320:          SCALEA = .TRUE.
321:          CSCALE = SMLNUM
322:       ELSE IF( ANRM.GT.BIGNUM ) THEN
323:          SCALEA = .TRUE.
324:          CSCALE = BIGNUM
325:       END IF
326:       IF( SCALEA )
327:      $   CALL SLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
328: *
329: *     Permute the matrix to make it more nearly triangular
330: *     (RWorkspace: need N)
331: *
332:       IBAL = 1
333:       CALL SGEBAL( 'P', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
334: *
335: *     Reduce to upper Hessenberg form
336: *     (RWorkspace: need 3*N, prefer 2*N+N*NB)
337: *
338:       ITAU = N + IBAL
339:       IWRK = N + ITAU
340:       CALL SGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
341:      $             LWORK-IWRK+1, IERR )
342: *
343:       IF( WANTVS ) THEN
344: *
345: *        Copy Householder vectors to VS
346: *
347:          CALL SLACPY( 'L', N, N, A, LDA, VS, LDVS )
348: *
349: *        Generate orthogonal matrix in VS
350: *        (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
351: *
352:          CALL SORGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
353:      $                LWORK-IWRK+1, IERR )
354:       END IF
355: *
356:       SDIM = 0
357: *
358: *     Perform QR iteration, accumulating Schur vectors in VS if desired
359: *     (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
360: *
361:       IWRK = ITAU
362:       CALL SHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, WR, WI, VS, LDVS,
363:      $             WORK( IWRK ), LWORK-IWRK+1, IEVAL )
364:       IF( IEVAL.GT.0 )
365:      $   INFO = IEVAL
366: *
367: *     Sort eigenvalues if desired
368: *
369:       IF( WANTST .AND. INFO.EQ.0 ) THEN
370:          IF( SCALEA ) THEN
371:             CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WR, N, IERR )
372:             CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WI, N, IERR )
373:          END IF
374:          DO 10 I = 1, N
375:             BWORK( I ) = SELECT( WR( I ), WI( I ) )
376:    10    CONTINUE
377: *
378: *        Reorder eigenvalues, transform Schur vectors, and compute
379: *        reciprocal condition numbers
380: *        (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
381: *                     otherwise, need N )
382: *        (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
383: *                     otherwise, need 0 )
384: *
385:          CALL STRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, WR, WI,
386:      $                SDIM, RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1,
387:      $                IWORK, LIWORK, ICOND )
388:          IF( .NOT.WANTSN )
389:      $      MAXWRK = MAX( MAXWRK, N+2*SDIM*( N-SDIM ) )
390:          IF( ICOND.EQ.-15 ) THEN
391: *
392: *           Not enough real workspace
393: *
394:             INFO = -16
395:          ELSE IF( ICOND.EQ.-17 ) THEN
396: *
397: *           Not enough integer workspace
398: *
399:             INFO = -18
400:          ELSE IF( ICOND.GT.0 ) THEN
401: *
402: *           STRSEN failed to reorder or to restore standard Schur form
403: *
404:             INFO = ICOND + N
405:          END IF
406:       END IF
407: *
408:       IF( WANTVS ) THEN
409: *
410: *        Undo balancing
411: *        (RWorkspace: need N)
412: *
413:          CALL SGEBAK( 'P', 'R', N, ILO, IHI, WORK( IBAL ), N, VS, LDVS,
414:      $                IERR )
415:       END IF
416: *
417:       IF( SCALEA ) THEN
418: *
419: *        Undo scaling for the Schur form of A
420: *
421:          CALL SLASCL( 'H', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
422:          CALL SCOPY( N, A, LDA+1, WR, 1 )
423:          IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN
424:             DUM( 1 ) = RCONDV
425:             CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
426:             RCONDV = DUM( 1 )
427:          END IF
428:          IF( CSCALE.EQ.SMLNUM ) THEN
429: *
430: *           If scaling back towards underflow, adjust WI if an
431: *           offdiagonal element of a 2-by-2 block in the Schur form
432: *           underflows.
433: *
434:             IF( IEVAL.GT.0 ) THEN
435:                I1 = IEVAL + 1
436:                I2 = IHI - 1
437:                CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
438:      $                      IERR )
439:             ELSE IF( WANTST ) THEN
440:                I1 = 1
441:                I2 = N - 1
442:             ELSE
443:                I1 = ILO
444:                I2 = IHI - 1
445:             END IF
446:             INXT = I1 - 1
447:             DO 20 I = I1, I2
448:                IF( I.LT.INXT )
449:      $            GO TO 20
450:                IF( WI( I ).EQ.ZERO ) THEN
451:                   INXT = I + 1
452:                ELSE
453:                   IF( A( I+1, I ).EQ.ZERO ) THEN
454:                      WI( I ) = ZERO
455:                      WI( I+1 ) = ZERO
456:                   ELSE IF( A( I+1, I ).NE.ZERO .AND. A( I, I+1 ).EQ.
457:      $                     ZERO ) THEN
458:                      WI( I ) = ZERO
459:                      WI( I+1 ) = ZERO
460:                      IF( I.GT.1 )
461:      $                  CALL SSWAP( I-1, A( 1, I ), 1, A( 1, I+1 ), 1 )
462:                      IF( N.GT.I+1 )
463:      $                  CALL SSWAP( N-I-1, A( I, I+2 ), LDA,
464:      $                              A( I+1, I+2 ), LDA )
465:                      CALL SSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
466:                      A( I, I+1 ) = A( I+1, I )
467:                      A( I+1, I ) = ZERO
468:                   END IF
469:                   INXT = I + 2
470:                END IF
471:    20       CONTINUE
472:          END IF
473:          CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-IEVAL, 1,
474:      $                WI( IEVAL+1 ), MAX( N-IEVAL, 1 ), IERR )
475:       END IF
476: *
477:       IF( WANTST .AND. INFO.EQ.0 ) THEN
478: *
479: *        Check if reordering successful
480: *
481:          LASTSL = .TRUE.
482:          LST2SL = .TRUE.
483:          SDIM = 0
484:          IP = 0
485:          DO 30 I = 1, N
486:             CURSL = SELECT( WR( I ), WI( I ) )
487:             IF( WI( I ).EQ.ZERO ) THEN
488:                IF( CURSL )
489:      $            SDIM = SDIM + 1
490:                IP = 0
491:                IF( CURSL .AND. .NOT.LASTSL )
492:      $            INFO = N + 2
493:             ELSE
494:                IF( IP.EQ.1 ) THEN
495: *
496: *                 Last eigenvalue of conjugate pair
497: *
498:                   CURSL = CURSL .OR. LASTSL
499:                   LASTSL = CURSL
500:                   IF( CURSL )
501:      $               SDIM = SDIM + 2
502:                   IP = -1
503:                   IF( CURSL .AND. .NOT.LST2SL )
504:      $               INFO = N + 2
505:                ELSE
506: *
507: *                 First eigenvalue of conjugate pair
508: *
509:                   IP = 1
510:                END IF
511:             END IF
512:             LST2SL = LASTSL
513:             LASTSL = CURSL
514:    30    CONTINUE
515:       END IF
516: *
517:       WORK( 1 ) = MAXWRK
518:       IF( WANTSV .OR. WANTSB ) THEN
519:          IWORK( 1 ) = SDIM*(N-SDIM)
520:       ELSE
521:          IWORK( 1 ) = 1
522:       END IF
523: *
524:       RETURN
525: *
526: *     End of SGEESX
527: *
528:       END
529: