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