1 SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
2 $ LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER HOWMNY, SIDE
11 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
12 * ..
13 * .. Array Arguments ..
14 LOGICAL SELECT( * )
15 DOUBLE PRECISION RWORK( * )
16 COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
17 $ VR( LDVR, * ), WORK( * )
18 * ..
19 *
20 *
21 * Purpose
22 * =======
23 *
24 * ZTGEVC computes some or all of the right and/or left eigenvectors of
25 * a pair of complex matrices (S,P), where S and P are upper triangular.
26 * Matrix pairs of this type are produced by the generalized Schur
27 * factorization of a complex matrix pair (A,B):
28 *
29 * A = Q*S*Z**H, B = Q*P*Z**H
30 *
31 * as computed by ZGGHRD + ZHGEQZ.
32 *
33 * The right eigenvector x and the left eigenvector y of (S,P)
34 * corresponding to an eigenvalue w are defined by:
35 *
36 * S*x = w*P*x, (y**H)*S = w*(y**H)*P,
37 *
38 * where y**H denotes the conjugate tranpose of y.
39 * The eigenvalues are not input to this routine, but are computed
40 * directly from the diagonal elements of S and P.
41 *
42 * This routine returns the matrices X and/or Y of right and left
43 * eigenvectors of (S,P), or the products Z*X and/or Q*Y,
44 * where Z and Q are input matrices.
45 * If Q and Z are the unitary factors from the generalized Schur
46 * factorization of a matrix pair (A,B), then Z*X and Q*Y
47 * are the matrices of right and left eigenvectors of (A,B).
48 *
49 * Arguments
50 * =========
51 *
52 * SIDE (input) CHARACTER*1
53 * = 'R': compute right eigenvectors only;
54 * = 'L': compute left eigenvectors only;
55 * = 'B': compute both right and left eigenvectors.
56 *
57 * HOWMNY (input) CHARACTER*1
58 * = 'A': compute all right and/or left eigenvectors;
59 * = 'B': compute all right and/or left eigenvectors,
60 * backtransformed by the matrices in VR and/or VL;
61 * = 'S': compute selected right and/or left eigenvectors,
62 * specified by the logical array SELECT.
63 *
64 * SELECT (input) LOGICAL array, dimension (N)
65 * If HOWMNY='S', SELECT specifies the eigenvectors to be
66 * computed. The eigenvector corresponding to the j-th
67 * eigenvalue is computed if SELECT(j) = .TRUE..
68 * Not referenced if HOWMNY = 'A' or 'B'.
69 *
70 * N (input) INTEGER
71 * The order of the matrices S and P. N >= 0.
72 *
73 * S (input) COMPLEX*16 array, dimension (LDS,N)
74 * The upper triangular matrix S from a generalized Schur
75 * factorization, as computed by ZHGEQZ.
76 *
77 * LDS (input) INTEGER
78 * The leading dimension of array S. LDS >= max(1,N).
79 *
80 * P (input) COMPLEX*16 array, dimension (LDP,N)
81 * The upper triangular matrix P from a generalized Schur
82 * factorization, as computed by ZHGEQZ. P must have real
83 * diagonal elements.
84 *
85 * LDP (input) INTEGER
86 * The leading dimension of array P. LDP >= max(1,N).
87 *
88 * VL (input/output) COMPLEX*16 array, dimension (LDVL,MM)
89 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
90 * contain an N-by-N matrix Q (usually the unitary matrix Q
91 * of left Schur vectors returned by ZHGEQZ).
92 * On exit, if SIDE = 'L' or 'B', VL contains:
93 * if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
94 * if HOWMNY = 'B', the matrix Q*Y;
95 * if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
96 * SELECT, stored consecutively in the columns of
97 * VL, in the same order as their eigenvalues.
98 * Not referenced if SIDE = 'R'.
99 *
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*16 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 ZHGEQZ).
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*16 array, dimension (2*N)
129 *
130 * RWORK (workspace) DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, ONE
140 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
141 COMPLEX*16 CZERO, CONE
142 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
143 $ CONE = ( 1.0D+0, 0.0D+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 DOUBLE PRECISION ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
151 $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
152 $ SCALE, SMALL, TEMP, ULP, XMAX
153 COMPLEX*16 BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
154 * ..
155 * .. External Functions ..
156 LOGICAL LSAME
157 DOUBLE PRECISION DLAMCH
158 COMPLEX*16 ZLADIV
159 EXTERNAL LSAME, DLAMCH, ZLADIV
160 * ..
161 * .. External Subroutines ..
162 EXTERNAL DLABAD, XERBLA, ZGEMV
163 * ..
164 * .. Intrinsic Functions ..
165 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
166 * ..
167 * .. Statement Functions ..
168 DOUBLE PRECISION ABS1
169 * ..
170 * .. Statement Function definitions ..
171 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( 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( 'ZTGEVC', -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( DIMAG( 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( 'ZTGEVC', -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 = DLAMCH( 'Safe minimum' )
269 BIG = ONE / SAFMIN
270 CALL DLABAD( SAFMIN, BIG )
271 ULP = DLAMCH( 'Epsilon' )*DLAMCH( '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( DBLE( 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( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
333 SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
334 SBETA = ( TEMP*DBLE( 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 + DCONJG( S( JR, J ) )*WORK( JR )
402 SUMB = SUMB + DCONJG( P( JR, J ) )*WORK( JR )
403 80 CONTINUE
404 SUM = ACOEFF*SUMA - DCONJG( 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 = DCONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
411 IF( ABS1( D ).LE.DMIN )
412 $ D = DCMPLX( 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 ) = ZLADIV( -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 ZGEMV( '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( DBLE( 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( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
499 SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
500 SBETA = ( TEMP*DBLE( 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 = DCMPLX( 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 ) = ZLADIV( -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 ZGEMV( '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 ZTGEVC
633 *
634 END
2 $ LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER HOWMNY, SIDE
11 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
12 * ..
13 * .. Array Arguments ..
14 LOGICAL SELECT( * )
15 DOUBLE PRECISION RWORK( * )
16 COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
17 $ VR( LDVR, * ), WORK( * )
18 * ..
19 *
20 *
21 * Purpose
22 * =======
23 *
24 * ZTGEVC computes some or all of the right and/or left eigenvectors of
25 * a pair of complex matrices (S,P), where S and P are upper triangular.
26 * Matrix pairs of this type are produced by the generalized Schur
27 * factorization of a complex matrix pair (A,B):
28 *
29 * A = Q*S*Z**H, B = Q*P*Z**H
30 *
31 * as computed by ZGGHRD + ZHGEQZ.
32 *
33 * The right eigenvector x and the left eigenvector y of (S,P)
34 * corresponding to an eigenvalue w are defined by:
35 *
36 * S*x = w*P*x, (y**H)*S = w*(y**H)*P,
37 *
38 * where y**H denotes the conjugate tranpose of y.
39 * The eigenvalues are not input to this routine, but are computed
40 * directly from the diagonal elements of S and P.
41 *
42 * This routine returns the matrices X and/or Y of right and left
43 * eigenvectors of (S,P), or the products Z*X and/or Q*Y,
44 * where Z and Q are input matrices.
45 * If Q and Z are the unitary factors from the generalized Schur
46 * factorization of a matrix pair (A,B), then Z*X and Q*Y
47 * are the matrices of right and left eigenvectors of (A,B).
48 *
49 * Arguments
50 * =========
51 *
52 * SIDE (input) CHARACTER*1
53 * = 'R': compute right eigenvectors only;
54 * = 'L': compute left eigenvectors only;
55 * = 'B': compute both right and left eigenvectors.
56 *
57 * HOWMNY (input) CHARACTER*1
58 * = 'A': compute all right and/or left eigenvectors;
59 * = 'B': compute all right and/or left eigenvectors,
60 * backtransformed by the matrices in VR and/or VL;
61 * = 'S': compute selected right and/or left eigenvectors,
62 * specified by the logical array SELECT.
63 *
64 * SELECT (input) LOGICAL array, dimension (N)
65 * If HOWMNY='S', SELECT specifies the eigenvectors to be
66 * computed. The eigenvector corresponding to the j-th
67 * eigenvalue is computed if SELECT(j) = .TRUE..
68 * Not referenced if HOWMNY = 'A' or 'B'.
69 *
70 * N (input) INTEGER
71 * The order of the matrices S and P. N >= 0.
72 *
73 * S (input) COMPLEX*16 array, dimension (LDS,N)
74 * The upper triangular matrix S from a generalized Schur
75 * factorization, as computed by ZHGEQZ.
76 *
77 * LDS (input) INTEGER
78 * The leading dimension of array S. LDS >= max(1,N).
79 *
80 * P (input) COMPLEX*16 array, dimension (LDP,N)
81 * The upper triangular matrix P from a generalized Schur
82 * factorization, as computed by ZHGEQZ. P must have real
83 * diagonal elements.
84 *
85 * LDP (input) INTEGER
86 * The leading dimension of array P. LDP >= max(1,N).
87 *
88 * VL (input/output) COMPLEX*16 array, dimension (LDVL,MM)
89 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
90 * contain an N-by-N matrix Q (usually the unitary matrix Q
91 * of left Schur vectors returned by ZHGEQZ).
92 * On exit, if SIDE = 'L' or 'B', VL contains:
93 * if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
94 * if HOWMNY = 'B', the matrix Q*Y;
95 * if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
96 * SELECT, stored consecutively in the columns of
97 * VL, in the same order as their eigenvalues.
98 * Not referenced if SIDE = 'R'.
99 *
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*16 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 ZHGEQZ).
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*16 array, dimension (2*N)
129 *
130 * RWORK (workspace) DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, ONE
140 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
141 COMPLEX*16 CZERO, CONE
142 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
143 $ CONE = ( 1.0D+0, 0.0D+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 DOUBLE PRECISION ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
151 $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
152 $ SCALE, SMALL, TEMP, ULP, XMAX
153 COMPLEX*16 BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
154 * ..
155 * .. External Functions ..
156 LOGICAL LSAME
157 DOUBLE PRECISION DLAMCH
158 COMPLEX*16 ZLADIV
159 EXTERNAL LSAME, DLAMCH, ZLADIV
160 * ..
161 * .. External Subroutines ..
162 EXTERNAL DLABAD, XERBLA, ZGEMV
163 * ..
164 * .. Intrinsic Functions ..
165 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
166 * ..
167 * .. Statement Functions ..
168 DOUBLE PRECISION ABS1
169 * ..
170 * .. Statement Function definitions ..
171 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( 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( 'ZTGEVC', -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( DIMAG( 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( 'ZTGEVC', -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 = DLAMCH( 'Safe minimum' )
269 BIG = ONE / SAFMIN
270 CALL DLABAD( SAFMIN, BIG )
271 ULP = DLAMCH( 'Epsilon' )*DLAMCH( '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( DBLE( 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( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
333 SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
334 SBETA = ( TEMP*DBLE( 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 + DCONJG( S( JR, J ) )*WORK( JR )
402 SUMB = SUMB + DCONJG( P( JR, J ) )*WORK( JR )
403 80 CONTINUE
404 SUM = ACOEFF*SUMA - DCONJG( 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 = DCONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
411 IF( ABS1( D ).LE.DMIN )
412 $ D = DCMPLX( 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 ) = ZLADIV( -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 ZGEMV( '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( DBLE( 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( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
499 SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
500 SBETA = ( TEMP*DBLE( 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 = DCMPLX( 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 ) = ZLADIV( -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 ZGEMV( '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 ZTGEVC
633 *
634 END