OpenRadioss 2025.1.11
OpenRadioss project
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*> \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
cmplx
float cmplx[2]
Definition pblas.h:136
slabad
subroutine slabad(small, large)
SLABAD
Definition slabad.f:74
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
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
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:158
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21