001:       SUBROUTINE SGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
002:      $                  BETA, VL, LDVL, VR, 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, LDB, LDVL, LDVR, LWORK, N
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
015:      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
016:      $                   VR( LDVR, * ), WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  This routine is deprecated and has been replaced by routine SGGEV.
023: *
024: *  SGEGV computes the eigenvalues and, optionally, the left and/or right
025: *  eigenvectors of a real matrix pair (A,B).
026: *  Given two square matrices A and B,
027: *  the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
028: *  eigenvalues lambda and corresponding (non-zero) eigenvectors x such
029: *  that
030: *
031: *     A*x = lambda*B*x.
032: *
033: *  An alternate form is to find the eigenvalues mu and corresponding
034: *  eigenvectors y such that
035: *
036: *     mu*A*y = B*y.
037: *
038: *  These two forms are equivalent with mu = 1/lambda and x = y if
039: *  neither lambda nor mu is zero.  In order to deal with the case that
040: *  lambda or mu is zero or small, two values alpha and beta are returned
041: *  for each eigenvalue, such that lambda = alpha/beta and
042: *  mu = beta/alpha.
043: *
044: *  The vectors x and y in the above equations are right eigenvectors of
045: *  the matrix pair (A,B).  Vectors u and v satisfying
046: *
047: *     u**H*A = lambda*u**H*B  or  mu*v**H*A = v**H*B
048: *
049: *  are left eigenvectors of (A,B).
050: *
051: *  Note: this routine performs "full balancing" on A and B -- see
052: *  "Further Details", below.
053: *
054: *  Arguments
055: *  =========
056: *
057: *  JOBVL   (input) CHARACTER*1
058: *          = 'N':  do not compute the left generalized eigenvectors;
059: *          = 'V':  compute the left generalized eigenvectors (returned
060: *                  in VL).
061: *
062: *  JOBVR   (input) CHARACTER*1
063: *          = 'N':  do not compute the right generalized eigenvectors;
064: *          = 'V':  compute the right generalized eigenvectors (returned
065: *                  in VR).
066: *
067: *  N       (input) INTEGER
068: *          The order of the matrices A, B, VL, and VR.  N >= 0.
069: *
070: *  A       (input/output) REAL array, dimension (LDA, N)
071: *          On entry, the matrix A.
072: *          If JOBVL = 'V' or JOBVR = 'V', then on exit A
073: *          contains the real Schur form of A from the generalized Schur
074: *          factorization of the pair (A,B) after balancing.
075: *          If no eigenvectors were computed, then only the diagonal
076: *          blocks from the Schur form will be correct.  See SGGHRD and
077: *          SHGEQZ for details.
078: *
079: *  LDA     (input) INTEGER
080: *          The leading dimension of A.  LDA >= max(1,N).
081: *
082: *  B       (input/output) REAL array, dimension (LDB, N)
083: *          On entry, the matrix B.
084: *          If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
085: *          upper triangular matrix obtained from B in the generalized
086: *          Schur factorization of the pair (A,B) after balancing.
087: *          If no eigenvectors were computed, then only those elements of
088: *          B corresponding to the diagonal blocks from the Schur form of
089: *          A will be correct.  See SGGHRD and SHGEQZ for details.
090: *
091: *  LDB     (input) INTEGER
092: *          The leading dimension of B.  LDB >= max(1,N).
093: *
094: *  ALPHAR  (output) REAL array, dimension (N)
095: *          The real parts of each scalar alpha defining an eigenvalue of
096: *          GNEP.
097: *
098: *  ALPHAI  (output) REAL array, dimension (N)
099: *          The imaginary parts of each scalar alpha defining an
100: *          eigenvalue of GNEP.  If ALPHAI(j) is zero, then the j-th
101: *          eigenvalue is real; if positive, then the j-th and
102: *          (j+1)-st eigenvalues are a complex conjugate pair, with
103: *          ALPHAI(j+1) = -ALPHAI(j).
104: *
105: *  BETA    (output) REAL array, dimension (N)
106: *          The scalars beta that define the eigenvalues of GNEP.
107: *          
108: *          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
109: *          beta = BETA(j) represent the j-th eigenvalue of the matrix
110: *          pair (A,B), in one of the forms lambda = alpha/beta or
111: *          mu = beta/alpha.  Since either lambda or mu may overflow,
112: *          they should not, in general, be computed.
113: *
114: *  VL      (output) REAL array, dimension (LDVL,N)
115: *          If JOBVL = 'V', the left eigenvectors u(j) are stored
116: *          in the columns of VL, in the same order as their eigenvalues.
117: *          If the j-th eigenvalue is real, then u(j) = VL(:,j).
118: *          If the j-th and (j+1)-st eigenvalues form a complex conjugate
119: *          pair, then
120: *             u(j) = VL(:,j) + i*VL(:,j+1)
121: *          and
122: *            u(j+1) = VL(:,j) - i*VL(:,j+1).
123: *
124: *          Each eigenvector is scaled so that its largest component has
125: *          abs(real part) + abs(imag. part) = 1, except for eigenvectors
126: *          corresponding to an eigenvalue with alpha = beta = 0, which
127: *          are set to zero.
128: *          Not referenced if JOBVL = 'N'.
129: *
130: *  LDVL    (input) INTEGER
131: *          The leading dimension of the matrix VL. LDVL >= 1, and
132: *          if JOBVL = 'V', LDVL >= N.
133: *
134: *  VR      (output) REAL array, dimension (LDVR,N)
135: *          If JOBVR = 'V', the right eigenvectors x(j) are stored
136: *          in the columns of VR, in the same order as their eigenvalues.
137: *          If the j-th eigenvalue is real, then x(j) = VR(:,j).
138: *          If the j-th and (j+1)-st eigenvalues form a complex conjugate
139: *          pair, then
140: *            x(j) = VR(:,j) + i*VR(:,j+1)
141: *          and
142: *            x(j+1) = VR(:,j) - i*VR(:,j+1).
143: *
144: *          Each eigenvector is scaled so that its largest component has
145: *          abs(real part) + abs(imag. part) = 1, except for eigenvalues
146: *          corresponding to an eigenvalue with alpha = beta = 0, which
147: *          are set to zero.
148: *          Not referenced if JOBVR = 'N'.
149: *
150: *  LDVR    (input) INTEGER
151: *          The leading dimension of the matrix VR. LDVR >= 1, and
152: *          if JOBVR = 'V', LDVR >= N.
153: *
154: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
155: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
156: *
157: *  LWORK   (input) INTEGER
158: *          The dimension of the array WORK.  LWORK >= max(1,8*N).
159: *          For good performance, LWORK must generally be larger.
160: *          To compute the optimal value of LWORK, call ILAENV to get
161: *          blocksizes (for SGEQRF, SORMQR, and SORGQR.)  Then compute:
162: *          NB  -- MAX of the blocksizes for SGEQRF, SORMQR, and SORGQR;
163: *          The optimal LWORK is:
164: *              2*N + MAX( 6*N, N*(NB+1) ).
165: *
166: *          If LWORK = -1, then a workspace query is assumed; the routine
167: *          only calculates the optimal size of the WORK array, returns
168: *          this value as the first entry of the WORK array, and no error
169: *          message related to LWORK is issued by XERBLA.
170: *
171: *  INFO    (output) INTEGER
172: *          = 0:  successful exit
173: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
174: *          = 1,...,N:
175: *                The QZ iteration failed.  No eigenvectors have been
176: *                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
177: *                should be correct for j=INFO+1,...,N.
178: *          > N:  errors that usually indicate LAPACK problems:
179: *                =N+1: error return from SGGBAL
180: *                =N+2: error return from SGEQRF
181: *                =N+3: error return from SORMQR
182: *                =N+4: error return from SORGQR
183: *                =N+5: error return from SGGHRD
184: *                =N+6: error return from SHGEQZ (other than failed
185: *                                                iteration)
186: *                =N+7: error return from STGEVC
187: *                =N+8: error return from SGGBAK (computing VL)
188: *                =N+9: error return from SGGBAK (computing VR)
189: *                =N+10: error return from SLASCL (various calls)
190: *
191: *  Further Details
192: *  ===============
193: *
194: *  Balancing
195: *  ---------
196: *
197: *  This driver calls SGGBAL to both permute and scale rows and columns
198: *  of A and B.  The permutations PL and PR are chosen so that PL*A*PR
199: *  and PL*B*R will be upper triangular except for the diagonal blocks
200: *  A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
201: *  possible.  The diagonal scaling matrices DL and DR are chosen so
202: *  that the pair  DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
203: *  one (except for the elements that start out zero.)
204: *
205: *  After the eigenvalues and eigenvectors of the balanced matrices
206: *  have been computed, SGGBAK transforms the eigenvectors back to what
207: *  they would have been (in perfect arithmetic) if they had not been
208: *  balanced.
209: *
210: *  Contents of A and B on Exit
211: *  -------- -- - --- - -- ----
212: *
213: *  If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
214: *  both), then on exit the arrays A and B will contain the real Schur
215: *  form[*] of the "balanced" versions of A and B.  If no eigenvectors
216: *  are computed, then only the diagonal blocks will be correct.
217: *
218: *  [*] See SHGEQZ, SGEGS, or read the book "Matrix Computations",
219: *      by Golub & van Loan, pub. by Johns Hopkins U. Press.
220: *
221: *  =====================================================================
222: *
223: *     .. Parameters ..
224:       REAL               ZERO, ONE
225:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
226: *     ..
227: *     .. Local Scalars ..
228:       LOGICAL            ILIMIT, ILV, ILVL, ILVR, LQUERY
229:       CHARACTER          CHTEMP
230:       INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
231:      $                   IN, IRIGHT, IROWS, ITAU, IWORK, JC, JR, LOPT,
232:      $                   LWKMIN, LWKOPT, NB, NB1, NB2, NB3
233:       REAL               ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
234:      $                   BNRM1, BNRM2, EPS, ONEPLS, SAFMAX, SAFMIN,
235:      $                   SALFAI, SALFAR, SBETA, SCALE, TEMP
236: *     ..
237: *     .. Local Arrays ..
238:       LOGICAL            LDUMMA( 1 )
239: *     ..
240: *     .. External Subroutines ..
241:       EXTERNAL           SGEQRF, SGGBAK, SGGBAL, SGGHRD, SHGEQZ, SLACPY,
242:      $                   SLASCL, SLASET, SORGQR, SORMQR, STGEVC, XERBLA
243: *     ..
244: *     .. External Functions ..
245:       LOGICAL            LSAME
246:       INTEGER            ILAENV
247:       REAL               SLAMCH, SLANGE
248:       EXTERNAL           ILAENV, LSAME, SLAMCH, SLANGE
249: *     ..
250: *     .. Intrinsic Functions ..
251:       INTRINSIC          ABS, INT, MAX
252: *     ..
253: *     .. Executable Statements ..
254: *
255: *     Decode the input arguments
256: *
257:       IF( LSAME( JOBVL, 'N' ) ) THEN
258:          IJOBVL = 1
259:          ILVL = .FALSE.
260:       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
261:          IJOBVL = 2
262:          ILVL = .TRUE.
263:       ELSE
264:          IJOBVL = -1
265:          ILVL = .FALSE.
266:       END IF
267: *
268:       IF( LSAME( JOBVR, 'N' ) ) THEN
269:          IJOBVR = 1
270:          ILVR = .FALSE.
271:       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
272:          IJOBVR = 2
273:          ILVR = .TRUE.
274:       ELSE
275:          IJOBVR = -1
276:          ILVR = .FALSE.
277:       END IF
278:       ILV = ILVL .OR. ILVR
279: *
280: *     Test the input arguments
281: *
282:       LWKMIN = MAX( 8*N, 1 )
283:       LWKOPT = LWKMIN
284:       WORK( 1 ) = LWKOPT
285:       LQUERY = ( LWORK.EQ.-1 )
286:       INFO = 0
287:       IF( IJOBVL.LE.0 ) THEN
288:          INFO = -1
289:       ELSE IF( IJOBVR.LE.0 ) THEN
290:          INFO = -2
291:       ELSE IF( N.LT.0 ) THEN
292:          INFO = -3
293:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
294:          INFO = -5
295:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
296:          INFO = -7
297:       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
298:          INFO = -12
299:       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
300:          INFO = -14
301:       ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
302:          INFO = -16
303:       END IF
304: *
305:       IF( INFO.EQ.0 ) THEN
306:          NB1 = ILAENV( 1, 'SGEQRF', ' ', N, N, -1, -1 )
307:          NB2 = ILAENV( 1, 'SORMQR', ' ', N, N, N, -1 )
308:          NB3 = ILAENV( 1, 'SORGQR', ' ', N, N, N, -1 )
309:          NB = MAX( NB1, NB2, NB3 )
310:          LOPT = 2*N + MAX( 6*N, N*(NB+1) )
311:          WORK( 1 ) = LOPT
312:       END IF
313: *
314:       IF( INFO.NE.0 ) THEN
315:          CALL XERBLA( 'SGEGV ', -INFO )
316:          RETURN
317:       ELSE IF( LQUERY ) THEN
318:          RETURN
319:       END IF
320: *
321: *     Quick return if possible
322: *
323:       IF( N.EQ.0 )
324:      $   RETURN
325: *
326: *     Get machine constants
327: *
328:       EPS = SLAMCH( 'E' )*SLAMCH( 'B' )
329:       SAFMIN = SLAMCH( 'S' )
330:       SAFMIN = SAFMIN + SAFMIN
331:       SAFMAX = ONE / SAFMIN
332:       ONEPLS = ONE + ( 4*EPS )
333: *
334: *     Scale A
335: *
336:       ANRM = SLANGE( 'M', N, N, A, LDA, WORK )
337:       ANRM1 = ANRM
338:       ANRM2 = ONE
339:       IF( ANRM.LT.ONE ) THEN
340:          IF( SAFMAX*ANRM.LT.ONE ) THEN
341:             ANRM1 = SAFMIN
342:             ANRM2 = SAFMAX*ANRM
343:          END IF
344:       END IF
345: *
346:       IF( ANRM.GT.ZERO ) THEN
347:          CALL SLASCL( 'G', -1, -1, ANRM, ONE, N, N, A, LDA, IINFO )
348:          IF( IINFO.NE.0 ) THEN
349:             INFO = N + 10
350:             RETURN
351:          END IF
352:       END IF
353: *
354: *     Scale B
355: *
356:       BNRM = SLANGE( 'M', N, N, B, LDB, WORK )
357:       BNRM1 = BNRM
358:       BNRM2 = ONE
359:       IF( BNRM.LT.ONE ) THEN
360:          IF( SAFMAX*BNRM.LT.ONE ) THEN
361:             BNRM1 = SAFMIN
362:             BNRM2 = SAFMAX*BNRM
363:          END IF
364:       END IF
365: *
366:       IF( BNRM.GT.ZERO ) THEN
367:          CALL SLASCL( 'G', -1, -1, BNRM, ONE, N, N, B, LDB, IINFO )
368:          IF( IINFO.NE.0 ) THEN
369:             INFO = N + 10
370:             RETURN
371:          END IF
372:       END IF
373: *
374: *     Permute the matrix to make it more nearly triangular
375: *     Workspace layout:  (8*N words -- "work" requires 6*N words)
376: *        left_permutation, right_permutation, work...
377: *
378:       ILEFT = 1
379:       IRIGHT = N + 1
380:       IWORK = IRIGHT + N
381:       CALL SGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
382:      $             WORK( IRIGHT ), WORK( IWORK ), IINFO )
383:       IF( IINFO.NE.0 ) THEN
384:          INFO = N + 1
385:          GO TO 120
386:       END IF
387: *
388: *     Reduce B to triangular form, and initialize VL and/or VR
389: *     Workspace layout:  ("work..." must have at least N words)
390: *        left_permutation, right_permutation, tau, work...
391: *
392:       IROWS = IHI + 1 - ILO
393:       IF( ILV ) THEN
394:          ICOLS = N + 1 - ILO
395:       ELSE
396:          ICOLS = IROWS
397:       END IF
398:       ITAU = IWORK
399:       IWORK = ITAU + IROWS
400:       CALL SGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
401:      $             WORK( IWORK ), LWORK+1-IWORK, IINFO )
402:       IF( IINFO.GE.0 )
403:      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
404:       IF( IINFO.NE.0 ) THEN
405:          INFO = N + 2
406:          GO TO 120
407:       END IF
408: *
409:       CALL SORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
410:      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
411:      $             LWORK+1-IWORK, IINFO )
412:       IF( IINFO.GE.0 )
413:      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
414:       IF( IINFO.NE.0 ) THEN
415:          INFO = N + 3
416:          GO TO 120
417:       END IF
418: *
419:       IF( ILVL ) THEN
420:          CALL SLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
421:          CALL SLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
422:      $                VL( ILO+1, ILO ), LDVL )
423:          CALL SORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
424:      $                WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
425:      $                IINFO )
426:          IF( IINFO.GE.0 )
427:      $      LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
428:          IF( IINFO.NE.0 ) THEN
429:             INFO = N + 4
430:             GO TO 120
431:          END IF
432:       END IF
433: *
434:       IF( ILVR )
435:      $   CALL SLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
436: *
437: *     Reduce to generalized Hessenberg form
438: *
439:       IF( ILV ) THEN
440: *
441: *        Eigenvectors requested -- work on whole matrix.
442: *
443:          CALL SGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
444:      $                LDVL, VR, LDVR, IINFO )
445:       ELSE
446:          CALL SGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
447:      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IINFO )
448:       END IF
449:       IF( IINFO.NE.0 ) THEN
450:          INFO = N + 5
451:          GO TO 120
452:       END IF
453: *
454: *     Perform QZ algorithm
455: *     Workspace layout:  ("work..." must have at least 1 word)
456: *        left_permutation, right_permutation, work...
457: *
458:       IWORK = ITAU
459:       IF( ILV ) THEN
460:          CHTEMP = 'S'
461:       ELSE
462:          CHTEMP = 'E'
463:       END IF
464:       CALL SHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
465:      $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
466:      $             WORK( IWORK ), LWORK+1-IWORK, IINFO )
467:       IF( IINFO.GE.0 )
468:      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
469:       IF( IINFO.NE.0 ) THEN
470:          IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
471:             INFO = IINFO
472:          ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
473:             INFO = IINFO - N
474:          ELSE
475:             INFO = N + 6
476:          END IF
477:          GO TO 120
478:       END IF
479: *
480:       IF( ILV ) THEN
481: *
482: *        Compute Eigenvectors  (STGEVC requires 6*N words of workspace)
483: *
484:          IF( ILVL ) THEN
485:             IF( ILVR ) THEN
486:                CHTEMP = 'B'
487:             ELSE
488:                CHTEMP = 'L'
489:             END IF
490:          ELSE
491:             CHTEMP = 'R'
492:          END IF
493: *
494:          CALL STGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
495:      $                VR, LDVR, N, IN, WORK( IWORK ), IINFO )
496:          IF( IINFO.NE.0 ) THEN
497:             INFO = N + 7
498:             GO TO 120
499:          END IF
500: *
501: *        Undo balancing on VL and VR, rescale
502: *
503:          IF( ILVL ) THEN
504:             CALL SGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
505:      $                   WORK( IRIGHT ), N, VL, LDVL, IINFO )
506:             IF( IINFO.NE.0 ) THEN
507:                INFO = N + 8
508:                GO TO 120
509:             END IF
510:             DO 50 JC = 1, N
511:                IF( ALPHAI( JC ).LT.ZERO )
512:      $            GO TO 50
513:                TEMP = ZERO
514:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
515:                   DO 10 JR = 1, N
516:                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
517:    10             CONTINUE
518:                ELSE
519:                   DO 20 JR = 1, N
520:                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
521:      $                      ABS( VL( JR, JC+1 ) ) )
522:    20             CONTINUE
523:                END IF
524:                IF( TEMP.LT.SAFMIN )
525:      $            GO TO 50
526:                TEMP = ONE / TEMP
527:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
528:                   DO 30 JR = 1, N
529:                      VL( JR, JC ) = VL( JR, JC )*TEMP
530:    30             CONTINUE
531:                ELSE
532:                   DO 40 JR = 1, N
533:                      VL( JR, JC ) = VL( JR, JC )*TEMP
534:                      VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
535:    40             CONTINUE
536:                END IF
537:    50       CONTINUE
538:          END IF
539:          IF( ILVR ) THEN
540:             CALL SGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
541:      $                   WORK( IRIGHT ), N, VR, LDVR, IINFO )
542:             IF( IINFO.NE.0 ) THEN
543:                INFO = N + 9
544:                GO TO 120
545:             END IF
546:             DO 100 JC = 1, N
547:                IF( ALPHAI( JC ).LT.ZERO )
548:      $            GO TO 100
549:                TEMP = ZERO
550:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
551:                   DO 60 JR = 1, N
552:                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
553:    60             CONTINUE
554:                ELSE
555:                   DO 70 JR = 1, N
556:                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
557:      $                      ABS( VR( JR, JC+1 ) ) )
558:    70             CONTINUE
559:                END IF
560:                IF( TEMP.LT.SAFMIN )
561:      $            GO TO 100
562:                TEMP = ONE / TEMP
563:                IF( ALPHAI( JC ).EQ.ZERO ) THEN
564:                   DO 80 JR = 1, N
565:                      VR( JR, JC ) = VR( JR, JC )*TEMP
566:    80             CONTINUE
567:                ELSE
568:                   DO 90 JR = 1, N
569:                      VR( JR, JC ) = VR( JR, JC )*TEMP
570:                      VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
571:    90             CONTINUE
572:                END IF
573:   100       CONTINUE
574:          END IF
575: *
576: *        End of eigenvector calculation
577: *
578:       END IF
579: *
580: *     Undo scaling in alpha, beta
581: *
582: *     Note: this does not give the alpha and beta for the unscaled
583: *     problem.
584: *
585: *     Un-scaling is limited to avoid underflow in alpha and beta
586: *     if they are significant.
587: *
588:       DO 110 JC = 1, N
589:          ABSAR = ABS( ALPHAR( JC ) )
590:          ABSAI = ABS( ALPHAI( JC ) )
591:          ABSB = ABS( BETA( JC ) )
592:          SALFAR = ANRM*ALPHAR( JC )
593:          SALFAI = ANRM*ALPHAI( JC )
594:          SBETA = BNRM*BETA( JC )
595:          ILIMIT = .FALSE.
596:          SCALE = ONE
597: *
598: *        Check for significant underflow in ALPHAI
599: *
600:          IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE.
601:      $       MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN
602:             ILIMIT = .TRUE.
603:             SCALE = ( ONEPLS*SAFMIN / ANRM1 ) /
604:      $              MAX( ONEPLS*SAFMIN, ANRM2*ABSAI )
605: *
606:          ELSE IF( SALFAI.EQ.ZERO ) THEN
607: *
608: *           If insignificant underflow in ALPHAI, then make the
609: *           conjugate eigenvalue real.
610: *
611:             IF( ALPHAI( JC ).LT.ZERO .AND. JC.GT.1 ) THEN
612:                ALPHAI( JC-1 ) = ZERO
613:             ELSE IF( ALPHAI( JC ).GT.ZERO .AND. JC.LT.N ) THEN
614:                ALPHAI( JC+1 ) = ZERO
615:             END IF
616:          END IF
617: *
618: *        Check for significant underflow in ALPHAR
619: *
620:          IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE.
621:      $       MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN
622:             ILIMIT = .TRUE.
623:             SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / ANRM1 ) /
624:      $              MAX( ONEPLS*SAFMIN, ANRM2*ABSAR ) )
625:          END IF
626: *
627: *        Check for significant underflow in BETA
628: *
629:          IF( ABS( SBETA ).LT.SAFMIN .AND. ABSB.GE.
630:      $       MAX( SAFMIN, EPS*ABSAR, EPS*ABSAI ) ) THEN
631:             ILIMIT = .TRUE.
632:             SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / BNRM1 ) /
633:      $              MAX( ONEPLS*SAFMIN, BNRM2*ABSB ) )
634:          END IF
635: *
636: *        Check for possible overflow when limiting scaling
637: *
638:          IF( ILIMIT ) THEN
639:             TEMP = ( SCALE*SAFMIN )*MAX( ABS( SALFAR ), ABS( SALFAI ),
640:      $             ABS( SBETA ) )
641:             IF( TEMP.GT.ONE )
642:      $         SCALE = SCALE / TEMP
643:             IF( SCALE.LT.ONE )
644:      $         ILIMIT = .FALSE.
645:          END IF
646: *
647: *        Recompute un-scaled ALPHAR, ALPHAI, BETA if necessary.
648: *
649:          IF( ILIMIT ) THEN
650:             SALFAR = ( SCALE*ALPHAR( JC ) )*ANRM
651:             SALFAI = ( SCALE*ALPHAI( JC ) )*ANRM
652:             SBETA = ( SCALE*BETA( JC ) )*BNRM
653:          END IF
654:          ALPHAR( JC ) = SALFAR
655:          ALPHAI( JC ) = SALFAI
656:          BETA( JC ) = SBETA
657:   110 CONTINUE
658: *
659:   120 CONTINUE
660:       WORK( 1 ) = LWKOPT
661: *
662:       RETURN
663: *
664: *     End of SGEGV
665: *
666:       END
667: