LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dtgevc()

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

DTGEVC

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

Purpose:
!>
!> DTGEVC computes some or all of the right and/or left eigenvectors of
!> a pair of real matrices (S,P), where S is a quasi-triangular matrix
!> and P is upper triangular.  Matrix pairs of this type are produced by
!> the generalized Schur factorization of a matrix pair (A,B):
!>
!>    A = Q*S*Z**T,  B = Q*P*Z**T
!>
!> as computed by DGGHRD + DHGEQZ.
!>
!> 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 blocks 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 orthogonal 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.  If w(j) is a real eigenvalue, the corresponding
!>          real eigenvector is computed if SELECT(j) is .TRUE..
!>          If w(j) and w(j+1) are the real and imaginary parts of a
!>          complex eigenvalue, the corresponding complex eigenvector
!>          is computed if either SELECT(j) or SELECT(j+1) is .TRUE.,
!>          and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is
!>          set to .FALSE..
!>          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 DOUBLE PRECISION array, dimension (LDS,N)
!>          The upper quasi-triangular matrix S from a generalized Schur
!>          factorization, as computed by DHGEQZ.
!>
[in]LDS
!>          LDS is INTEGER
!>          The leading dimension of array S.  LDS >= max(1,N).
!>
[in]P
!>          P is DOUBLE PRECISION array, dimension (LDP,N)
!>          The upper triangular matrix P from a generalized Schur
!>          factorization, as computed by DHGEQZ.
!>          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
!>          of S must be in positive diagonal form.
!>
[in]LDP
!>          LDP is INTEGER
!>          The leading dimension of array P.  LDP >= max(1,N).
!>
[in,out]VL
!>          VL is DOUBLE PRECISION 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 orthogonal matrix Q
!>          of left Schur vectors returned by DHGEQZ).
!>          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.
!>
!>          A complex eigenvector corresponding to a complex eigenvalue
!>          is stored in two consecutive columns, the first holding the
!>          real part, and the second the imaginary part.
!>
!>          Not referenced if SIDE = 'R'.
!>
[in]LDVL
!>          LDVL is INTEGER
!>          The leading dimension of array VL.  LDVL >= 1, and if
!>          SIDE = 'L' or 'B', LDVL >= N.
!>
[in,out]VR
!>          VR is DOUBLE PRECISION 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 orthogonal matrix Z
!>          of right Schur vectors returned by DHGEQZ).
!>
!>          On exit, if SIDE = 'R' or 'B', VR contains:
!>          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
!>          if HOWMNY = 'B' or 'b', the matrix Z*X;
!>          if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
!>                      specified by SELECT, stored consecutively in the
!>                      columns of VR, in the same order as their
!>                      eigenvalues.
!>
!>          A complex eigenvector corresponding to a complex eigenvalue
!>          is stored in two consecutive columns, the first holding the
!>          real part and the second the imaginary part.
!>
!>          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 real eigenvector occupies one
!>          column and each selected complex eigenvector occupies two
!>          columns.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (6*N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  the 2-by-2 block (INFO:INFO+1) does not have a complex
!>                eigenvalue.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  Allocation of workspace:
!>  ---------- -- ---------
!>
!>     WORK( j ) = 1-norm of j-th column of A, above the diagonal
!>     WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
!>     WORK( 2*N+1:3*N ) = real part of eigenvector
!>     WORK( 3*N+1:4*N ) = imaginary part of eigenvector
!>     WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
!>     WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
!>
!>  Rowwise vs. columnwise solution methods:
!>  ------- --  ---------- -------- -------
!>
!>  Finding a generalized eigenvector consists basically of solving the
!>  singular triangular system
!>
!>   (A - w B) x = 0     (for right) or:   (A - w B)**H y = 0  (for left)
!>
!>  Consider finding the i-th right eigenvector (assume all eigenvalues
!>  are real). The equation to be solved is:
!>       n                   i
!>  0 = sum  C(j,k) v(k)  = sum  C(j,k) v(k)     for j = i,. . .,1
!>      k=j                 k=j
!>
!>  where  C = (A - w B)  (The components v(i+1:n) are 0.)
!>
!>  The  method is:
!>
!>  (1)  v(i) := 1
!>  for j = i-1,. . .,1:
!>                          i
!>      (2) compute  s = - sum C(j,k) v(k)   and
!>                        k=j+1
!>
!>      (3) v(j) := s / C(j,j)
!>
!>  Step 2 is sometimes called the  step, since it is an
!>  inner product between the j-th row and the portion of the eigenvector
!>  that has been computed so far.
!>
!>  The  method consists basically in doing the sums
!>  for all the rows in parallel.  As each v(j) is computed, the
!>  contribution of v(j) times the j-th column of C is added to the
!>  partial sums.  Since FORTRAN arrays are stored columnwise, this has
!>  the advantage that at each step, the elements of C that are accessed
!>  are adjacent to one another, whereas with the rowwise method, the
!>  elements accessed at a step are spaced LDS (and LDP) words apart.
!>
!>  When finding left eigenvectors, the matrix in question is the
!>  transpose of the one in storage, so the rowwise method then
!>  actually accesses columns of A and B at each step, and so is the
!>  preferred method.
!>

Definition at line 291 of file dtgevc.f.

293*
294* -- LAPACK computational routine --
295* -- LAPACK is a software package provided by Univ. of Tennessee, --
296* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
297*
298* .. Scalar Arguments ..
299 CHARACTER HOWMNY, SIDE
300 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
301* ..
302* .. Array Arguments ..
303 LOGICAL SELECT( * )
304 DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
305 $ VR( LDVR, * ), WORK( * )
306* ..
307*
308*
309* =====================================================================
310*
311* .. Parameters ..
312 DOUBLE PRECISION ZERO, ONE, SAFETY
313 parameter( zero = 0.0d+0, one = 1.0d+0,
314 $ safety = 1.0d+2 )
315* ..
316* .. Local Scalars ..
317 LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
318 $ ILBBAD, ILCOMP, ILCPLX, LSA, LSB
319 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
320 $ J, JA, JC, JE, JR, JW, NA, NW
321 DOUBLE PRECISION ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
322 $ BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A,
323 $ CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA,
324 $ CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE,
325 $ SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX,
326 $ XSCALE
327* ..
328* .. Local Arrays ..
329 DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
330 $ SUMP( 2, 2 )
331* ..
332* .. External Functions ..
333 LOGICAL LSAME
334 DOUBLE PRECISION DLAMCH
335 EXTERNAL lsame, dlamch
336* ..
337* .. External Subroutines ..
338 EXTERNAL dgemv, dlacpy, dlag2, dlaln2,
339 $ xerbla
340* ..
341* .. Intrinsic Functions ..
342 INTRINSIC abs, max, min
343* ..
344* .. Executable Statements ..
345*
346* Decode and Test the input parameters
347*
348 IF( lsame( howmny, 'A' ) ) THEN
349 ihwmny = 1
350 ilall = .true.
351 ilback = .false.
352 ELSE IF( lsame( howmny, 'S' ) ) THEN
353 ihwmny = 2
354 ilall = .false.
355 ilback = .false.
356 ELSE IF( lsame( howmny, 'B' ) ) THEN
357 ihwmny = 3
358 ilall = .true.
359 ilback = .true.
360 ELSE
361 ihwmny = -1
362 ilall = .true.
363 END IF
364*
365 IF( lsame( side, 'R' ) ) THEN
366 iside = 1
367 compl = .false.
368 compr = .true.
369 ELSE IF( lsame( side, 'L' ) ) THEN
370 iside = 2
371 compl = .true.
372 compr = .false.
373 ELSE IF( lsame( side, 'B' ) ) THEN
374 iside = 3
375 compl = .true.
376 compr = .true.
377 ELSE
378 iside = -1
379 END IF
380*
381 info = 0
382 IF( iside.LT.0 ) THEN
383 info = -1
384 ELSE IF( ihwmny.LT.0 ) THEN
385 info = -2
386 ELSE IF( n.LT.0 ) THEN
387 info = -4
388 ELSE IF( lds.LT.max( 1, n ) ) THEN
389 info = -6
390 ELSE IF( ldp.LT.max( 1, n ) ) THEN
391 info = -8
392 END IF
393 IF( info.NE.0 ) THEN
394 CALL xerbla( 'DTGEVC', -info )
395 RETURN
396 END IF
397*
398* Count the number of eigenvectors to be computed
399*
400 IF( .NOT.ilall ) THEN
401 im = 0
402 ilcplx = .false.
403 DO 10 j = 1, n
404 IF( ilcplx ) THEN
405 ilcplx = .false.
406 GO TO 10
407 END IF
408 IF( j.LT.n ) THEN
409 IF( s( j+1, j ).NE.zero )
410 $ ilcplx = .true.
411 END IF
412 IF( ilcplx ) THEN
413 IF( SELECT( j ) .OR. SELECT( j+1 ) )
414 $ im = im + 2
415 ELSE
416 IF( SELECT( j ) )
417 $ im = im + 1
418 END IF
419 10 CONTINUE
420 ELSE
421 im = n
422 END IF
423*
424* Check 2-by-2 diagonal blocks of A, B
425*
426 ilabad = .false.
427 ilbbad = .false.
428 DO 20 j = 1, n - 1
429 IF( s( j+1, j ).NE.zero ) THEN
430 IF( p( j, j ).EQ.zero .OR. p( j+1, j+1 ).EQ.zero .OR.
431 $ p( j, j+1 ).NE.zero )ilbbad = .true.
432 IF( j.LT.n-1 ) THEN
433 IF( s( j+2, j+1 ).NE.zero )
434 $ ilabad = .true.
435 END IF
436 END IF
437 20 CONTINUE
438*
439 IF( ilabad ) THEN
440 info = -5
441 ELSE IF( ilbbad ) THEN
442 info = -7
443 ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 ) THEN
444 info = -10
445 ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 ) THEN
446 info = -12
447 ELSE IF( mm.LT.im ) THEN
448 info = -13
449 END IF
450 IF( info.NE.0 ) THEN
451 CALL xerbla( 'DTGEVC', -info )
452 RETURN
453 END IF
454*
455* Quick return if possible
456*
457 m = im
458 IF( n.EQ.0 )
459 $ RETURN
460*
461* Machine Constants
462*
463 safmin = dlamch( 'Safe minimum' )
464 big = one / safmin
465 ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
466 small = safmin*n / ulp
467 big = one / small
468 bignum = one / ( safmin*n )
469*
470* Compute the 1-norm of each column of the strictly upper triangular
471* part (i.e., excluding all elements belonging to the diagonal
472* blocks) of A and B to check for possible overflow in the
473* triangular solver.
474*
475 anorm = abs( s( 1, 1 ) )
476 IF( n.GT.1 )
477 $ anorm = anorm + abs( s( 2, 1 ) )
478 bnorm = abs( p( 1, 1 ) )
479 work( 1 ) = zero
480 work( n+1 ) = zero
481*
482 DO 50 j = 2, n
483 temp = zero
484 temp2 = zero
485 IF( s( j, j-1 ).EQ.zero ) THEN
486 iend = j - 1
487 ELSE
488 iend = j - 2
489 END IF
490 DO 30 i = 1, iend
491 temp = temp + abs( s( i, j ) )
492 temp2 = temp2 + abs( p( i, j ) )
493 30 CONTINUE
494 work( j ) = temp
495 work( n+j ) = temp2
496 DO 40 i = iend + 1, min( j+1, n )
497 temp = temp + abs( s( i, j ) )
498 temp2 = temp2 + abs( p( i, j ) )
499 40 CONTINUE
500 anorm = max( anorm, temp )
501 bnorm = max( bnorm, temp2 )
502 50 CONTINUE
503*
504 ascale = one / max( anorm, safmin )
505 bscale = one / max( bnorm, safmin )
506*
507* Left eigenvectors
508*
509 IF( compl ) THEN
510 ieig = 0
511*
512* Main loop over eigenvalues
513*
514 ilcplx = .false.
515 DO 220 je = 1, n
516*
517* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
518* (b) this would be the second of a complex pair.
519* Check for complex eigenvalue, so as to be sure of which
520* entry(-ies) of SELECT to look at.
521*
522 IF( ilcplx ) THEN
523 ilcplx = .false.
524 GO TO 220
525 END IF
526 nw = 1
527 IF( je.LT.n ) THEN
528 IF( s( je+1, je ).NE.zero ) THEN
529 ilcplx = .true.
530 nw = 2
531 END IF
532 END IF
533 IF( ilall ) THEN
534 ilcomp = .true.
535 ELSE IF( ilcplx ) THEN
536 ilcomp = SELECT( je ) .OR. SELECT( je+1 )
537 ELSE
538 ilcomp = SELECT( je )
539 END IF
540 IF( .NOT.ilcomp )
541 $ GO TO 220
542*
543* Decide if (a) singular pencil, (b) real eigenvalue, or
544* (c) complex eigenvalue.
545*
546 IF( .NOT.ilcplx ) THEN
547 IF( abs( s( je, je ) ).LE.safmin .AND.
548 $ abs( p( je, je ) ).LE.safmin ) THEN
549*
550* Singular matrix pencil -- return unit eigenvector
551*
552 ieig = ieig + 1
553 DO 60 jr = 1, n
554 vl( jr, ieig ) = zero
555 60 CONTINUE
556 vl( ieig, ieig ) = one
557 GO TO 220
558 END IF
559 END IF
560*
561* Clear vector
562*
563 DO 70 jr = 1, nw*n
564 work( 2*n+jr ) = zero
565 70 CONTINUE
566* T
567* Compute coefficients in ( a A - b B ) y = 0
568* a is ACOEF
569* b is BCOEFR + i*BCOEFI
570*
571 IF( .NOT.ilcplx ) THEN
572*
573* Real eigenvalue
574*
575 temp = one / max( abs( s( je, je ) )*ascale,
576 $ abs( p( je, je ) )*bscale, safmin )
577 salfar = ( temp*s( je, je ) )*ascale
578 sbeta = ( temp*p( je, je ) )*bscale
579 acoef = sbeta*ascale
580 bcoefr = salfar*bscale
581 bcoefi = zero
582*
583* Scale to avoid underflow
584*
585 scale = one
586 lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
587 lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
588 $ small
589 IF( lsa )
590 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
591 IF( lsb )
592 $ scale = max( scale, ( small / abs( salfar ) )*
593 $ min( bnorm, big ) )
594 IF( lsa .OR. lsb ) THEN
595 scale = min( scale, one /
596 $ ( safmin*max( one, abs( acoef ),
597 $ abs( bcoefr ) ) ) )
598 IF( lsa ) THEN
599 acoef = ascale*( scale*sbeta )
600 ELSE
601 acoef = scale*acoef
602 END IF
603 IF( lsb ) THEN
604 bcoefr = bscale*( scale*salfar )
605 ELSE
606 bcoefr = scale*bcoefr
607 END IF
608 END IF
609 acoefa = abs( acoef )
610 bcoefa = abs( bcoefr )
611*
612* First component is 1
613*
614 work( 2*n+je ) = one
615 xmax = one
616 ELSE
617*
618* Complex eigenvalue
619*
620 CALL dlag2( s( je, je ), lds, p( je, je ), ldp,
621 $ safmin*safety, acoef, temp, bcoefr, temp2,
622 $ bcoefi )
623 bcoefi = -bcoefi
624 IF( bcoefi.EQ.zero ) THEN
625 info = je
626 RETURN
627 END IF
628*
629* Scale to avoid over/underflow
630*
631 acoefa = abs( acoef )
632 bcoefa = abs( bcoefr ) + abs( bcoefi )
633 scale = one
634 IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
635 $ scale = ( safmin / ulp ) / acoefa
636 IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
637 $ scale = max( scale, ( safmin / ulp ) / bcoefa )
638 IF( safmin*acoefa.GT.ascale )
639 $ scale = ascale / ( safmin*acoefa )
640 IF( safmin*bcoefa.GT.bscale )
641 $ scale = min( scale, bscale / ( safmin*bcoefa ) )
642 IF( scale.NE.one ) THEN
643 acoef = scale*acoef
644 acoefa = abs( acoef )
645 bcoefr = scale*bcoefr
646 bcoefi = scale*bcoefi
647 bcoefa = abs( bcoefr ) + abs( bcoefi )
648 END IF
649*
650* Compute first two components of eigenvector
651*
652 temp = acoef*s( je+1, je )
653 temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
654 temp2i = -bcoefi*p( je, je )
655 IF( abs( temp ).GT.abs( temp2r )+abs( temp2i ) ) THEN
656 work( 2*n+je ) = one
657 work( 3*n+je ) = zero
658 work( 2*n+je+1 ) = -temp2r / temp
659 work( 3*n+je+1 ) = -temp2i / temp
660 ELSE
661 work( 2*n+je+1 ) = one
662 work( 3*n+je+1 ) = zero
663 temp = acoef*s( je, je+1 )
664 work( 2*n+je ) = ( bcoefr*p( je+1, je+1 )-acoef*
665 $ s( je+1, je+1 ) ) / temp
666 work( 3*n+je ) = bcoefi*p( je+1, je+1 ) / temp
667 END IF
668 xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
669 $ abs( work( 2*n+je+1 ) )+abs( work( 3*n+je+1 ) ) )
670 END IF
671*
672 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
673*
674* T
675* Triangular solve of (a A - b B) y = 0
676*
677* T
678* (rowwise in (a A - b B) , or columnwise in (a A - b B) )
679*
680 il2by2 = .false.
681*
682 DO 160 j = je + nw, n
683 IF( il2by2 ) THEN
684 il2by2 = .false.
685 GO TO 160
686 END IF
687*
688 na = 1
689 bdiag( 1 ) = p( j, j )
690 IF( j.LT.n ) THEN
691 IF( s( j+1, j ).NE.zero ) THEN
692 il2by2 = .true.
693 bdiag( 2 ) = p( j+1, j+1 )
694 na = 2
695 END IF
696 END IF
697*
698* Check whether scaling is necessary for dot products
699*
700 xscale = one / max( one, xmax )
701 temp = max( work( j ), work( n+j ),
702 $ acoefa*work( j )+bcoefa*work( n+j ) )
703 IF( il2by2 )
704 $ temp = max( temp, work( j+1 ), work( n+j+1 ),
705 $ acoefa*work( j+1 )+bcoefa*work( n+j+1 ) )
706 IF( temp.GT.bignum*xscale ) THEN
707 DO 90 jw = 0, nw - 1
708 DO 80 jr = je, j - 1
709 work( ( jw+2 )*n+jr ) = xscale*
710 $ work( ( jw+2 )*n+jr )
711 80 CONTINUE
712 90 CONTINUE
713 xmax = xmax*xscale
714 END IF
715*
716* Compute dot products
717*
718* j-1
719* SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
720* k=je
721*
722* To reduce the op count, this is done as
723*
724* _ j-1 _ j-1
725* a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) )
726* k=je k=je
727*
728* which may cause underflow problems if A or B are close
729* to underflow. (E.g., less than SMALL.)
730*
731*
732 DO 120 jw = 1, nw
733 DO 110 ja = 1, na
734 sums( ja, jw ) = zero
735 sump( ja, jw ) = zero
736*
737 DO 100 jr = je, j - 1
738 sums( ja, jw ) = sums( ja, jw ) +
739 $ s( jr, j+ja-1 )*
740 $ work( ( jw+1 )*n+jr )
741 sump( ja, jw ) = sump( ja, jw ) +
742 $ p( jr, j+ja-1 )*
743 $ work( ( jw+1 )*n+jr )
744 100 CONTINUE
745 110 CONTINUE
746 120 CONTINUE
747*
748 DO 130 ja = 1, na
749 IF( ilcplx ) THEN
750 sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
751 $ bcoefr*sump( ja, 1 ) -
752 $ bcoefi*sump( ja, 2 )
753 sum( ja, 2 ) = -acoef*sums( ja, 2 ) +
754 $ bcoefr*sump( ja, 2 ) +
755 $ bcoefi*sump( ja, 1 )
756 ELSE
757 sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
758 $ bcoefr*sump( ja, 1 )
759 END IF
760 130 CONTINUE
761*
762* T
763* Solve ( a A - b B ) y = SUM(,)
764* with scaling and perturbation of the denominator
765*
766 CALL dlaln2( .true., na, nw, dmin, acoef, s( j, j ),
767 $ lds,
768 $ bdiag( 1 ), bdiag( 2 ), sum, 2, bcoefr,
769 $ bcoefi, work( 2*n+j ), n, scale, temp,
770 $ iinfo )
771 IF( scale.LT.one ) THEN
772 DO 150 jw = 0, nw - 1
773 DO 140 jr = je, j - 1
774 work( ( jw+2 )*n+jr ) = scale*
775 $ work( ( jw+2 )*n+jr )
776 140 CONTINUE
777 150 CONTINUE
778 xmax = scale*xmax
779 END IF
780 xmax = max( xmax, temp )
781 160 CONTINUE
782*
783* Copy eigenvector to VL, back transforming if
784* HOWMNY='B'.
785*
786 ieig = ieig + 1
787 IF( ilback ) THEN
788 DO 170 jw = 0, nw - 1
789 CALL dgemv( 'N', n, n+1-je, one, vl( 1, je ), ldvl,
790 $ work( ( jw+2 )*n+je ), 1, zero,
791 $ work( ( jw+4 )*n+1 ), 1 )
792 170 CONTINUE
793 CALL dlacpy( ' ', n, nw, work( 4*n+1 ), n, vl( 1,
794 $ je ),
795 $ ldvl )
796 ibeg = 1
797 ELSE
798 CALL dlacpy( ' ', n, nw, work( 2*n+1 ), n, vl( 1,
799 $ ieig ),
800 $ ldvl )
801 ibeg = je
802 END IF
803*
804* Scale eigenvector
805*
806 xmax = zero
807 IF( ilcplx ) THEN
808 DO 180 j = ibeg, n
809 xmax = max( xmax, abs( vl( j, ieig ) )+
810 $ abs( vl( j, ieig+1 ) ) )
811 180 CONTINUE
812 ELSE
813 DO 190 j = ibeg, n
814 xmax = max( xmax, abs( vl( j, ieig ) ) )
815 190 CONTINUE
816 END IF
817*
818 IF( xmax.GT.safmin ) THEN
819 xscale = one / xmax
820*
821 DO 210 jw = 0, nw - 1
822 DO 200 jr = ibeg, n
823 vl( jr, ieig+jw ) = xscale*vl( jr, ieig+jw )
824 200 CONTINUE
825 210 CONTINUE
826 END IF
827 ieig = ieig + nw - 1
828*
829 220 CONTINUE
830 END IF
831*
832* Right eigenvectors
833*
834 IF( compr ) THEN
835 ieig = im + 1
836*
837* Main loop over eigenvalues
838*
839 ilcplx = .false.
840 DO 500 je = n, 1, -1
841*
842* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
843* (b) this would be the second of a complex pair.
844* Check for complex eigenvalue, so as to be sure of which
845* entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
846* or SELECT(JE-1).
847* If this is a complex pair, the 2-by-2 diagonal block
848* corresponding to the eigenvalue is in rows/columns JE-1:JE
849*
850 IF( ilcplx ) THEN
851 ilcplx = .false.
852 GO TO 500
853 END IF
854 nw = 1
855 IF( je.GT.1 ) THEN
856 IF( s( je, je-1 ).NE.zero ) THEN
857 ilcplx = .true.
858 nw = 2
859 END IF
860 END IF
861 IF( ilall ) THEN
862 ilcomp = .true.
863 ELSE IF( ilcplx ) THEN
864 ilcomp = SELECT( je ) .OR. SELECT( je-1 )
865 ELSE
866 ilcomp = SELECT( je )
867 END IF
868 IF( .NOT.ilcomp )
869 $ GO TO 500
870*
871* Decide if (a) singular pencil, (b) real eigenvalue, or
872* (c) complex eigenvalue.
873*
874 IF( .NOT.ilcplx ) THEN
875 IF( abs( s( je, je ) ).LE.safmin .AND.
876 $ abs( p( je, je ) ).LE.safmin ) THEN
877*
878* Singular matrix pencil -- unit eigenvector
879*
880 ieig = ieig - 1
881 DO 230 jr = 1, n
882 vr( jr, ieig ) = zero
883 230 CONTINUE
884 vr( ieig, ieig ) = one
885 GO TO 500
886 END IF
887 END IF
888*
889* Clear vector
890*
891 DO 250 jw = 0, nw - 1
892 DO 240 jr = 1, n
893 work( ( jw+2 )*n+jr ) = zero
894 240 CONTINUE
895 250 CONTINUE
896*
897* Compute coefficients in ( a A - b B ) x = 0
898* a is ACOEF
899* b is BCOEFR + i*BCOEFI
900*
901 IF( .NOT.ilcplx ) THEN
902*
903* Real eigenvalue
904*
905 temp = one / max( abs( s( je, je ) )*ascale,
906 $ abs( p( je, je ) )*bscale, safmin )
907 salfar = ( temp*s( je, je ) )*ascale
908 sbeta = ( temp*p( je, je ) )*bscale
909 acoef = sbeta*ascale
910 bcoefr = salfar*bscale
911 bcoefi = zero
912*
913* Scale to avoid underflow
914*
915 scale = one
916 lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
917 lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
918 $ small
919 IF( lsa )
920 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
921 IF( lsb )
922 $ scale = max( scale, ( small / abs( salfar ) )*
923 $ min( bnorm, big ) )
924 IF( lsa .OR. lsb ) THEN
925 scale = min( scale, one /
926 $ ( safmin*max( one, abs( acoef ),
927 $ abs( bcoefr ) ) ) )
928 IF( lsa ) THEN
929 acoef = ascale*( scale*sbeta )
930 ELSE
931 acoef = scale*acoef
932 END IF
933 IF( lsb ) THEN
934 bcoefr = bscale*( scale*salfar )
935 ELSE
936 bcoefr = scale*bcoefr
937 END IF
938 END IF
939 acoefa = abs( acoef )
940 bcoefa = abs( bcoefr )
941*
942* First component is 1
943*
944 work( 2*n+je ) = one
945 xmax = one
946*
947* Compute contribution from column JE of A and B to sum
948* (See "Further Details", above.)
949*
950 DO 260 jr = 1, je - 1
951 work( 2*n+jr ) = bcoefr*p( jr, je ) -
952 $ acoef*s( jr, je )
953 260 CONTINUE
954 ELSE
955*
956* Complex eigenvalue
957*
958 CALL dlag2( s( je-1, je-1 ), lds, p( je-1, je-1 ),
959 $ ldp,
960 $ safmin*safety, acoef, temp, bcoefr, temp2,
961 $ bcoefi )
962 IF( bcoefi.EQ.zero ) THEN
963 info = je - 1
964 RETURN
965 END IF
966*
967* Scale to avoid over/underflow
968*
969 acoefa = abs( acoef )
970 bcoefa = abs( bcoefr ) + abs( bcoefi )
971 scale = one
972 IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
973 $ scale = ( safmin / ulp ) / acoefa
974 IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
975 $ scale = max( scale, ( safmin / ulp ) / bcoefa )
976 IF( safmin*acoefa.GT.ascale )
977 $ scale = ascale / ( safmin*acoefa )
978 IF( safmin*bcoefa.GT.bscale )
979 $ scale = min( scale, bscale / ( safmin*bcoefa ) )
980 IF( scale.NE.one ) THEN
981 acoef = scale*acoef
982 acoefa = abs( acoef )
983 bcoefr = scale*bcoefr
984 bcoefi = scale*bcoefi
985 bcoefa = abs( bcoefr ) + abs( bcoefi )
986 END IF
987*
988* Compute first two components of eigenvector
989* and contribution to sums
990*
991 temp = acoef*s( je, je-1 )
992 temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
993 temp2i = -bcoefi*p( je, je )
994 IF( abs( temp ).GE.abs( temp2r )+abs( temp2i ) ) THEN
995 work( 2*n+je ) = one
996 work( 3*n+je ) = zero
997 work( 2*n+je-1 ) = -temp2r / temp
998 work( 3*n+je-1 ) = -temp2i / temp
999 ELSE
1000 work( 2*n+je-1 ) = one
1001 work( 3*n+je-1 ) = zero
1002 temp = acoef*s( je-1, je )
1003 work( 2*n+je ) = ( bcoefr*p( je-1, je-1 )-acoef*
1004 $ s( je-1, je-1 ) ) / temp
1005 work( 3*n+je ) = bcoefi*p( je-1, je-1 ) / temp
1006 END IF
1007*
1008 xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
1009 $ abs( work( 2*n+je-1 ) )+abs( work( 3*n+je-1 ) ) )
1010*
1011* Compute contribution from columns JE and JE-1
1012* of A and B to the sums.
1013*
1014 creala = acoef*work( 2*n+je-1 )
1015 cimaga = acoef*work( 3*n+je-1 )
1016 crealb = bcoefr*work( 2*n+je-1 ) -
1017 $ bcoefi*work( 3*n+je-1 )
1018 cimagb = bcoefi*work( 2*n+je-1 ) +
1019 $ bcoefr*work( 3*n+je-1 )
1020 cre2a = acoef*work( 2*n+je )
1021 cim2a = acoef*work( 3*n+je )
1022 cre2b = bcoefr*work( 2*n+je ) - bcoefi*work( 3*n+je )
1023 cim2b = bcoefi*work( 2*n+je ) + bcoefr*work( 3*n+je )
1024 DO 270 jr = 1, je - 2
1025 work( 2*n+jr ) = -creala*s( jr, je-1 ) +
1026 $ crealb*p( jr, je-1 ) -
1027 $ cre2a*s( jr, je ) + cre2b*p( jr, je )
1028 work( 3*n+jr ) = -cimaga*s( jr, je-1 ) +
1029 $ cimagb*p( jr, je-1 ) -
1030 $ cim2a*s( jr, je ) + cim2b*p( jr, je )
1031 270 CONTINUE
1032 END IF
1033*
1034 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
1035*
1036* Columnwise triangular solve of (a A - b B) x = 0
1037*
1038 il2by2 = .false.
1039 DO 370 j = je - nw, 1, -1
1040*
1041* If a 2-by-2 block, is in position j-1:j, wait until
1042* next iteration to process it (when it will be j:j+1)
1043*
1044 IF( .NOT.il2by2 .AND. j.GT.1 ) THEN
1045 IF( s( j, j-1 ).NE.zero ) THEN
1046 il2by2 = .true.
1047 GO TO 370
1048 END IF
1049 END IF
1050 bdiag( 1 ) = p( j, j )
1051 IF( il2by2 ) THEN
1052 na = 2
1053 bdiag( 2 ) = p( j+1, j+1 )
1054 ELSE
1055 na = 1
1056 END IF
1057*
1058* Compute x(j) (and x(j+1), if 2-by-2 block)
1059*
1060 CALL dlaln2( .false., na, nw, dmin, acoef, s( j, j ),
1061 $ lds, bdiag( 1 ), bdiag( 2 ), work( 2*n+j ),
1062 $ n, bcoefr, bcoefi, sum, 2, scale, temp,
1063 $ iinfo )
1064 IF( scale.LT.one ) THEN
1065*
1066 DO 290 jw = 0, nw - 1
1067 DO 280 jr = 1, je
1068 work( ( jw+2 )*n+jr ) = scale*
1069 $ work( ( jw+2 )*n+jr )
1070 280 CONTINUE
1071 290 CONTINUE
1072 END IF
1073 xmax = max( scale*xmax, temp )
1074*
1075 DO 310 jw = 1, nw
1076 DO 300 ja = 1, na
1077 work( ( jw+1 )*n+j+ja-1 ) = sum( ja, jw )
1078 300 CONTINUE
1079 310 CONTINUE
1080*
1081* w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
1082*
1083 IF( j.GT.1 ) THEN
1084*
1085* Check whether scaling is necessary for sum.
1086*
1087 xscale = one / max( one, xmax )
1088 temp = acoefa*work( j ) + bcoefa*work( n+j )
1089 IF( il2by2 )
1090 $ temp = max( temp, acoefa*work( j+1 )+bcoefa*
1091 $ work( n+j+1 ) )
1092 temp = max( temp, acoefa, bcoefa )
1093 IF( temp.GT.bignum*xscale ) THEN
1094*
1095 DO 330 jw = 0, nw - 1
1096 DO 320 jr = 1, je
1097 work( ( jw+2 )*n+jr ) = xscale*
1098 $ work( ( jw+2 )*n+jr )
1099 320 CONTINUE
1100 330 CONTINUE
1101 xmax = xmax*xscale
1102 END IF
1103*
1104* Compute the contributions of the off-diagonals of
1105* column j (and j+1, if 2-by-2 block) of A and B to the
1106* sums.
1107*
1108*
1109 DO 360 ja = 1, na
1110 IF( ilcplx ) THEN
1111 creala = acoef*work( 2*n+j+ja-1 )
1112 cimaga = acoef*work( 3*n+j+ja-1 )
1113 crealb = bcoefr*work( 2*n+j+ja-1 ) -
1114 $ bcoefi*work( 3*n+j+ja-1 )
1115 cimagb = bcoefi*work( 2*n+j+ja-1 ) +
1116 $ bcoefr*work( 3*n+j+ja-1 )
1117 DO 340 jr = 1, j - 1
1118 work( 2*n+jr ) = work( 2*n+jr ) -
1119 $ creala*s( jr, j+ja-1 ) +
1120 $ crealb*p( jr, j+ja-1 )
1121 work( 3*n+jr ) = work( 3*n+jr ) -
1122 $ cimaga*s( jr, j+ja-1 ) +
1123 $ cimagb*p( jr, j+ja-1 )
1124 340 CONTINUE
1125 ELSE
1126 creala = acoef*work( 2*n+j+ja-1 )
1127 crealb = bcoefr*work( 2*n+j+ja-1 )
1128 DO 350 jr = 1, j - 1
1129 work( 2*n+jr ) = work( 2*n+jr ) -
1130 $ creala*s( jr, j+ja-1 ) +
1131 $ crealb*p( jr, j+ja-1 )
1132 350 CONTINUE
1133 END IF
1134 360 CONTINUE
1135 END IF
1136*
1137 il2by2 = .false.
1138 370 CONTINUE
1139*
1140* Copy eigenvector to VR, back transforming if
1141* HOWMNY='B'.
1142*
1143 ieig = ieig - nw
1144 IF( ilback ) THEN
1145*
1146 DO 410 jw = 0, nw - 1
1147 DO 380 jr = 1, n
1148 work( ( jw+4 )*n+jr ) = work( ( jw+2 )*n+1 )*
1149 $ vr( jr, 1 )
1150 380 CONTINUE
1151*
1152* A series of compiler directives to defeat
1153* vectorization for the next loop
1154*
1155*
1156 DO 400 jc = 2, je
1157 DO 390 jr = 1, n
1158 work( ( jw+4 )*n+jr ) = work( ( jw+4 )*n+jr ) +
1159 $ work( ( jw+2 )*n+jc )*vr( jr, jc )
1160 390 CONTINUE
1161 400 CONTINUE
1162 410 CONTINUE
1163*
1164 DO 430 jw = 0, nw - 1
1165 DO 420 jr = 1, n
1166 vr( jr, ieig+jw ) = work( ( jw+4 )*n+jr )
1167 420 CONTINUE
1168 430 CONTINUE
1169*
1170 iend = n
1171 ELSE
1172 DO 450 jw = 0, nw - 1
1173 DO 440 jr = 1, n
1174 vr( jr, ieig+jw ) = work( ( jw+2 )*n+jr )
1175 440 CONTINUE
1176 450 CONTINUE
1177*
1178 iend = je
1179 END IF
1180*
1181* Scale eigenvector
1182*
1183 xmax = zero
1184 IF( ilcplx ) THEN
1185 DO 460 j = 1, iend
1186 xmax = max( xmax, abs( vr( j, ieig ) )+
1187 $ abs( vr( j, ieig+1 ) ) )
1188 460 CONTINUE
1189 ELSE
1190 DO 470 j = 1, iend
1191 xmax = max( xmax, abs( vr( j, ieig ) ) )
1192 470 CONTINUE
1193 END IF
1194*
1195 IF( xmax.GT.safmin ) THEN
1196 xscale = one / xmax
1197 DO 490 jw = 0, nw - 1
1198 DO 480 jr = 1, iend
1199 vr( jr, ieig+jw ) = xscale*vr( jr, ieig+jw )
1200 480 CONTINUE
1201 490 CONTINUE
1202 END IF
1203 500 CONTINUE
1204 END IF
1205*
1206 RETURN
1207*
1208* End of DTGEVC
1209*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
dlag2
subroutine dlag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition dlag2.f:154
dlaln2
subroutine dlaln2(ltrans, na, nw, smin, ca, a, lda, d1, d2, b, ldb, wr, wi, x, ldx, scale, xnorm, info)
DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
Definition dlaln2.f:216
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: