001:       SUBROUTINE SGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
002:      $                  LDVR, WORK, LWORK, 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          JOBVL, JOBVR
011:       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
015:      $                   WI( * ), WORK( * ), WR( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  SGEEV computes for an N-by-N real nonsymmetric matrix A, the
022: *  eigenvalues and, optionally, the left and/or right eigenvectors.
023: *
024: *  The right eigenvector v(j) of A satisfies
025: *                   A * v(j) = lambda(j) * v(j)
026: *  where lambda(j) is its eigenvalue.
027: *  The left eigenvector u(j) of A satisfies
028: *                u(j)**H * A = lambda(j) * u(j)**H
029: *  where u(j)**H denotes the conjugate transpose of u(j).
030: *
031: *  The computed eigenvectors are normalized to have Euclidean norm
032: *  equal to 1 and largest component real.
033: *
034: *  Arguments
035: *  =========
036: *
037: *  JOBVL   (input) CHARACTER*1
038: *          = 'N': left eigenvectors of A are not computed;
039: *          = 'V': left eigenvectors of A are computed.
040: *
041: *  JOBVR   (input) CHARACTER*1
042: *          = 'N': right eigenvectors of A are not computed;
043: *          = 'V': right eigenvectors of A are computed.
044: *
045: *  N       (input) INTEGER
046: *          The order of the matrix A. N >= 0.
047: *
048: *  A       (input/output) REAL array, dimension (LDA,N)
049: *          On entry, the N-by-N matrix A.
050: *          On exit, A has been overwritten.
051: *
052: *  LDA     (input) INTEGER
053: *          The leading dimension of the array A.  LDA >= max(1,N).
054: *
055: *  WR      (output) REAL array, dimension (N)
056: *  WI      (output) REAL array, dimension (N)
057: *          WR and WI contain the real and imaginary parts,
058: *          respectively, of the computed eigenvalues.  Complex
059: *          conjugate pairs of eigenvalues appear consecutively
060: *          with the eigenvalue having the positive imaginary part
061: *          first.
062: *
063: *  VL      (output) REAL array, dimension (LDVL,N)
064: *          If JOBVL = 'V', the left eigenvectors u(j) are stored one
065: *          after another in the columns of VL, in the same order
066: *          as their eigenvalues.
067: *          If JOBVL = 'N', VL is not referenced.
068: *          If the j-th eigenvalue is real, then u(j) = VL(:,j),
069: *          the j-th column of VL.
070: *          If the j-th and (j+1)-st eigenvalues form a complex
071: *          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
072: *          u(j+1) = VL(:,j) - i*VL(:,j+1).
073: *
074: *  LDVL    (input) INTEGER
075: *          The leading dimension of the array VL.  LDVL >= 1; if
076: *          JOBVL = 'V', LDVL >= N.
077: *
078: *  VR      (output) REAL array, dimension (LDVR,N)
079: *          If JOBVR = 'V', the right eigenvectors v(j) are stored one
080: *          after another in the columns of VR, in the same order
081: *          as their eigenvalues.
082: *          If JOBVR = 'N', VR is not referenced.
083: *          If the j-th eigenvalue is real, then v(j) = VR(:,j),
084: *          the j-th column of VR.
085: *          If the j-th and (j+1)-st eigenvalues form a complex
086: *          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
087: *          v(j+1) = VR(:,j) - i*VR(:,j+1).
088: *
089: *  LDVR    (input) INTEGER
090: *          The leading dimension of the array VR.  LDVR >= 1; if
091: *          JOBVR = 'V', LDVR >= N.
092: *
093: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
094: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
095: *
096: *  LWORK   (input) INTEGER
097: *          The dimension of the array WORK.  LWORK >= max(1,3*N), and
098: *          if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.  For good
099: *          performance, LWORK must generally be larger.
100: *
101: *          If LWORK = -1, then a workspace query is assumed; the routine
102: *          only calculates the optimal size of the WORK array, returns
103: *          this value as the first entry of the WORK array, and no error
104: *          message related to LWORK is issued by XERBLA.
105: *
106: *  INFO    (output) INTEGER
107: *          = 0:  successful exit
108: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
109: *          > 0:  if INFO = i, the QR algorithm failed to compute all the
110: *                eigenvalues, and no eigenvectors have been computed;
111: *                elements i+1:N of WR and WI contain eigenvalues which
112: *                have converged.
113: *
114: *  =====================================================================
115: *
116: *     .. Parameters ..
117:       REAL               ZERO, ONE
118:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
119: *     ..
120: *     .. Local Scalars ..
121:       LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR
122:       CHARACTER          SIDE
123:       INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
124:      $                   MAXWRK, MINWRK, NOUT
125:       REAL               ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
126:      $                   SN
127: *     ..
128: *     .. Local Arrays ..
129:       LOGICAL            SELECT( 1 )
130:       REAL               DUM( 1 )
131: *     ..
132: *     .. External Subroutines ..
133:       EXTERNAL           SGEBAK, SGEBAL, SGEHRD, SHSEQR, SLABAD, SLACPY,
134:      $                   SLARTG, SLASCL, SORGHR, SROT, SSCAL, STREVC,
135:      $                   XERBLA
136: *     ..
137: *     .. External Functions ..
138:       LOGICAL            LSAME
139:       INTEGER            ILAENV, ISAMAX
140:       REAL               SLAMCH, SLANGE, SLAPY2, SNRM2
141:       EXTERNAL           LSAME, ILAENV, ISAMAX, SLAMCH, SLANGE, SLAPY2,
142:      $                   SNRM2
143: *     ..
144: *     .. Intrinsic Functions ..
145:       INTRINSIC          MAX, SQRT
146: *     ..
147: *     .. Executable Statements ..
148: *
149: *     Test the input arguments
150: *
151:       INFO = 0
152:       LQUERY = ( LWORK.EQ.-1 )
153:       WANTVL = LSAME( JOBVL, 'V' )
154:       WANTVR = LSAME( JOBVR, 'V' )
155:       IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
156:          INFO = -1
157:       ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
158:          INFO = -2
159:       ELSE IF( N.LT.0 ) THEN
160:          INFO = -3
161:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
162:          INFO = -5
163:       ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
164:          INFO = -9
165:       ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
166:          INFO = -11
167:       END IF
168: *
169: *     Compute workspace
170: *      (Note: Comments in the code beginning "Workspace:" describe the
171: *       minimal amount of workspace needed at that point in the code,
172: *       as well as the preferred amount for good performance.
173: *       NB refers to the optimal block size for the immediately
174: *       following subroutine, as returned by ILAENV.
175: *       HSWORK refers to the workspace preferred by SHSEQR, as
176: *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
177: *       the worst case.)
178: *
179:       IF( INFO.EQ.0 ) THEN
180:          IF( N.EQ.0 ) THEN
181:             MINWRK = 1
182:             MAXWRK = 1
183:          ELSE
184:             MAXWRK = 2*N + N*ILAENV( 1, 'SGEHRD', ' ', N, 1, N, 0 )
185:             IF( WANTVL ) THEN
186:                MINWRK = 4*N
187:                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
188:      $                       'SORGHR', ' ', N, 1, N, -1 ) )
189:                CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
190:      $                WORK, -1, INFO )
191:                HSWORK = WORK( 1 )
192:                MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
193:                MAXWRK = MAX( MAXWRK, 4*N )
194:             ELSE IF( WANTVR ) THEN
195:                MINWRK = 4*N
196:                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
197:      $                       'SORGHR', ' ', N, 1, N, -1 ) )
198:                CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
199:      $                WORK, -1, INFO )
200:                HSWORK = WORK( 1 )
201:                MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
202:                MAXWRK = MAX( MAXWRK, 4*N )
203:             ELSE 
204:                MINWRK = 3*N
205:                CALL SHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
206:      $                WORK, -1, INFO )
207:                HSWORK = WORK( 1 )
208:                MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
209:             END IF
210:             MAXWRK = MAX( MAXWRK, MINWRK )
211:          END IF
212:          WORK( 1 ) = MAXWRK
213: *
214:          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
215:             INFO = -13
216:          END IF
217:       END IF
218: *
219:       IF( INFO.NE.0 ) THEN
220:          CALL XERBLA( 'SGEEV ', -INFO )
221:          RETURN
222:       ELSE IF( LQUERY ) THEN
223:          RETURN
224:       END IF
225: *
226: *     Quick return if possible
227: *
228:       IF( N.EQ.0 )
229:      $   RETURN
230: *
231: *     Get machine constants
232: *
233:       EPS = SLAMCH( 'P' )
234:       SMLNUM = SLAMCH( 'S' )
235:       BIGNUM = ONE / SMLNUM
236:       CALL SLABAD( SMLNUM, BIGNUM )
237:       SMLNUM = SQRT( SMLNUM ) / EPS
238:       BIGNUM = ONE / SMLNUM
239: *
240: *     Scale A if max element outside range [SMLNUM,BIGNUM]
241: *
242:       ANRM = SLANGE( 'M', N, N, A, LDA, DUM )
243:       SCALEA = .FALSE.
244:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
245:          SCALEA = .TRUE.
246:          CSCALE = SMLNUM
247:       ELSE IF( ANRM.GT.BIGNUM ) THEN
248:          SCALEA = .TRUE.
249:          CSCALE = BIGNUM
250:       END IF
251:       IF( SCALEA )
252:      $   CALL SLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
253: *
254: *     Balance the matrix
255: *     (Workspace: need N)
256: *
257:       IBAL = 1
258:       CALL SGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
259: *
260: *     Reduce to upper Hessenberg form
261: *     (Workspace: need 3*N, prefer 2*N+N*NB)
262: *
263:       ITAU = IBAL + N
264:       IWRK = ITAU + N
265:       CALL SGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
266:      $             LWORK-IWRK+1, IERR )
267: *
268:       IF( WANTVL ) THEN
269: *
270: *        Want left eigenvectors
271: *        Copy Householder vectors to VL
272: *
273:          SIDE = 'L'
274:          CALL SLACPY( 'L', N, N, A, LDA, VL, LDVL )
275: *
276: *        Generate orthogonal matrix in VL
277: *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
278: *
279:          CALL SORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
280:      $                LWORK-IWRK+1, IERR )
281: *
282: *        Perform QR iteration, accumulating Schur vectors in VL
283: *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
284: *
285:          IWRK = ITAU
286:          CALL SHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
287:      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
288: *
289:          IF( WANTVR ) THEN
290: *
291: *           Want left and right eigenvectors
292: *           Copy Schur vectors to VR
293: *
294:             SIDE = 'B'
295:             CALL SLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
296:          END IF
297: *
298:       ELSE IF( WANTVR ) THEN
299: *
300: *        Want right eigenvectors
301: *        Copy Householder vectors to VR
302: *
303:          SIDE = 'R'
304:          CALL SLACPY( 'L', N, N, A, LDA, VR, LDVR )
305: *
306: *        Generate orthogonal matrix in VR
307: *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
308: *
309:          CALL SORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
310:      $                LWORK-IWRK+1, IERR )
311: *
312: *        Perform QR iteration, accumulating Schur vectors in VR
313: *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
314: *
315:          IWRK = ITAU
316:          CALL SHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
317:      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
318: *
319:       ELSE
320: *
321: *        Compute eigenvalues only
322: *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
323: *
324:          IWRK = ITAU
325:          CALL SHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
326:      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
327:       END IF
328: *
329: *     If INFO > 0 from SHSEQR, then quit
330: *
331:       IF( INFO.GT.0 )
332:      $   GO TO 50
333: *
334:       IF( WANTVL .OR. WANTVR ) THEN
335: *
336: *        Compute left and/or right eigenvectors
337: *        (Workspace: need 4*N)
338: *
339:          CALL STREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
340:      $                N, NOUT, WORK( IWRK ), IERR )
341:       END IF
342: *
343:       IF( WANTVL ) THEN
344: *
345: *        Undo balancing of left eigenvectors
346: *        (Workspace: need N)
347: *
348:          CALL SGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
349:      $                IERR )
350: *
351: *        Normalize left eigenvectors and make largest component real
352: *
353:          DO 20 I = 1, N
354:             IF( WI( I ).EQ.ZERO ) THEN
355:                SCL = ONE / SNRM2( N, VL( 1, I ), 1 )
356:                CALL SSCAL( N, SCL, VL( 1, I ), 1 )
357:             ELSE IF( WI( I ).GT.ZERO ) THEN
358:                SCL = ONE / SLAPY2( SNRM2( N, VL( 1, I ), 1 ),
359:      $               SNRM2( N, VL( 1, I+1 ), 1 ) )
360:                CALL SSCAL( N, SCL, VL( 1, I ), 1 )
361:                CALL SSCAL( N, SCL, VL( 1, I+1 ), 1 )
362:                DO 10 K = 1, N
363:                   WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
364:    10          CONTINUE
365:                K = ISAMAX( N, WORK( IWRK ), 1 )
366:                CALL SLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
367:                CALL SROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
368:                VL( K, I+1 ) = ZERO
369:             END IF
370:    20    CONTINUE
371:       END IF
372: *
373:       IF( WANTVR ) THEN
374: *
375: *        Undo balancing of right eigenvectors
376: *        (Workspace: need N)
377: *
378:          CALL SGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
379:      $                IERR )
380: *
381: *        Normalize right eigenvectors and make largest component real
382: *
383:          DO 40 I = 1, N
384:             IF( WI( I ).EQ.ZERO ) THEN
385:                SCL = ONE / SNRM2( N, VR( 1, I ), 1 )
386:                CALL SSCAL( N, SCL, VR( 1, I ), 1 )
387:             ELSE IF( WI( I ).GT.ZERO ) THEN
388:                SCL = ONE / SLAPY2( SNRM2( N, VR( 1, I ), 1 ),
389:      $               SNRM2( N, VR( 1, I+1 ), 1 ) )
390:                CALL SSCAL( N, SCL, VR( 1, I ), 1 )
391:                CALL SSCAL( N, SCL, VR( 1, I+1 ), 1 )
392:                DO 30 K = 1, N
393:                   WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
394:    30          CONTINUE
395:                K = ISAMAX( N, WORK( IWRK ), 1 )
396:                CALL SLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
397:                CALL SROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
398:                VR( K, I+1 ) = ZERO
399:             END IF
400:    40    CONTINUE
401:       END IF
402: *
403: *     Undo scaling if necessary
404: *
405:    50 CONTINUE
406:       IF( SCALEA ) THEN
407:          CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
408:      $                MAX( N-INFO, 1 ), IERR )
409:          CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
410:      $                MAX( N-INFO, 1 ), IERR )
411:          IF( INFO.GT.0 ) THEN
412:             CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
413:      $                   IERR )
414:             CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
415:      $                   IERR )
416:          END IF
417:       END IF
418: *
419:       WORK( 1 ) = MAXWRK
420:       RETURN
421: *
422: *     End of SGEEV
423: *
424:       END
425: