LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
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

Download ZTGEVC + dependencies [TGZ] [ZIP] [TXT]

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.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 215 of file ztgevc.f.

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 DOUBLE PRECISION RWORK( * )
229 COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
230 $ VR( LDVR, * ), WORK( * )
231* ..
232*
233*
234* =====================================================================
235*
236* .. Parameters ..
237 DOUBLE PRECISION ZERO, ONE
238 parameter( zero = 0.0d+0, one = 1.0d+0 )
239 COMPLEX*16 CZERO, CONE
240 parameter( czero = ( 0.0d+0, 0.0d+0 ),
241 $ cone = ( 1.0d+0, 0.0d+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 DOUBLE PRECISION ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
249 $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
250 $ SCALE, SMALL, TEMP, ULP, XMAX
251 COMPLEX*16 BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
252* ..
253* .. External Functions ..
254 LOGICAL LSAME
255 DOUBLE PRECISION DLAMCH
256 COMPLEX*16 ZLADIV
257 EXTERNAL lsame, dlamch, zladiv
258* ..
259* .. External Subroutines ..
260 EXTERNAL xerbla, zgemv
261* ..
262* .. Intrinsic Functions ..
263 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min
264* ..
265* .. Statement Functions ..
266 DOUBLE PRECISION ABS1
267* ..
268* .. Statement Function definitions ..
269 abs1( x ) = abs( dble( x ) ) + abs( dimag( 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( 'ZTGEVC', -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( dimag( 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( 'ZTGEVC', -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 = dlamch( 'Safe minimum' )
367 big = one / safmin
368 ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
369 small = safmin*n / ulp
370 big = one / small
371 bignum = one / ( safmin*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( dble( 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( dble( p( je, je ) ) )*bscale, safmin )
430 salpha = ( temp*s( je, je ) )*ascale
431 sbeta = ( temp*dble( 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 + dconjg( s( jr, j ) )*work( jr )
499 sumb = sumb + dconjg( p( jr, j ) )*work( jr )
500 80 CONTINUE
501 sum = acoeff*suma - dconjg( 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 = dconjg( acoeff*s( j, j )-bcoeff*p( j, j ) )
508 IF( abs1( d ).LE.dmin )
509 $ d = dcmplx( 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 ) = zladiv( -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 zgemv( '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( dble( 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( dble( p( je, je ) ) )*bscale, safmin )
597 salpha = ( temp*s( je, je ) )*ascale
598 sbeta = ( temp*dble( 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 = dcmplx( 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 ) = zladiv( -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 zgemv( '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 ZTGEVC
731*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
zladiv
complex *16 function zladiv(x, y)
ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition zladiv.f:62
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: