LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ 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 tranpose 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 Q (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.
Date
December 2016

Definition at line 221 of file ztgevc.f.

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