LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

◆ ztgevc()

 subroutine ztgevc ( character side, character howmny, logical, dimension( * ) select, integer n, complex*16, dimension( lds, * ) s, integer lds, complex*16, dimension( ldp, * ) p, integer ldp, complex*16, dimension( ldvl, * ) vl, integer ldvl, complex*16, dimension( ldvr, * ) vr, integer ldvr, integer mm, integer m, complex*16, dimension( * ) work, double precision, dimension( * ) rwork, integer info )

ZTGEVC

Purpose:
``` ZTGEVC computes some or all of the right and/or left eigenvectors of
a pair of complex matrices (S,P), where S and P are upper triangular.
Matrix pairs of this type are produced by the generalized Schur
factorization of a complex matrix pair (A,B):

A = Q*S*Z**H,  B = Q*P*Z**H

as computed by ZGGHRD + ZHGEQZ.

The right eigenvector x and the left eigenvector y of (S,P)
corresponding to an eigenvalue w are defined by:

S*x = w*P*x,  (y**H)*S = w*(y**H)*P,

where y**H denotes the conjugate transpose of y.
The eigenvalues are not input to this routine, but are computed
directly from the diagonal elements of S and P.

This routine returns the matrices X and/or Y of right and left
eigenvectors of (S,P), or the products Z*X and/or Q*Y,
where Z and Q are input matrices.
If Q and Z are the unitary factors from the generalized Schur
factorization of a matrix pair (A,B), then Z*X and Q*Y
are the matrices of right and left eigenvectors of (A,B).```
Parameters
 [in] SIDE ``` SIDE is CHARACTER*1 = 'R': compute right eigenvectors only; = 'L': compute left eigenvectors only; = 'B': compute both right and left eigenvectors.``` [in] HOWMNY ``` HOWMNY is CHARACTER*1 = 'A': compute all right and/or left eigenvectors; = 'B': compute all right and/or left eigenvectors, backtransformed by the matrices in VR and/or VL; = 'S': compute selected right and/or left eigenvectors, specified by the logical array SELECT.``` [in] SELECT ``` SELECT is LOGICAL array, dimension (N) If HOWMNY='S', SELECT specifies the eigenvectors to be computed. The eigenvector corresponding to the j-th eigenvalue is computed if SELECT(j) = .TRUE.. Not referenced if HOWMNY = 'A' or 'B'.``` [in] N ``` N is INTEGER The order of the matrices S and P. N >= 0.``` [in] S ``` S is COMPLEX*16 array, dimension (LDS,N) The upper triangular matrix S from a generalized Schur factorization, as computed by ZHGEQZ.``` [in] LDS ``` LDS is INTEGER The leading dimension of array S. LDS >= max(1,N).``` [in] P ``` P is COMPLEX*16 array, dimension (LDP,N) The upper triangular matrix P from a generalized Schur factorization, as computed by ZHGEQZ. P must have real diagonal elements.``` [in] LDP ``` LDP is INTEGER The leading dimension of array P. LDP >= max(1,N).``` [in,out] VL ``` VL is COMPLEX*16 array, dimension (LDVL,MM) On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must contain an N-by-N matrix Q (usually the unitary matrix Q of left Schur vectors returned by ZHGEQZ). On exit, if SIDE = 'L' or 'B', VL contains: if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P); if HOWMNY = 'B', the matrix Q*Y; if HOWMNY = 'S', the left eigenvectors of (S,P) specified by SELECT, stored consecutively in the columns of VL, in the same order as their eigenvalues. Not referenced if SIDE = 'R'.``` [in] LDVL ``` LDVL is INTEGER The leading dimension of array VL. LDVL >= 1, and if SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.``` [in,out] VR ``` VR is COMPLEX*16 array, dimension (LDVR,MM) On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must contain an N-by-N matrix Z (usually the unitary matrix Z of right Schur vectors returned by ZHGEQZ). On exit, if SIDE = 'R' or 'B', VR contains: if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P); if HOWMNY = 'B', the matrix Z*X; if HOWMNY = 'S', the right eigenvectors of (S,P) specified by SELECT, stored consecutively in the columns of VR, in the same order as their eigenvalues. Not referenced if SIDE = 'L'.``` [in] LDVR ``` LDVR is INTEGER The leading dimension of the array VR. LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.``` [in] MM ``` MM is INTEGER The number of columns in the arrays VL and/or VR. MM >= M.``` [out] M ``` M is INTEGER The number of columns in the arrays VL and/or VR actually used to store the eigenvectors. If HOWMNY = 'A' or 'B', M is set to N. Each selected eigenvector occupies one column.``` [out] WORK ` WORK is COMPLEX*16 array, dimension (2*N)` [out] RWORK ` RWORK is DOUBLE PRECISION array, dimension (2*N)` [out] INFO ``` INFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value.```

Definition at line 217 of file ztgevc.f.

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 DOUBLE PRECISION RWORK( * )
231 COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
232 \$ VR( LDVR, * ), WORK( * )
233* ..
234*
235*
236* =====================================================================
237*
238* .. Parameters ..
239 DOUBLE PRECISION ZERO, ONE
240 parameter( zero = 0.0d+0, one = 1.0d+0 )
241 COMPLEX*16 CZERO, CONE
242 parameter( czero = ( 0.0d+0, 0.0d+0 ),
243 \$ cone = ( 1.0d+0, 0.0d+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 DOUBLE PRECISION ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
251 \$ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
252 \$ SCALE, SMALL, TEMP, ULP, XMAX
253 COMPLEX*16 BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
254* ..
255* .. External Functions ..
256 LOGICAL LSAME
257 DOUBLE PRECISION DLAMCH
260* ..
261* .. External Subroutines ..
262 EXTERNAL xerbla, zgemv
263* ..
264* .. Intrinsic Functions ..
265 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min
266* ..
267* .. Statement Functions ..
268 DOUBLE PRECISION ABS1
269* ..
270* .. Statement Function definitions ..
271 abs1( x ) = abs( dble( x ) ) + abs( dimag( 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( 'ZTGEVC', -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*
341 DO 20 j = 1, n
342 IF( dimag( p( j, j ) ).NE.zero )
344 20 CONTINUE
345*
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( 'ZTGEVC', -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 = dlamch( 'Safe minimum' )
369 big = one / safmin
370 ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
371 small = safmin*n / ulp
372 big = one / small
373 bignum = one / ( safmin*n )
374*
375* Compute the 1-norm of each column of the strictly upper triangular
376* part of A and B to check for possible overflow in the triangular
377* solver.
378*
379 anorm = abs1( s( 1, 1 ) )
380 bnorm = abs1( p( 1, 1 ) )
381 rwork( 1 ) = zero
382 rwork( n+1 ) = zero
383 DO 40 j = 2, n
384 rwork( j ) = zero
385 rwork( n+j ) = zero
386 DO 30 i = 1, j - 1
387 rwork( j ) = rwork( j ) + abs1( s( i, j ) )
388 rwork( n+j ) = rwork( n+j ) + abs1( p( i, j ) )
389 30 CONTINUE
390 anorm = max( anorm, rwork( j )+abs1( s( j, j ) ) )
391 bnorm = max( bnorm, rwork( n+j )+abs1( p( j, j ) ) )
392 40 CONTINUE
393*
394 ascale = one / max( anorm, safmin )
395 bscale = one / max( bnorm, safmin )
396*
397* Left eigenvectors
398*
399 IF( compl ) THEN
400 ieig = 0
401*
402* Main loop over eigenvalues
403*
404 DO 140 je = 1, n
405 IF( ilall ) THEN
406 ilcomp = .true.
407 ELSE
408 ilcomp = SELECT( je )
409 END IF
410 IF( ilcomp ) THEN
411 ieig = ieig + 1
412*
413 IF( abs1( s( je, je ) ).LE.safmin .AND.
414 \$ abs( dble( p( je, je ) ) ).LE.safmin ) THEN
415*
416* Singular matrix pencil -- return unit eigenvector
417*
418 DO 50 jr = 1, n
419 vl( jr, ieig ) = czero
420 50 CONTINUE
421 vl( ieig, ieig ) = cone
422 GO TO 140
423 END IF
424*
425* Non-singular eigenvalue:
426* Compute coefficients a and b in
427* H
428* y ( a A - b B ) = 0
429*
430 temp = one / max( abs1( s( je, je ) )*ascale,
431 \$ abs( dble( p( je, je ) ) )*bscale, safmin )
432 salpha = ( temp*s( je, je ) )*ascale
433 sbeta = ( temp*dble( p( je, je ) ) )*bscale
434 acoeff = sbeta*ascale
435 bcoeff = salpha*bscale
436*
437* Scale to avoid underflow
438*
439 lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
440 lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
441 \$ small
442*
443 scale = one
444 IF( lsa )
445 \$ scale = ( small / abs( sbeta ) )*min( anorm, big )
446 IF( lsb )
447 \$ scale = max( scale, ( small / abs1( salpha ) )*
448 \$ min( bnorm, big ) )
449 IF( lsa .OR. lsb ) THEN
450 scale = min( scale, one /
451 \$ ( safmin*max( one, abs( acoeff ),
452 \$ abs1( bcoeff ) ) ) )
453 IF( lsa ) THEN
454 acoeff = ascale*( scale*sbeta )
455 ELSE
456 acoeff = scale*acoeff
457 END IF
458 IF( lsb ) THEN
459 bcoeff = bscale*( scale*salpha )
460 ELSE
461 bcoeff = scale*bcoeff
462 END IF
463 END IF
464*
465 acoefa = abs( acoeff )
466 bcoefa = abs1( bcoeff )
467 xmax = one
468 DO 60 jr = 1, n
469 work( jr ) = czero
470 60 CONTINUE
471 work( je ) = cone
472 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
473*
474* H
475* Triangular solve of (a A - b B) y = 0
476*
477* H
478* (rowwise in (a A - b B) , or columnwise in a A - b B)
479*
480 DO 100 j = je + 1, n
481*
482* Compute
483* j-1
484* SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
485* k=je
486* (Scale if necessary)
487*
488 temp = one / xmax
489 IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GT.bignum*
490 \$ temp ) THEN
491 DO 70 jr = je, j - 1
492 work( jr ) = temp*work( jr )
493 70 CONTINUE
494 xmax = one
495 END IF
496 suma = czero
497 sumb = czero
498*
499 DO 80 jr = je, j - 1
500 suma = suma + dconjg( s( jr, j ) )*work( jr )
501 sumb = sumb + dconjg( p( jr, j ) )*work( jr )
502 80 CONTINUE
503 sum = acoeff*suma - dconjg( bcoeff )*sumb
504*
505* Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
506*
507* with scaling and perturbation of the denominator
508*
509 d = dconjg( acoeff*s( j, j )-bcoeff*p( j, j ) )
510 IF( abs1( d ).LE.dmin )
511 \$ d = dcmplx( dmin )
512*
513 IF( abs1( d ).LT.one ) THEN
514 IF( abs1( sum ).GE.bignum*abs1( d ) ) THEN
515 temp = one / abs1( sum )
516 DO 90 jr = je, j - 1
517 work( jr ) = temp*work( jr )
518 90 CONTINUE
519 xmax = temp*xmax
520 sum = temp*sum
521 END IF
522 END IF
523 work( j ) = zladiv( -sum, d )
524 xmax = max( xmax, abs1( work( j ) ) )
525 100 CONTINUE
526*
527* Back transform eigenvector if HOWMNY='B'.
528*
529 IF( ilback ) THEN
530 CALL zgemv( 'N', n, n+1-je, cone, vl( 1, je ), ldvl,
531 \$ work( je ), 1, czero, work( n+1 ), 1 )
532 isrc = 2
533 ibeg = 1
534 ELSE
535 isrc = 1
536 ibeg = je
537 END IF
538*
539* Copy and scale eigenvector into column of VL
540*
541 xmax = zero
542 DO 110 jr = ibeg, n
543 xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
544 110 CONTINUE
545*
546 IF( xmax.GT.safmin ) THEN
547 temp = one / xmax
548 DO 120 jr = ibeg, n
549 vl( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
550 120 CONTINUE
551 ELSE
552 ibeg = n + 1
553 END IF
554*
555 DO 130 jr = 1, ibeg - 1
556 vl( jr, ieig ) = czero
557 130 CONTINUE
558*
559 END IF
560 140 CONTINUE
561 END IF
562*
563* Right eigenvectors
564*
565 IF( compr ) THEN
566 ieig = im + 1
567*
568* Main loop over eigenvalues
569*
570 DO 250 je = n, 1, -1
571 IF( ilall ) THEN
572 ilcomp = .true.
573 ELSE
574 ilcomp = SELECT( je )
575 END IF
576 IF( ilcomp ) THEN
577 ieig = ieig - 1
578*
579 IF( abs1( s( je, je ) ).LE.safmin .AND.
580 \$ abs( dble( p( je, je ) ) ).LE.safmin ) THEN
581*
582* Singular matrix pencil -- return unit eigenvector
583*
584 DO 150 jr = 1, n
585 vr( jr, ieig ) = czero
586 150 CONTINUE
587 vr( ieig, ieig ) = cone
588 GO TO 250
589 END IF
590*
591* Non-singular eigenvalue:
592* Compute coefficients a and b in
593*
594* ( a A - b B ) x = 0
595*
596 temp = one / max( abs1( s( je, je ) )*ascale,
597 \$ abs( dble( p( je, je ) ) )*bscale, safmin )
598 salpha = ( temp*s( je, je ) )*ascale
599 sbeta = ( temp*dble( p( je, je ) ) )*bscale
600 acoeff = sbeta*ascale
601 bcoeff = salpha*bscale
602*
603* Scale to avoid underflow
604*
605 lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
606 lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
607 \$ small
608*
609 scale = one
610 IF( lsa )
611 \$ scale = ( small / abs( sbeta ) )*min( anorm, big )
612 IF( lsb )
613 \$ scale = max( scale, ( small / abs1( salpha ) )*
614 \$ min( bnorm, big ) )
615 IF( lsa .OR. lsb ) THEN
616 scale = min( scale, one /
617 \$ ( safmin*max( one, abs( acoeff ),
618 \$ abs1( bcoeff ) ) ) )
619 IF( lsa ) THEN
620 acoeff = ascale*( scale*sbeta )
621 ELSE
622 acoeff = scale*acoeff
623 END IF
624 IF( lsb ) THEN
625 bcoeff = bscale*( scale*salpha )
626 ELSE
627 bcoeff = scale*bcoeff
628 END IF
629 END IF
630*
631 acoefa = abs( acoeff )
632 bcoefa = abs1( bcoeff )
633 xmax = one
634 DO 160 jr = 1, n
635 work( jr ) = czero
636 160 CONTINUE
637 work( je ) = cone
638 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
639*
640* Triangular solve of (a A - b B) x = 0 (columnwise)
641*
642* WORK(1:j-1) contains sums w,
643* WORK(j+1:JE) contains x
644*
645 DO 170 jr = 1, je - 1
646 work( jr ) = acoeff*s( jr, je ) - bcoeff*p( jr, je )
647 170 CONTINUE
648 work( je ) = cone
649*
650 DO 210 j = je - 1, 1, -1
651*
652* Form x(j) := - w(j) / d
653* with scaling and perturbation of the denominator
654*
655 d = acoeff*s( j, j ) - bcoeff*p( j, j )
656 IF( abs1( d ).LE.dmin )
657 \$ d = dcmplx( dmin )
658*
659 IF( abs1( d ).LT.one ) THEN
660 IF( abs1( work( j ) ).GE.bignum*abs1( d ) ) THEN
661 temp = one / abs1( work( j ) )
662 DO 180 jr = 1, je
663 work( jr ) = temp*work( jr )
664 180 CONTINUE
665 END IF
666 END IF
667*
668 work( j ) = zladiv( -work( j ), d )
669*
670 IF( j.GT.1 ) THEN
671*
672* w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
673*
674 IF( abs1( work( j ) ).GT.one ) THEN
675 temp = one / abs1( work( j ) )
676 IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GE.
677 \$ bignum*temp ) THEN
678 DO 190 jr = 1, je
679 work( jr ) = temp*work( jr )
680 190 CONTINUE
681 END IF
682 END IF
683*
684 ca = acoeff*work( j )
685 cb = bcoeff*work( j )
686 DO 200 jr = 1, j - 1
687 work( jr ) = work( jr ) + ca*s( jr, j ) -
688 \$ cb*p( jr, j )
689 200 CONTINUE
690 END IF
691 210 CONTINUE
692*
693* Back transform eigenvector if HOWMNY='B'.
694*
695 IF( ilback ) THEN
696 CALL zgemv( 'N', n, je, cone, vr, ldvr, work, 1,
697 \$ czero, work( n+1 ), 1 )
698 isrc = 2
699 iend = n
700 ELSE
701 isrc = 1
702 iend = je
703 END IF
704*
705* Copy and scale eigenvector into column of VR
706*
707 xmax = zero
708 DO 220 jr = 1, iend
709 xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
710 220 CONTINUE
711*
712 IF( xmax.GT.safmin ) THEN
713 temp = one / xmax
714 DO 230 jr = 1, iend
715 vr( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
716 230 CONTINUE
717 ELSE
718 iend = 0
719 END IF
720*
721 DO 240 jr = iend + 1, n
722 vr( jr, ieig ) = czero
723 240 CONTINUE
724*
725 END IF
726 250 CONTINUE
727 END IF
728*
729 RETURN
730*
731* End of ZTGEVC
732*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160