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