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