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