LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
ctgevc.f
Go to the documentation of this file.
1 *> \brief \b CTGEVC
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CTGEVC + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgevc.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgevc.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgevc.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
22 * LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER HOWMNY, SIDE
26 * INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
27 * ..
28 * .. Array Arguments ..
29 * LOGICAL SELECT( * )
30 * REAL RWORK( * )
31 * COMPLEX P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
32 * $ VR( LDVR, * ), WORK( * )
33 * ..
34 *
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> CTGEVC computes some or all of the right and/or left eigenvectors of
43 *> a pair of complex matrices (S,P), where S and P are upper triangular.
44 *> Matrix pairs of this type are produced by the generalized Schur
45 *> factorization of a complex matrix pair (A,B):
46 *>
47 *> A = Q*S*Z**H, B = Q*P*Z**H
48 *>
49 *> as computed by CGGHRD + CHGEQZ.
50 *>
51 *> The right eigenvector x and the left eigenvector y of (S,P)
52 *> corresponding to an eigenvalue w are defined by:
53 *>
54 *> S*x = w*P*x, (y**H)*S = w*(y**H)*P,
55 *>
56 *> where y**H denotes the conjugate tranpose of y.
57 *> The eigenvalues are not input to this routine, but are computed
58 *> directly from the diagonal elements of S and P.
59 *>
60 *> This routine returns the matrices X and/or Y of right and left
61 *> eigenvectors of (S,P), or the products Z*X and/or Q*Y,
62 *> where Z and Q are input matrices.
63 *> If Q and Z are the unitary factors from the generalized Schur
64 *> factorization of a matrix pair (A,B), then Z*X and Q*Y
65 *> are the matrices of right and left eigenvectors of (A,B).
66 *> \endverbatim
67 *
68 * Arguments:
69 * ==========
70 *
71 *> \param[in] SIDE
72 *> \verbatim
73 *> SIDE is CHARACTER*1
74 *> = 'R': compute right eigenvectors only;
75 *> = 'L': compute left eigenvectors only;
76 *> = 'B': compute both right and left eigenvectors.
77 *> \endverbatim
78 *>
79 *> \param[in] HOWMNY
80 *> \verbatim
81 *> HOWMNY is CHARACTER*1
82 *> = 'A': compute all right and/or left eigenvectors;
83 *> = 'B': compute all right and/or left eigenvectors,
84 *> backtransformed by the matrices in VR and/or VL;
85 *> = 'S': compute selected right and/or left eigenvectors,
86 *> specified by the logical array SELECT.
87 *> \endverbatim
88 *>
89 *> \param[in] SELECT
90 *> \verbatim
91 *> SELECT is LOGICAL array, dimension (N)
92 *> If HOWMNY='S', SELECT specifies the eigenvectors to be
93 *> computed. The eigenvector corresponding to the j-th
94 *> eigenvalue is computed if SELECT(j) = .TRUE..
95 *> Not referenced if HOWMNY = 'A' or 'B'.
96 *> \endverbatim
97 *>
98 *> \param[in] N
99 *> \verbatim
100 *> N is INTEGER
101 *> The order of the matrices S and P. N >= 0.
102 *> \endverbatim
103 *>
104 *> \param[in] S
105 *> \verbatim
106 *> S is COMPLEX array, dimension (LDS,N)
107 *> The upper triangular matrix S from a generalized Schur
108 *> factorization, as computed by CHGEQZ.
109 *> \endverbatim
110 *>
111 *> \param[in] LDS
112 *> \verbatim
113 *> LDS is INTEGER
114 *> The leading dimension of array S. LDS >= max(1,N).
115 *> \endverbatim
116 *>
117 *> \param[in] P
118 *> \verbatim
119 *> P is COMPLEX array, dimension (LDP,N)
120 *> The upper triangular matrix P from a generalized Schur
121 *> factorization, as computed by CHGEQZ. P must have real
122 *> diagonal elements.
123 *> \endverbatim
124 *>
125 *> \param[in] LDP
126 *> \verbatim
127 *> LDP is INTEGER
128 *> The leading dimension of array P. LDP >= max(1,N).
129 *> \endverbatim
130 *>
131 *> \param[in,out] VL
132 *> \verbatim
133 *> VL is COMPLEX array, dimension (LDVL,MM)
134 *> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
135 *> contain an N-by-N matrix Q (usually the unitary matrix Q
136 *> of left Schur vectors returned by CHGEQZ).
137 *> On exit, if SIDE = 'L' or 'B', VL contains:
138 *> if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
139 *> if HOWMNY = 'B', the matrix Q*Y;
140 *> if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
141 *> SELECT, stored consecutively in the columns of
142 *> VL, in the same order as their eigenvalues.
143 *> Not referenced if SIDE = 'R'.
144 *> \endverbatim
145 *>
146 *> \param[in] LDVL
147 *> \verbatim
148 *> LDVL is INTEGER
149 *> The leading dimension of array VL. LDVL >= 1, and if
150 *> SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
151 *> \endverbatim
152 *>
153 *> \param[in,out] VR
154 *> \verbatim
155 *> VR is COMPLEX array, dimension (LDVR,MM)
156 *> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
157 *> contain an N-by-N matrix Q (usually the unitary matrix Z
158 *> of right Schur vectors returned by CHGEQZ).
159 *> On exit, if SIDE = 'R' or 'B', VR contains:
160 *> if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
161 *> if HOWMNY = 'B', the matrix Z*X;
162 *> if HOWMNY = 'S', the right eigenvectors of (S,P) specified by
163 *> SELECT, stored consecutively in the columns of
164 *> VR, in the same order as their eigenvalues.
165 *> Not referenced if SIDE = 'L'.
166 *> \endverbatim
167 *>
168 *> \param[in] LDVR
169 *> \verbatim
170 *> LDVR is INTEGER
171 *> The leading dimension of the array VR. LDVR >= 1, and if
172 *> SIDE = 'R' or 'B', LDVR >= N.
173 *> \endverbatim
174 *>
175 *> \param[in] MM
176 *> \verbatim
177 *> MM is INTEGER
178 *> The number of columns in the arrays VL and/or VR. MM >= M.
179 *> \endverbatim
180 *>
181 *> \param[out] M
182 *> \verbatim
183 *> M is INTEGER
184 *> The number of columns in the arrays VL and/or VR actually
185 *> used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
186 *> is set to N. Each selected eigenvector occupies one column.
187 *> \endverbatim
188 *>
189 *> \param[out] WORK
190 *> \verbatim
191 *> WORK is COMPLEX array, dimension (2*N)
192 *> \endverbatim
193 *>
194 *> \param[out] RWORK
195 *> \verbatim
196 *> RWORK is REAL array, dimension (2*N)
197 *> \endverbatim
198 *>
199 *> \param[out] INFO
200 *> \verbatim
201 *> INFO is INTEGER
202 *> = 0: successful exit.
203 *> < 0: if INFO = -i, the i-th argument had an illegal value.
204 *> \endverbatim
205 *
206 * Authors:
207 * ========
208 *
209 *> \author Univ. of Tennessee
210 *> \author Univ. of California Berkeley
211 *> \author Univ. of Colorado Denver
212 *> \author NAG Ltd.
213 *
214 *> \date November 2011
215 *
216 *> \ingroup complexGEcomputational
217 *
218 * =====================================================================
219  SUBROUTINE ctgevc( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
220  $ ldvl, vr, ldvr, mm, m, work, rwork, info )
221 *
222 * -- LAPACK computational routine (version 3.4.0) --
223 * -- LAPACK is a software package provided by Univ. of Tennessee, --
224 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
225 * November 2011
226 *
227 * .. Scalar Arguments ..
228  CHARACTER HOWMNY, SIDE
229  INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
230 * ..
231 * .. Array Arguments ..
232  LOGICAL SELECT( * )
233  REAL RWORK( * )
234  COMPLEX P( ldp, * ), S( lds, * ), VL( ldvl, * ),
235  $ vr( ldvr, * ), work( * )
236 * ..
237 *
238 *
239 * =====================================================================
240 *
241 * .. Parameters ..
242  REAL ZERO, ONE
243  parameter ( zero = 0.0e+0, one = 1.0e+0 )
244  COMPLEX CZERO, CONE
245  parameter ( czero = ( 0.0e+0, 0.0e+0 ),
246  $ cone = ( 1.0e+0, 0.0e+0 ) )
247 * ..
248 * .. Local Scalars ..
249  LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
250  $ lsa, lsb
251  INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
252  $ j, je, jr
253  REAL ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
254  $ bignum, bnorm, bscale, dmin, safmin, sbeta,
255  $ scale, small, temp, ulp, xmax
256  COMPLEX BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
257 * ..
258 * .. External Functions ..
259  LOGICAL LSAME
260  REAL SLAMCH
261  COMPLEX CLADIV
262  EXTERNAL lsame, slamch, cladiv
263 * ..
264 * .. External Subroutines ..
265  EXTERNAL cgemv, slabad, xerbla
266 * ..
267 * .. Intrinsic Functions ..
268  INTRINSIC abs, aimag, cmplx, conjg, max, min, real
269 * ..
270 * .. Statement Functions ..
271  REAL ABS1
272 * ..
273 * .. Statement Function definitions ..
274  abs1( x ) = abs( REAL( X ) ) + abs( AIMAG( x ) )
275 * ..
276 * .. Executable Statements ..
277 *
278 * Decode and Test the input parameters
279 *
280  IF( lsame( howmny, 'A' ) ) THEN
281  ihwmny = 1
282  ilall = .true.
283  ilback = .false.
284  ELSE IF( lsame( howmny, 'S' ) ) THEN
285  ihwmny = 2
286  ilall = .false.
287  ilback = .false.
288  ELSE IF( lsame( howmny, 'B' ) ) THEN
289  ihwmny = 3
290  ilall = .true.
291  ilback = .true.
292  ELSE
293  ihwmny = -1
294  END IF
295 *
296  IF( lsame( side, 'R' ) ) THEN
297  iside = 1
298  compl = .false.
299  compr = .true.
300  ELSE IF( lsame( side, 'L' ) ) THEN
301  iside = 2
302  compl = .true.
303  compr = .false.
304  ELSE IF( lsame( side, 'B' ) ) THEN
305  iside = 3
306  compl = .true.
307  compr = .true.
308  ELSE
309  iside = -1
310  END IF
311 *
312  info = 0
313  IF( iside.LT.0 ) THEN
314  info = -1
315  ELSE IF( ihwmny.LT.0 ) THEN
316  info = -2
317  ELSE IF( n.LT.0 ) THEN
318  info = -4
319  ELSE IF( lds.LT.max( 1, n ) ) THEN
320  info = -6
321  ELSE IF( ldp.LT.max( 1, n ) ) THEN
322  info = -8
323  END IF
324  IF( info.NE.0 ) THEN
325  CALL xerbla( 'CTGEVC', -info )
326  RETURN
327  END IF
328 *
329 * Count the number of eigenvectors
330 *
331  IF( .NOT.ilall ) THEN
332  im = 0
333  DO 10 j = 1, n
334  IF( SELECT( j ) )
335  $ im = im + 1
336  10 CONTINUE
337  ELSE
338  im = n
339  END IF
340 *
341 * Check diagonal of B
342 *
343  ilbbad = .false.
344  DO 20 j = 1, n
345  IF( aimag( p( j, j ) ).NE.zero )
346  $ ilbbad = .true.
347  20 CONTINUE
348 *
349  IF( ilbbad ) THEN
350  info = -7
351  ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 ) THEN
352  info = -10
353  ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 ) THEN
354  info = -12
355  ELSE IF( mm.LT.im ) THEN
356  info = -13
357  END IF
358  IF( info.NE.0 ) THEN
359  CALL xerbla( 'CTGEVC', -info )
360  RETURN
361  END IF
362 *
363 * Quick return if possible
364 *
365  m = im
366  IF( n.EQ.0 )
367  $ RETURN
368 *
369 * Machine Constants
370 *
371  safmin = slamch( 'Safe minimum' )
372  big = one / safmin
373  CALL slabad( safmin, big )
374  ulp = slamch( 'Epsilon' )*slamch( 'Base' )
375  small = safmin*n / ulp
376  big = one / small
377  bignum = one / ( safmin*n )
378 *
379 * Compute the 1-norm of each column of the strictly upper triangular
380 * part of A and B to check for possible overflow in the triangular
381 * solver.
382 *
383  anorm = abs1( s( 1, 1 ) )
384  bnorm = abs1( p( 1, 1 ) )
385  rwork( 1 ) = zero
386  rwork( n+1 ) = zero
387  DO 40 j = 2, n
388  rwork( j ) = zero
389  rwork( n+j ) = zero
390  DO 30 i = 1, j - 1
391  rwork( j ) = rwork( j ) + abs1( s( i, j ) )
392  rwork( n+j ) = rwork( n+j ) + abs1( p( i, j ) )
393  30 CONTINUE
394  anorm = max( anorm, rwork( j )+abs1( s( j, j ) ) )
395  bnorm = max( bnorm, rwork( n+j )+abs1( p( j, j ) ) )
396  40 CONTINUE
397 *
398  ascale = one / max( anorm, safmin )
399  bscale = one / max( bnorm, safmin )
400 *
401 * Left eigenvectors
402 *
403  IF( compl ) THEN
404  ieig = 0
405 *
406 * Main loop over eigenvalues
407 *
408  DO 140 je = 1, n
409  IF( ilall ) THEN
410  ilcomp = .true.
411  ELSE
412  ilcomp = SELECT( je )
413  END IF
414  IF( ilcomp ) THEN
415  ieig = ieig + 1
416 *
417  IF( abs1( s( je, je ) ).LE.safmin .AND.
418  $ abs( REAL( P( JE, JE ) ) ).LE.safmin ) then
419 *
420 * Singular matrix pencil -- return unit eigenvector
421 *
422  DO 50 jr = 1, n
423  vl( jr, ieig ) = czero
424  50 CONTINUE
425  vl( ieig, ieig ) = cone
426  GO TO 140
427  END IF
428 *
429 * Non-singular eigenvalue:
430 * Compute coefficients a and b in
431 * H
432 * y ( a A - b B ) = 0
433 *
434  temp = one / max( abs1( s( je, je ) )*ascale,
435  $ abs( REAL( P( JE, JE ) ) )*bscale, safmin )
436  salpha = ( temp*s( je, je ) )*ascale
437  sbeta = ( temp*REAL( P( JE, JE ) ) )*bscale
438  acoeff = sbeta*ascale
439  bcoeff = salpha*bscale
440 *
441 * Scale to avoid underflow
442 *
443  lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
444  lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
445  $ small
446 *
447  scale = one
448  IF( lsa )
449  $ scale = ( small / abs( sbeta ) )*min( anorm, big )
450  IF( lsb )
451  $ scale = max( scale, ( small / abs1( salpha ) )*
452  $ min( bnorm, big ) )
453  IF( lsa .OR. lsb ) THEN
454  scale = min( scale, one /
455  $ ( safmin*max( one, abs( acoeff ),
456  $ abs1( bcoeff ) ) ) )
457  IF( lsa ) THEN
458  acoeff = ascale*( scale*sbeta )
459  ELSE
460  acoeff = scale*acoeff
461  END IF
462  IF( lsb ) THEN
463  bcoeff = bscale*( scale*salpha )
464  ELSE
465  bcoeff = scale*bcoeff
466  END IF
467  END IF
468 *
469  acoefa = abs( acoeff )
470  bcoefa = abs1( bcoeff )
471  xmax = one
472  DO 60 jr = 1, n
473  work( jr ) = czero
474  60 CONTINUE
475  work( je ) = cone
476  dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
477 *
478 * H
479 * Triangular solve of (a A - b B) y = 0
480 *
481 * H
482 * (rowwise in (a A - b B) , or columnwise in a A - b B)
483 *
484  DO 100 j = je + 1, n
485 *
486 * Compute
487 * j-1
488 * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
489 * k=je
490 * (Scale if necessary)
491 *
492  temp = one / xmax
493  IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GT.bignum*
494  $ temp ) THEN
495  DO 70 jr = je, j - 1
496  work( jr ) = temp*work( jr )
497  70 CONTINUE
498  xmax = one
499  END IF
500  suma = czero
501  sumb = czero
502 *
503  DO 80 jr = je, j - 1
504  suma = suma + conjg( s( jr, j ) )*work( jr )
505  sumb = sumb + conjg( p( jr, j ) )*work( jr )
506  80 CONTINUE
507  sum = acoeff*suma - conjg( bcoeff )*sumb
508 *
509 * Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
510 *
511 * with scaling and perturbation of the denominator
512 *
513  d = conjg( acoeff*s( j, j )-bcoeff*p( j, j ) )
514  IF( abs1( d ).LE.dmin )
515  $ d = cmplx( dmin )
516 *
517  IF( abs1( d ).LT.one ) THEN
518  IF( abs1( sum ).GE.bignum*abs1( d ) ) THEN
519  temp = one / abs1( sum )
520  DO 90 jr = je, j - 1
521  work( jr ) = temp*work( jr )
522  90 CONTINUE
523  xmax = temp*xmax
524  sum = temp*sum
525  END IF
526  END IF
527  work( j ) = cladiv( -sum, d )
528  xmax = max( xmax, abs1( work( j ) ) )
529  100 CONTINUE
530 *
531 * Back transform eigenvector if HOWMNY='B'.
532 *
533  IF( ilback ) THEN
534  CALL cgemv( 'N', n, n+1-je, cone, vl( 1, je ), ldvl,
535  $ work( je ), 1, czero, work( n+1 ), 1 )
536  isrc = 2
537  ibeg = 1
538  ELSE
539  isrc = 1
540  ibeg = je
541  END IF
542 *
543 * Copy and scale eigenvector into column of VL
544 *
545  xmax = zero
546  DO 110 jr = ibeg, n
547  xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
548  110 CONTINUE
549 *
550  IF( xmax.GT.safmin ) THEN
551  temp = one / xmax
552  DO 120 jr = ibeg, n
553  vl( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
554  120 CONTINUE
555  ELSE
556  ibeg = n + 1
557  END IF
558 *
559  DO 130 jr = 1, ibeg - 1
560  vl( jr, ieig ) = czero
561  130 CONTINUE
562 *
563  END IF
564  140 CONTINUE
565  END IF
566 *
567 * Right eigenvectors
568 *
569  IF( compr ) THEN
570  ieig = im + 1
571 *
572 * Main loop over eigenvalues
573 *
574  DO 250 je = n, 1, -1
575  IF( ilall ) THEN
576  ilcomp = .true.
577  ELSE
578  ilcomp = SELECT( je )
579  END IF
580  IF( ilcomp ) THEN
581  ieig = ieig - 1
582 *
583  IF( abs1( s( je, je ) ).LE.safmin .AND.
584  $ abs( REAL( P( JE, JE ) ) ).LE.safmin ) then
585 *
586 * Singular matrix pencil -- return unit eigenvector
587 *
588  DO 150 jr = 1, n
589  vr( jr, ieig ) = czero
590  150 CONTINUE
591  vr( ieig, ieig ) = cone
592  GO TO 250
593  END IF
594 *
595 * Non-singular eigenvalue:
596 * Compute coefficients a and b in
597 *
598 * ( a A - b B ) x = 0
599 *
600  temp = one / max( abs1( s( je, je ) )*ascale,
601  $ abs( REAL( P( JE, JE ) ) )*bscale, safmin )
602  salpha = ( temp*s( je, je ) )*ascale
603  sbeta = ( temp*REAL( P( JE, JE ) ) )*bscale
604  acoeff = sbeta*ascale
605  bcoeff = salpha*bscale
606 *
607 * Scale to avoid underflow
608 *
609  lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
610  lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
611  $ small
612 *
613  scale = one
614  IF( lsa )
615  $ scale = ( small / abs( sbeta ) )*min( anorm, big )
616  IF( lsb )
617  $ scale = max( scale, ( small / abs1( salpha ) )*
618  $ min( bnorm, big ) )
619  IF( lsa .OR. lsb ) THEN
620  scale = min( scale, one /
621  $ ( safmin*max( one, abs( acoeff ),
622  $ abs1( bcoeff ) ) ) )
623  IF( lsa ) THEN
624  acoeff = ascale*( scale*sbeta )
625  ELSE
626  acoeff = scale*acoeff
627  END IF
628  IF( lsb ) THEN
629  bcoeff = bscale*( scale*salpha )
630  ELSE
631  bcoeff = scale*bcoeff
632  END IF
633  END IF
634 *
635  acoefa = abs( acoeff )
636  bcoefa = abs1( bcoeff )
637  xmax = one
638  DO 160 jr = 1, n
639  work( jr ) = czero
640  160 CONTINUE
641  work( je ) = cone
642  dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
643 *
644 * Triangular solve of (a A - b B) x = 0 (columnwise)
645 *
646 * WORK(1:j-1) contains sums w,
647 * WORK(j+1:JE) contains x
648 *
649  DO 170 jr = 1, je - 1
650  work( jr ) = acoeff*s( jr, je ) - bcoeff*p( jr, je )
651  170 CONTINUE
652  work( je ) = cone
653 *
654  DO 210 j = je - 1, 1, -1
655 *
656 * Form x(j) := - w(j) / d
657 * with scaling and perturbation of the denominator
658 *
659  d = acoeff*s( j, j ) - bcoeff*p( j, j )
660  IF( abs1( d ).LE.dmin )
661  $ d = cmplx( dmin )
662 *
663  IF( abs1( d ).LT.one ) THEN
664  IF( abs1( work( j ) ).GE.bignum*abs1( d ) ) THEN
665  temp = one / abs1( work( j ) )
666  DO 180 jr = 1, je
667  work( jr ) = temp*work( jr )
668  180 CONTINUE
669  END IF
670  END IF
671 *
672  work( j ) = cladiv( -work( j ), d )
673 *
674  IF( j.GT.1 ) THEN
675 *
676 * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
677 *
678  IF( abs1( work( j ) ).GT.one ) THEN
679  temp = one / abs1( work( j ) )
680  IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GE.
681  $ bignum*temp ) THEN
682  DO 190 jr = 1, je
683  work( jr ) = temp*work( jr )
684  190 CONTINUE
685  END IF
686  END IF
687 *
688  ca = acoeff*work( j )
689  cb = bcoeff*work( j )
690  DO 200 jr = 1, j - 1
691  work( jr ) = work( jr ) + ca*s( jr, j ) -
692  $ cb*p( jr, j )
693  200 CONTINUE
694  END IF
695  210 CONTINUE
696 *
697 * Back transform eigenvector if HOWMNY='B'.
698 *
699  IF( ilback ) THEN
700  CALL cgemv( 'N', n, je, cone, vr, ldvr, work, 1,
701  $ czero, work( n+1 ), 1 )
702  isrc = 2
703  iend = n
704  ELSE
705  isrc = 1
706  iend = je
707  END IF
708 *
709 * Copy and scale eigenvector into column of VR
710 *
711  xmax = zero
712  DO 220 jr = 1, iend
713  xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
714  220 CONTINUE
715 *
716  IF( xmax.GT.safmin ) THEN
717  temp = one / xmax
718  DO 230 jr = 1, iend
719  vr( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
720  230 CONTINUE
721  ELSE
722  iend = 0
723  END IF
724 *
725  DO 240 jr = iend + 1, n
726  vr( jr, ieig ) = czero
727  240 CONTINUE
728 *
729  END IF
730  250 CONTINUE
731  END IF
732 *
733  RETURN
734 *
735 * End of CTGEVC
736 *
737  END
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cgemv
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
ctgevc
subroutine ctgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
CTGEVC
Definition: ctgevc.f:221