001:       SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
002:      $                   LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
003: *
004: *  -- LAPACK 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          HOWMNY, SIDE
011:       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
012: *     ..
013: *     .. Array Arguments ..
014:       LOGICAL            SELECT( * )
015:       REAL               RWORK( * )
016:       COMPLEX            P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
017:      $                   VR( LDVR, * ), WORK( * )
018: *     ..
019: *
020: *
021: *  Purpose
022: *  =======
023: *
024: *  CTGEVC computes some or all of the right and/or left eigenvectors of
025: *  a pair of complex matrices (S,P), where S and P are upper triangular.
026: *  Matrix pairs of this type are produced by the generalized Schur
027: *  factorization of a complex matrix pair (A,B):
028: *  
029: *     A = Q*S*Z**H,  B = Q*P*Z**H
030: *  
031: *  as computed by CGGHRD + CHGEQZ.
032: *  
033: *  The right eigenvector x and the left eigenvector y of (S,P)
034: *  corresponding to an eigenvalue w are defined by:
035: *  
036: *     S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
037: *  
038: *  where y**H denotes the conjugate tranpose of y.
039: *  The eigenvalues are not input to this routine, but are computed
040: *  directly from the diagonal elements of S and P.
041: *  
042: *  This routine returns the matrices X and/or Y of right and left
043: *  eigenvectors of (S,P), or the products Z*X and/or Q*Y,
044: *  where Z and Q are input matrices.
045: *  If Q and Z are the unitary factors from the generalized Schur
046: *  factorization of a matrix pair (A,B), then Z*X and Q*Y
047: *  are the matrices of right and left eigenvectors of (A,B).
048: *
049: *  Arguments
050: *  =========
051: *
052: *  SIDE    (input) CHARACTER*1
053: *          = 'R': compute right eigenvectors only;
054: *          = 'L': compute left eigenvectors only;
055: *          = 'B': compute both right and left eigenvectors.
056: *
057: *  HOWMNY  (input) CHARACTER*1
058: *          = 'A': compute all right and/or left eigenvectors;
059: *          = 'B': compute all right and/or left eigenvectors,
060: *                 backtransformed by the matrices in VR and/or VL;
061: *          = 'S': compute selected right and/or left eigenvectors,
062: *                 specified by the logical array SELECT.
063: *
064: *  SELECT  (input) LOGICAL array, dimension (N)
065: *          If HOWMNY='S', SELECT specifies the eigenvectors to be
066: *          computed.  The eigenvector corresponding to the j-th
067: *          eigenvalue is computed if SELECT(j) = .TRUE..
068: *          Not referenced if HOWMNY = 'A' or 'B'.
069: *
070: *  N       (input) INTEGER
071: *          The order of the matrices S and P.  N >= 0.
072: *
073: *  S       (input) COMPLEX array, dimension (LDS,N)
074: *          The upper triangular matrix S from a generalized Schur
075: *          factorization, as computed by CHGEQZ.
076: *
077: *  LDS     (input) INTEGER
078: *          The leading dimension of array S.  LDS >= max(1,N).
079: *
080: *  P       (input) COMPLEX array, dimension (LDP,N)
081: *          The upper triangular matrix P from a generalized Schur
082: *          factorization, as computed by CHGEQZ.  P must have real
083: *          diagonal elements.
084: *
085: *  LDP     (input) INTEGER
086: *          The leading dimension of array P.  LDP >= max(1,N).
087: *
088: *  VL      (input/output) COMPLEX array, dimension (LDVL,MM)
089: *          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
090: *          contain an N-by-N matrix Q (usually the unitary matrix Q
091: *          of left Schur vectors returned by CHGEQZ).
092: *          On exit, if SIDE = 'L' or 'B', VL contains:
093: *          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
094: *          if HOWMNY = 'B', the matrix Q*Y;
095: *          if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
096: *                      SELECT, stored consecutively in the columns of
097: *                      VL, in the same order as their eigenvalues.
098: *          Not referenced if SIDE = 'R'.
099: *
100: *  LDVL    (input) INTEGER
101: *          The leading dimension of array VL.  LDVL >= 1, and if
102: *          SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
103: *
104: *  VR      (input/output) COMPLEX array, dimension (LDVR,MM)
105: *          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
106: *          contain an N-by-N matrix Q (usually the unitary matrix Z
107: *          of right Schur vectors returned by CHGEQZ).
108: *          On exit, if SIDE = 'R' or 'B', VR contains:
109: *          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
110: *          if HOWMNY = 'B', the matrix Z*X;
111: *          if HOWMNY = 'S', the right eigenvectors of (S,P) specified by
112: *                      SELECT, stored consecutively in the columns of
113: *                      VR, in the same order as their eigenvalues.
114: *          Not referenced if SIDE = 'L'.
115: *
116: *  LDVR    (input) INTEGER
117: *          The leading dimension of the array VR.  LDVR >= 1, and if
118: *          SIDE = 'R' or 'B', LDVR >= N.
119: *
120: *  MM      (input) INTEGER
121: *          The number of columns in the arrays VL and/or VR. MM >= M.
122: *
123: *  M       (output) INTEGER
124: *          The number of columns in the arrays VL and/or VR actually
125: *          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
126: *          is set to N.  Each selected eigenvector occupies one column.
127: *
128: *  WORK    (workspace) COMPLEX array, dimension (2*N)
129: *
130: *  RWORK   (workspace) REAL array, dimension (2*N)
131: *
132: *  INFO    (output) INTEGER
133: *          = 0:  successful exit.
134: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
135: *
136: *  =====================================================================
137: *
138: *     .. Parameters ..
139:       REAL               ZERO, ONE
140:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
141:       COMPLEX            CZERO, CONE
142:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
143:      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
144: *     ..
145: *     .. Local Scalars ..
146:       LOGICAL            COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
147:      $                   LSA, LSB
148:       INTEGER            I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
149:      $                   J, JE, JR
150:       REAL               ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
151:      $                   BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
152:      $                   SCALE, SMALL, TEMP, ULP, XMAX
153:       COMPLEX            BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
154: *     ..
155: *     .. External Functions ..
156:       LOGICAL            LSAME
157:       REAL               SLAMCH
158:       COMPLEX            CLADIV
159:       EXTERNAL           LSAME, SLAMCH, CLADIV
160: *     ..
161: *     .. External Subroutines ..
162:       EXTERNAL           CGEMV, SLABAD, XERBLA
163: *     ..
164: *     .. Intrinsic Functions ..
165:       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL
166: *     ..
167: *     .. Statement Functions ..
168:       REAL               ABS1
169: *     ..
170: *     .. Statement Function definitions ..
171:       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
172: *     ..
173: *     .. Executable Statements ..
174: *
175: *     Decode and Test the input parameters
176: *
177:       IF( LSAME( HOWMNY, 'A' ) ) THEN
178:          IHWMNY = 1
179:          ILALL = .TRUE.
180:          ILBACK = .FALSE.
181:       ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
182:          IHWMNY = 2
183:          ILALL = .FALSE.
184:          ILBACK = .FALSE.
185:       ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
186:          IHWMNY = 3
187:          ILALL = .TRUE.
188:          ILBACK = .TRUE.
189:       ELSE
190:          IHWMNY = -1
191:       END IF
192: *
193:       IF( LSAME( SIDE, 'R' ) ) THEN
194:          ISIDE = 1
195:          COMPL = .FALSE.
196:          COMPR = .TRUE.
197:       ELSE IF( LSAME( SIDE, 'L' ) ) THEN
198:          ISIDE = 2
199:          COMPL = .TRUE.
200:          COMPR = .FALSE.
201:       ELSE IF( LSAME( SIDE, 'B' ) ) THEN
202:          ISIDE = 3
203:          COMPL = .TRUE.
204:          COMPR = .TRUE.
205:       ELSE
206:          ISIDE = -1
207:       END IF
208: *
209:       INFO = 0
210:       IF( ISIDE.LT.0 ) THEN
211:          INFO = -1
212:       ELSE IF( IHWMNY.LT.0 ) THEN
213:          INFO = -2
214:       ELSE IF( N.LT.0 ) THEN
215:          INFO = -4
216:       ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
217:          INFO = -6
218:       ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
219:          INFO = -8
220:       END IF
221:       IF( INFO.NE.0 ) THEN
222:          CALL XERBLA( 'CTGEVC', -INFO )
223:          RETURN
224:       END IF
225: *
226: *     Count the number of eigenvectors
227: *
228:       IF( .NOT.ILALL ) THEN
229:          IM = 0
230:          DO 10 J = 1, N
231:             IF( SELECT( J ) )
232:      $         IM = IM + 1
233:    10    CONTINUE
234:       ELSE
235:          IM = N
236:       END IF
237: *
238: *     Check diagonal of B
239: *
240:       ILBBAD = .FALSE.
241:       DO 20 J = 1, N
242:          IF( AIMAG( P( J, J ) ).NE.ZERO )
243:      $      ILBBAD = .TRUE.
244:    20 CONTINUE
245: *
246:       IF( ILBBAD ) THEN
247:          INFO = -7
248:       ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
249:          INFO = -10
250:       ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
251:          INFO = -12
252:       ELSE IF( MM.LT.IM ) THEN
253:          INFO = -13
254:       END IF
255:       IF( INFO.NE.0 ) THEN
256:          CALL XERBLA( 'CTGEVC', -INFO )
257:          RETURN
258:       END IF
259: *
260: *     Quick return if possible
261: *
262:       M = IM
263:       IF( N.EQ.0 )
264:      $   RETURN
265: *
266: *     Machine Constants
267: *
268:       SAFMIN = SLAMCH( 'Safe minimum' )
269:       BIG = ONE / SAFMIN
270:       CALL SLABAD( SAFMIN, BIG )
271:       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
272:       SMALL = SAFMIN*N / ULP
273:       BIG = ONE / SMALL
274:       BIGNUM = ONE / ( SAFMIN*N )
275: *
276: *     Compute the 1-norm of each column of the strictly upper triangular
277: *     part of A and B to check for possible overflow in the triangular
278: *     solver.
279: *
280:       ANORM = ABS1( S( 1, 1 ) )
281:       BNORM = ABS1( P( 1, 1 ) )
282:       RWORK( 1 ) = ZERO
283:       RWORK( N+1 ) = ZERO
284:       DO 40 J = 2, N
285:          RWORK( J ) = ZERO
286:          RWORK( N+J ) = ZERO
287:          DO 30 I = 1, J - 1
288:             RWORK( J ) = RWORK( J ) + ABS1( S( I, J ) )
289:             RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) )
290:    30    CONTINUE
291:          ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) )
292:          BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) )
293:    40 CONTINUE
294: *
295:       ASCALE = ONE / MAX( ANORM, SAFMIN )
296:       BSCALE = ONE / MAX( BNORM, SAFMIN )
297: *
298: *     Left eigenvectors
299: *
300:       IF( COMPL ) THEN
301:          IEIG = 0
302: *
303: *        Main loop over eigenvalues
304: *
305:          DO 140 JE = 1, N
306:             IF( ILALL ) THEN
307:                ILCOMP = .TRUE.
308:             ELSE
309:                ILCOMP = SELECT( JE )
310:             END IF
311:             IF( ILCOMP ) THEN
312:                IEIG = IEIG + 1
313: *
314:                IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
315:      $             ABS( REAL( P( JE, JE ) ) ).LE.SAFMIN ) THEN
316: *
317: *                 Singular matrix pencil -- return unit eigenvector
318: *
319:                   DO 50 JR = 1, N
320:                      VL( JR, IEIG ) = CZERO
321:    50             CONTINUE
322:                   VL( IEIG, IEIG ) = CONE
323:                   GO TO 140
324:                END IF
325: *
326: *              Non-singular eigenvalue:
327: *              Compute coefficients  a  and  b  in
328: *                   H
329: *                 y  ( a A - b B ) = 0
330: *
331:                TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
332:      $                ABS( REAL( P( JE, JE ) ) )*BSCALE, SAFMIN )
333:                SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
334:                SBETA = ( TEMP*REAL( P( JE, JE ) ) )*BSCALE
335:                ACOEFF = SBETA*ASCALE
336:                BCOEFF = SALPHA*BSCALE
337: *
338: *              Scale to avoid underflow
339: *
340:                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
341:                LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
342:      $               SMALL
343: *
344:                SCALE = ONE
345:                IF( LSA )
346:      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
347:                IF( LSB )
348:      $            SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
349:      $                    MIN( BNORM, BIG ) )
350:                IF( LSA .OR. LSB ) THEN
351:                   SCALE = MIN( SCALE, ONE /
352:      $                    ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
353:      $                    ABS1( BCOEFF ) ) ) )
354:                   IF( LSA ) THEN
355:                      ACOEFF = ASCALE*( SCALE*SBETA )
356:                   ELSE
357:                      ACOEFF = SCALE*ACOEFF
358:                   END IF
359:                   IF( LSB ) THEN
360:                      BCOEFF = BSCALE*( SCALE*SALPHA )
361:                   ELSE
362:                      BCOEFF = SCALE*BCOEFF
363:                   END IF
364:                END IF
365: *
366:                ACOEFA = ABS( ACOEFF )
367:                BCOEFA = ABS1( BCOEFF )
368:                XMAX = ONE
369:                DO 60 JR = 1, N
370:                   WORK( JR ) = CZERO
371:    60          CONTINUE
372:                WORK( JE ) = CONE
373:                DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
374: *
375: *                                              H
376: *              Triangular solve of  (a A - b B)  y = 0
377: *
378: *                                      H
379: *              (rowwise in  (a A - b B) , or columnwise in a A - b B)
380: *
381:                DO 100 J = JE + 1, N
382: *
383: *                 Compute
384: *                       j-1
385: *                 SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
386: *                       k=je
387: *                 (Scale if necessary)
388: *
389:                   TEMP = ONE / XMAX
390:                   IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM*
391:      $                TEMP ) THEN
392:                      DO 70 JR = JE, J - 1
393:                         WORK( JR ) = TEMP*WORK( JR )
394:    70                CONTINUE
395:                      XMAX = ONE
396:                   END IF
397:                   SUMA = CZERO
398:                   SUMB = CZERO
399: *
400:                   DO 80 JR = JE, J - 1
401:                      SUMA = SUMA + CONJG( S( JR, J ) )*WORK( JR )
402:                      SUMB = SUMB + CONJG( P( JR, J ) )*WORK( JR )
403:    80             CONTINUE
404:                   SUM = ACOEFF*SUMA - CONJG( BCOEFF )*SUMB
405: *
406: *                 Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
407: *
408: *                 with scaling and perturbation of the denominator
409: *
410:                   D = CONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
411:                   IF( ABS1( D ).LE.DMIN )
412:      $               D = CMPLX( DMIN )
413: *
414:                   IF( ABS1( D ).LT.ONE ) THEN
415:                      IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN
416:                         TEMP = ONE / ABS1( SUM )
417:                         DO 90 JR = JE, J - 1
418:                            WORK( JR ) = TEMP*WORK( JR )
419:    90                   CONTINUE
420:                         XMAX = TEMP*XMAX
421:                         SUM = TEMP*SUM
422:                      END IF
423:                   END IF
424:                   WORK( J ) = CLADIV( -SUM, D )
425:                   XMAX = MAX( XMAX, ABS1( WORK( J ) ) )
426:   100          CONTINUE
427: *
428: *              Back transform eigenvector if HOWMNY='B'.
429: *
430:                IF( ILBACK ) THEN
431:                   CALL CGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL,
432:      $                        WORK( JE ), 1, CZERO, WORK( N+1 ), 1 )
433:                   ISRC = 2
434:                   IBEG = 1
435:                ELSE
436:                   ISRC = 1
437:                   IBEG = JE
438:                END IF
439: *
440: *              Copy and scale eigenvector into column of VL
441: *
442:                XMAX = ZERO
443:                DO 110 JR = IBEG, N
444:                   XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
445:   110          CONTINUE
446: *
447:                IF( XMAX.GT.SAFMIN ) THEN
448:                   TEMP = ONE / XMAX
449:                   DO 120 JR = IBEG, N
450:                      VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
451:   120             CONTINUE
452:                ELSE
453:                   IBEG = N + 1
454:                END IF
455: *
456:                DO 130 JR = 1, IBEG - 1
457:                   VL( JR, IEIG ) = CZERO
458:   130          CONTINUE
459: *
460:             END IF
461:   140    CONTINUE
462:       END IF
463: *
464: *     Right eigenvectors
465: *
466:       IF( COMPR ) THEN
467:          IEIG = IM + 1
468: *
469: *        Main loop over eigenvalues
470: *
471:          DO 250 JE = N, 1, -1
472:             IF( ILALL ) THEN
473:                ILCOMP = .TRUE.
474:             ELSE
475:                ILCOMP = SELECT( JE )
476:             END IF
477:             IF( ILCOMP ) THEN
478:                IEIG = IEIG - 1
479: *
480:                IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
481:      $             ABS( REAL( P( JE, JE ) ) ).LE.SAFMIN ) THEN
482: *
483: *                 Singular matrix pencil -- return unit eigenvector
484: *
485:                   DO 150 JR = 1, N
486:                      VR( JR, IEIG ) = CZERO
487:   150             CONTINUE
488:                   VR( IEIG, IEIG ) = CONE
489:                   GO TO 250
490:                END IF
491: *
492: *              Non-singular eigenvalue:
493: *              Compute coefficients  a  and  b  in
494: *
495: *              ( a A - b B ) x  = 0
496: *
497:                TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
498:      $                ABS( REAL( P( JE, JE ) ) )*BSCALE, SAFMIN )
499:                SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
500:                SBETA = ( TEMP*REAL( P( JE, JE ) ) )*BSCALE
501:                ACOEFF = SBETA*ASCALE
502:                BCOEFF = SALPHA*BSCALE
503: *
504: *              Scale to avoid underflow
505: *
506:                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
507:                LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
508:      $               SMALL
509: *
510:                SCALE = ONE
511:                IF( LSA )
512:      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
513:                IF( LSB )
514:      $            SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
515:      $                    MIN( BNORM, BIG ) )
516:                IF( LSA .OR. LSB ) THEN
517:                   SCALE = MIN( SCALE, ONE /
518:      $                    ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
519:      $                    ABS1( BCOEFF ) ) ) )
520:                   IF( LSA ) THEN
521:                      ACOEFF = ASCALE*( SCALE*SBETA )
522:                   ELSE
523:                      ACOEFF = SCALE*ACOEFF
524:                   END IF
525:                   IF( LSB ) THEN
526:                      BCOEFF = BSCALE*( SCALE*SALPHA )
527:                   ELSE
528:                      BCOEFF = SCALE*BCOEFF
529:                   END IF
530:                END IF
531: *
532:                ACOEFA = ABS( ACOEFF )
533:                BCOEFA = ABS1( BCOEFF )
534:                XMAX = ONE
535:                DO 160 JR = 1, N
536:                   WORK( JR ) = CZERO
537:   160          CONTINUE
538:                WORK( JE ) = CONE
539:                DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
540: *
541: *              Triangular solve of  (a A - b B) x = 0  (columnwise)
542: *
543: *              WORK(1:j-1) contains sums w,
544: *              WORK(j+1:JE) contains x
545: *
546:                DO 170 JR = 1, JE - 1
547:                   WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE )
548:   170          CONTINUE
549:                WORK( JE ) = CONE
550: *
551:                DO 210 J = JE - 1, 1, -1
552: *
553: *                 Form x(j) := - w(j) / d
554: *                 with scaling and perturbation of the denominator
555: *
556:                   D = ACOEFF*S( J, J ) - BCOEFF*P( J, J )
557:                   IF( ABS1( D ).LE.DMIN )
558:      $               D = CMPLX( DMIN )
559: *
560:                   IF( ABS1( D ).LT.ONE ) THEN
561:                      IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN
562:                         TEMP = ONE / ABS1( WORK( J ) )
563:                         DO 180 JR = 1, JE
564:                            WORK( JR ) = TEMP*WORK( JR )
565:   180                   CONTINUE
566:                      END IF
567:                   END IF
568: *
569:                   WORK( J ) = CLADIV( -WORK( J ), D )
570: *
571:                   IF( J.GT.1 ) THEN
572: *
573: *                    w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
574: *
575:                      IF( ABS1( WORK( J ) ).GT.ONE ) THEN
576:                         TEMP = ONE / ABS1( WORK( J ) )
577:                         IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE.
578:      $                      BIGNUM*TEMP ) THEN
579:                            DO 190 JR = 1, JE
580:                               WORK( JR ) = TEMP*WORK( JR )
581:   190                      CONTINUE
582:                         END IF
583:                      END IF
584: *
585:                      CA = ACOEFF*WORK( J )
586:                      CB = BCOEFF*WORK( J )
587:                      DO 200 JR = 1, J - 1
588:                         WORK( JR ) = WORK( JR ) + CA*S( JR, J ) -
589:      $                               CB*P( JR, J )
590:   200                CONTINUE
591:                   END IF
592:   210          CONTINUE
593: *
594: *              Back transform eigenvector if HOWMNY='B'.
595: *
596:                IF( ILBACK ) THEN
597:                   CALL CGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1,
598:      $                        CZERO, WORK( N+1 ), 1 )
599:                   ISRC = 2
600:                   IEND = N
601:                ELSE
602:                   ISRC = 1
603:                   IEND = JE
604:                END IF
605: *
606: *              Copy and scale eigenvector into column of VR
607: *
608:                XMAX = ZERO
609:                DO 220 JR = 1, IEND
610:                   XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
611:   220          CONTINUE
612: *
613:                IF( XMAX.GT.SAFMIN ) THEN
614:                   TEMP = ONE / XMAX
615:                   DO 230 JR = 1, IEND
616:                      VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
617:   230             CONTINUE
618:                ELSE
619:                   IEND = 0
620:                END IF
621: *
622:                DO 240 JR = IEND + 1, N
623:                   VR( JR, IEIG ) = CZERO
624:   240          CONTINUE
625: *
626:             END IF
627:   250    CONTINUE
628:       END IF
629: *
630:       RETURN
631: *
632: *     End of CTGEVC
633: *
634:       END
635: