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