LAPACK  3.10.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 *> \ingroup complexGEcomputational
215 *
216 * =====================================================================
217  SUBROUTINE ctgevc( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
218  $ LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
219 *
220 * -- LAPACK computational routine --
221 * -- LAPACK is a software package provided by Univ. of Tennessee, --
222 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
223 *
224 * .. Scalar Arguments ..
225  CHARACTER HOWMNY, SIDE
226  INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
227 * ..
228 * .. Array Arguments ..
229  LOGICAL SELECT( * )
230  REAL RWORK( * )
231  COMPLEX P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
232  $ vr( ldvr, * ), work( * )
233 * ..
234 *
235 *
236 * =====================================================================
237 *
238 * .. Parameters ..
239  REAL ZERO, ONE
240  parameter( zero = 0.0e+0, one = 1.0e+0 )
241  COMPLEX CZERO, CONE
242  parameter( czero = ( 0.0e+0, 0.0e+0 ),
243  $ cone = ( 1.0e+0, 0.0e+0 ) )
244 * ..
245 * .. Local Scalars ..
246  LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
247  $ lsa, lsb
248  INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
249  $ j, je, jr
250  REAL ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
251  $ bignum, bnorm, bscale, dmin, safmin, sbeta,
252  $ scale, small, temp, ulp, xmax
253  COMPLEX BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
254 * ..
255 * .. External Functions ..
256  LOGICAL LSAME
257  REAL SLAMCH
258  COMPLEX CLADIV
259  EXTERNAL lsame, slamch, cladiv
260 * ..
261 * .. External Subroutines ..
262  EXTERNAL cgemv, slabad, xerbla
263 * ..
264 * .. Intrinsic Functions ..
265  INTRINSIC abs, aimag, cmplx, conjg, max, min, real
266 * ..
267 * .. Statement Functions ..
268  REAL ABS1
269 * ..
270 * .. Statement Function definitions ..
271  abs1( x ) = abs( real( x ) ) + abs( aimag( x ) )
272 * ..
273 * .. Executable Statements ..
274 *
275 * Decode and Test the input parameters
276 *
277  IF( lsame( howmny, 'A' ) ) THEN
278  ihwmny = 1
279  ilall = .true.
280  ilback = .false.
281  ELSE IF( lsame( howmny, 'S' ) ) THEN
282  ihwmny = 2
283  ilall = .false.
284  ilback = .false.
285  ELSE IF( lsame( howmny, 'B' ) ) THEN
286  ihwmny = 3
287  ilall = .true.
288  ilback = .true.
289  ELSE
290  ihwmny = -1
291  END IF
292 *
293  IF( lsame( side, 'R' ) ) THEN
294  iside = 1
295  compl = .false.
296  compr = .true.
297  ELSE IF( lsame( side, 'L' ) ) THEN
298  iside = 2
299  compl = .true.
300  compr = .false.
301  ELSE IF( lsame( side, 'B' ) ) THEN
302  iside = 3
303  compl = .true.
304  compr = .true.
305  ELSE
306  iside = -1
307  END IF
308 *
309  info = 0
310  IF( iside.LT.0 ) THEN
311  info = -1
312  ELSE IF( ihwmny.LT.0 ) THEN
313  info = -2
314  ELSE IF( n.LT.0 ) THEN
315  info = -4
316  ELSE IF( lds.LT.max( 1, n ) ) THEN
317  info = -6
318  ELSE IF( ldp.LT.max( 1, n ) ) THEN
319  info = -8
320  END IF
321  IF( info.NE.0 ) THEN
322  CALL xerbla( 'CTGEVC', -info )
323  RETURN
324  END IF
325 *
326 * Count the number of eigenvectors
327 *
328  IF( .NOT.ilall ) THEN
329  im = 0
330  DO 10 j = 1, n
331  IF( SELECT( j ) )
332  $ im = im + 1
333  10 CONTINUE
334  ELSE
335  im = n
336  END IF
337 *
338 * Check diagonal of B
339 *
340  ilbbad = .false.
341  DO 20 j = 1, n
342  IF( aimag( p( j, j ) ).NE.zero )
343  $ ilbbad = .true.
344  20 CONTINUE
345 *
346  IF( ilbbad ) THEN
347  info = -7
348  ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 ) THEN
349  info = -10
350  ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 ) THEN
351  info = -12
352  ELSE IF( mm.LT.im ) THEN
353  info = -13
354  END IF
355  IF( info.NE.0 ) THEN
356  CALL xerbla( 'CTGEVC', -info )
357  RETURN
358  END IF
359 *
360 * Quick return if possible
361 *
362  m = im
363  IF( n.EQ.0 )
364  $ RETURN
365 *
366 * Machine Constants
367 *
368  safmin = slamch( 'Safe minimum' )
369  big = one / safmin
370  CALL slabad( safmin, big )
371  ulp = slamch( 'Epsilon' )*slamch( 'Base' )
372  small = safmin*n / ulp
373  big = one / small
374  bignum = one / ( safmin*n )
375 *
376 * Compute the 1-norm of each column of the strictly upper triangular
377 * part of A and B to check for possible overflow in the triangular
378 * solver.
379 *
380  anorm = abs1( s( 1, 1 ) )
381  bnorm = abs1( p( 1, 1 ) )
382  rwork( 1 ) = zero
383  rwork( n+1 ) = zero
384  DO 40 j = 2, n
385  rwork( j ) = zero
386  rwork( n+j ) = zero
387  DO 30 i = 1, j - 1
388  rwork( j ) = rwork( j ) + abs1( s( i, j ) )
389  rwork( n+j ) = rwork( n+j ) + abs1( p( i, j ) )
390  30 CONTINUE
391  anorm = max( anorm, rwork( j )+abs1( s( j, j ) ) )
392  bnorm = max( bnorm, rwork( n+j )+abs1( p( j, j ) ) )
393  40 CONTINUE
394 *
395  ascale = one / max( anorm, safmin )
396  bscale = one / max( bnorm, safmin )
397 *
398 * Left eigenvectors
399 *
400  IF( compl ) THEN
401  ieig = 0
402 *
403 * Main loop over eigenvalues
404 *
405  DO 140 je = 1, n
406  IF( ilall ) THEN
407  ilcomp = .true.
408  ELSE
409  ilcomp = SELECT( je )
410  END IF
411  IF( ilcomp ) THEN
412  ieig = ieig + 1
413 *
414  IF( abs1( s( je, je ) ).LE.safmin .AND.
415  $ abs( real( p( je, je ) ) ).LE.safmin ) THEN
416 *
417 * Singular matrix pencil -- return unit eigenvector
418 *
419  DO 50 jr = 1, n
420  vl( jr, ieig ) = czero
421  50 CONTINUE
422  vl( ieig, ieig ) = cone
423  GO TO 140
424  END IF
425 *
426 * Non-singular eigenvalue:
427 * Compute coefficients a and b in
428 * H
429 * y ( a A - b B ) = 0
430 *
431  temp = one / max( abs1( s( je, je ) )*ascale,
432  $ abs( real( p( je, je ) ) )*bscale, safmin )
433  salpha = ( temp*s( je, je ) )*ascale
434  sbeta = ( temp*real( p( je, je ) ) )*bscale
435  acoeff = sbeta*ascale
436  bcoeff = salpha*bscale
437 *
438 * Scale to avoid underflow
439 *
440  lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
441  lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
442  $ small
443 *
444  scale = one
445  IF( lsa )
446  $ scale = ( small / abs( sbeta ) )*min( anorm, big )
447  IF( lsb )
448  $ scale = max( scale, ( small / abs1( salpha ) )*
449  $ min( bnorm, big ) )
450  IF( lsa .OR. lsb ) THEN
451  scale = min( scale, one /
452  $ ( safmin*max( one, abs( acoeff ),
453  $ abs1( bcoeff ) ) ) )
454  IF( lsa ) THEN
455  acoeff = ascale*( scale*sbeta )
456  ELSE
457  acoeff = scale*acoeff
458  END IF
459  IF( lsb ) THEN
460  bcoeff = bscale*( scale*salpha )
461  ELSE
462  bcoeff = scale*bcoeff
463  END IF
464  END IF
465 *
466  acoefa = abs( acoeff )
467  bcoefa = abs1( bcoeff )
468  xmax = one
469  DO 60 jr = 1, n
470  work( jr ) = czero
471  60 CONTINUE
472  work( je ) = cone
473  dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
474 *
475 * H
476 * Triangular solve of (a A - b B) y = 0
477 *
478 * H
479 * (rowwise in (a A - b B) , or columnwise in a A - b B)
480 *
481  DO 100 j = je + 1, n
482 *
483 * Compute
484 * j-1
485 * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
486 * k=je
487 * (Scale if necessary)
488 *
489  temp = one / xmax
490  IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GT.bignum*
491  $ temp ) THEN
492  DO 70 jr = je, j - 1
493  work( jr ) = temp*work( jr )
494  70 CONTINUE
495  xmax = one
496  END IF
497  suma = czero
498  sumb = czero
499 *
500  DO 80 jr = je, j - 1
501  suma = suma + conjg( s( jr, j ) )*work( jr )
502  sumb = sumb + conjg( p( jr, j ) )*work( jr )
503  80 CONTINUE
504  sum = acoeff*suma - conjg( bcoeff )*sumb
505 *
506 * Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
507 *
508 * with scaling and perturbation of the denominator
509 *
510  d = conjg( acoeff*s( j, j )-bcoeff*p( j, j ) )
511  IF( abs1( d ).LE.dmin )
512  $ d = cmplx( dmin )
513 *
514  IF( abs1( d ).LT.one ) THEN
515  IF( abs1( sum ).GE.bignum*abs1( d ) ) THEN
516  temp = one / abs1( sum )
517  DO 90 jr = je, j - 1
518  work( jr ) = temp*work( jr )
519  90 CONTINUE
520  xmax = temp*xmax
521  sum = temp*sum
522  END IF
523  END IF
524  work( j ) = cladiv( -sum, d )
525  xmax = max( xmax, abs1( work( j ) ) )
526  100 CONTINUE
527 *
528 * Back transform eigenvector if HOWMNY='B'.
529 *
530  IF( ilback ) THEN
531  CALL cgemv( 'N', n, n+1-je, cone, vl( 1, je ), ldvl,
532  $ work( je ), 1, czero, work( n+1 ), 1 )
533  isrc = 2
534  ibeg = 1
535  ELSE
536  isrc = 1
537  ibeg = je
538  END IF
539 *
540 * Copy and scale eigenvector into column of VL
541 *
542  xmax = zero
543  DO 110 jr = ibeg, n
544  xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
545  110 CONTINUE
546 *
547  IF( xmax.GT.safmin ) THEN
548  temp = one / xmax
549  DO 120 jr = ibeg, n
550  vl( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
551  120 CONTINUE
552  ELSE
553  ibeg = n + 1
554  END IF
555 *
556  DO 130 jr = 1, ibeg - 1
557  vl( jr, ieig ) = czero
558  130 CONTINUE
559 *
560  END IF
561  140 CONTINUE
562  END IF
563 *
564 * Right eigenvectors
565 *
566  IF( compr ) THEN
567  ieig = im + 1
568 *
569 * Main loop over eigenvalues
570 *
571  DO 250 je = n, 1, -1
572  IF( ilall ) THEN
573  ilcomp = .true.
574  ELSE
575  ilcomp = SELECT( je )
576  END IF
577  IF( ilcomp ) THEN
578  ieig = ieig - 1
579 *
580  IF( abs1( s( je, je ) ).LE.safmin .AND.
581  $ abs( real( p( je, je ) ) ).LE.safmin ) THEN
582 *
583 * Singular matrix pencil -- return unit eigenvector
584 *
585  DO 150 jr = 1, n
586  vr( jr, ieig ) = czero
587  150 CONTINUE
588  vr( ieig, ieig ) = cone
589  GO TO 250
590  END IF
591 *
592 * Non-singular eigenvalue:
593 * Compute coefficients a and b in
594 *
595 * ( a A - b B ) x = 0
596 *
597  temp = one / max( abs1( s( je, je ) )*ascale,
598  $ abs( real( p( je, je ) ) )*bscale, safmin )
599  salpha = ( temp*s( je, je ) )*ascale
600  sbeta = ( temp*real( p( je, je ) ) )*bscale
601  acoeff = sbeta*ascale
602  bcoeff = salpha*bscale
603 *
604 * Scale to avoid underflow
605 *
606  lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
607  lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
608  $ small
609 *
610  scale = one
611  IF( lsa )
612  $ scale = ( small / abs( sbeta ) )*min( anorm, big )
613  IF( lsb )
614  $ scale = max( scale, ( small / abs1( salpha ) )*
615  $ min( bnorm, big ) )
616  IF( lsa .OR. lsb ) THEN
617  scale = min( scale, one /
618  $ ( safmin*max( one, abs( acoeff ),
619  $ abs1( bcoeff ) ) ) )
620  IF( lsa ) THEN
621  acoeff = ascale*( scale*sbeta )
622  ELSE
623  acoeff = scale*acoeff
624  END IF
625  IF( lsb ) THEN
626  bcoeff = bscale*( scale*salpha )
627  ELSE
628  bcoeff = scale*bcoeff
629  END IF
630  END IF
631 *
632  acoefa = abs( acoeff )
633  bcoefa = abs1( bcoeff )
634  xmax = one
635  DO 160 jr = 1, n
636  work( jr ) = czero
637  160 CONTINUE
638  work( je ) = cone
639  dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
640 *
641 * Triangular solve of (a A - b B) x = 0 (columnwise)
642 *
643 * WORK(1:j-1) contains sums w,
644 * WORK(j+1:JE) contains x
645 *
646  DO 170 jr = 1, je - 1
647  work( jr ) = acoeff*s( jr, je ) - bcoeff*p( jr, je )
648  170 CONTINUE
649  work( je ) = cone
650 *
651  DO 210 j = je - 1, 1, -1
652 *
653 * Form x(j) := - w(j) / d
654 * with scaling and perturbation of the denominator
655 *
656  d = acoeff*s( j, j ) - bcoeff*p( j, j )
657  IF( abs1( d ).LE.dmin )
658  $ d = cmplx( dmin )
659 *
660  IF( abs1( d ).LT.one ) THEN
661  IF( abs1( work( j ) ).GE.bignum*abs1( d ) ) THEN
662  temp = one / abs1( work( j ) )
663  DO 180 jr = 1, je
664  work( jr ) = temp*work( jr )
665  180 CONTINUE
666  END IF
667  END IF
668 *
669  work( j ) = cladiv( -work( j ), d )
670 *
671  IF( j.GT.1 ) THEN
672 *
673 * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
674 *
675  IF( abs1( work( j ) ).GT.one ) THEN
676  temp = one / abs1( work( j ) )
677  IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GE.
678  $ bignum*temp ) THEN
679  DO 190 jr = 1, je
680  work( jr ) = temp*work( jr )
681  190 CONTINUE
682  END IF
683  END IF
684 *
685  ca = acoeff*work( j )
686  cb = bcoeff*work( j )
687  DO 200 jr = 1, j - 1
688  work( jr ) = work( jr ) + ca*s( jr, j ) -
689  $ cb*p( jr, j )
690  200 CONTINUE
691  END IF
692  210 CONTINUE
693 *
694 * Back transform eigenvector if HOWMNY='B'.
695 *
696  IF( ilback ) THEN
697  CALL cgemv( 'N', n, je, cone, vr, ldvr, work, 1,
698  $ czero, work( n+1 ), 1 )
699  isrc = 2
700  iend = n
701  ELSE
702  isrc = 1
703  iend = je
704  END IF
705 *
706 * Copy and scale eigenvector into column of VR
707 *
708  xmax = zero
709  DO 220 jr = 1, iend
710  xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
711  220 CONTINUE
712 *
713  IF( xmax.GT.safmin ) THEN
714  temp = one / xmax
715  DO 230 jr = 1, iend
716  vr( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
717  230 CONTINUE
718  ELSE
719  iend = 0
720  END IF
721 *
722  DO 240 jr = iend + 1, n
723  vr( jr, ieig ) = czero
724  240 CONTINUE
725 *
726  END IF
727  250 CONTINUE
728  END IF
729 *
730  RETURN
731 *
732 * End of CTGEVC
733 *
734  END
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
cgemv
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
ctgevc
subroutine ctgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
CTGEVC
Definition: ctgevc.f:219