375 SUBROUTINE dget23( COMP, BALANC, JTYPE, THRESH, ISEED, NOUNIT, N,
376 $ a, lda, h, wr, wi, wr1, wi1, vl, ldvl, vr,
377 $ ldvr, lre, ldlre, rcondv, rcndv1, rcdvin,
378 $ rconde, rcnde1, rcdein, scale, scale1, result,
379 $ work, lwork, iwork, info )
389 INTEGER INFO, JTYPE, LDA, LDLRE, LDVL, LDVR, LWORK, N,
391 DOUBLE PRECISION THRESH
394 INTEGER ISEED( 4 ), IWORK( * )
395 DOUBLE PRECISION A( lda, * ), H( lda, * ), LRE( ldlre, * ),
396 $ rcdein( * ), rcdvin( * ), rcnde1( * ),
397 $ rcndv1( * ), rconde( * ), rcondv( * ),
398 $ result( 11 ), scale( * ), scale1( * ),
399 $ vl( ldvl, * ), vr( ldvr, * ), wi( * ),
400 $ wi1( * ), work( * ), wr( * ), wr1( * )
407 DOUBLE PRECISION ZERO, ONE, TWO
408 parameter ( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
409 DOUBLE PRECISION EPSIN
410 parameter ( epsin = 5.9605d-8 )
415 INTEGER I, IHI, IHI1, IINFO, ILO, ILO1, ISENS, ISENSM,
417 DOUBLE PRECISION ABNRM, ABNRM1, EPS, SMLNUM, TNRM, TOL, TOLIN,
418 $ ulp, ulpinv, v, vimin, vmax, vmx, vrmin, vrmx,
423 DOUBLE PRECISION DUM( 1 ), RES( 2 )
427 DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2
428 EXTERNAL lsame, dlamch, dlapy2, dnrm2
434 INTRINSIC abs, dble, max, min
437 DATA sens /
'N',
'V' /
443 nobal = lsame( balanc,
'N' )
444 balok = nobal .OR. lsame( balanc,
'P' ) .OR.
445 $ lsame( balanc,
'S' ) .OR. lsame( balanc,
'B' )
447 IF( .NOT.balok )
THEN
449 ELSE IF( thresh.LT.zero )
THEN
451 ELSE IF( nounit.LE.0 )
THEN
453 ELSE IF( n.LT.0 )
THEN
455 ELSE IF( lda.LT.1 .OR. lda.LT.n )
THEN
457 ELSE IF( ldvl.LT.1 .OR. ldvl.LT.n )
THEN
459 ELSE IF( ldvr.LT.1 .OR. ldvr.LT.n )
THEN
461 ELSE IF( ldlre.LT.1 .OR. ldlre.LT.n )
THEN
463 ELSE IF( lwork.LT.3*n .OR. ( comp .AND. lwork.LT.6*n+n*n ) )
THEN
468 CALL xerbla(
'DGET23', -info )
483 ulp = dlamch(
'Precision' )
484 smlnum = dlamch(
'S' )
489 IF( lwork.GE.6*n+n*n )
THEN
496 CALL dlacpy(
'F', n, n, a, lda, h, lda )
497 CALL dgeevx( balanc,
'V',
'V', sense, n, h, lda, wr, wi, vl, ldvl,
498 $ vr, ldvr, ilo, ihi, scale, abnrm, rconde, rcondv,
499 $ work, lwork, iwork, iinfo )
500 IF( iinfo.NE.0 )
THEN
502 IF( jtype.NE.22 )
THEN
503 WRITE( nounit, fmt = 9998 )
'DGEEVX1', iinfo, n, jtype,
506 WRITE( nounit, fmt = 9999 )
'DGEEVX1', iinfo, n, iseed( 1 )
514 CALL dget22(
'N',
'N',
'N', n, a, lda, vr, ldvr, wr, wi, work,
516 result( 1 ) = res( 1 )
520 CALL dget22(
'T',
'N',
'T', n, a, lda, vl, ldvl, wr, wi, work,
522 result( 2 ) = res( 1 )
528 IF( wi( j ).EQ.zero )
THEN
529 tnrm = dnrm2( n, vr( 1, j ), 1 )
530 ELSE IF( wi( j ).GT.zero )
THEN
531 tnrm = dlapy2( dnrm2( n, vr( 1, j ), 1 ),
532 $ dnrm2( n, vr( 1, j+1 ), 1 ) )
534 result( 3 ) = max( result( 3 ),
535 $ min( ulpinv, abs( tnrm-one ) / ulp ) )
536 IF( wi( j ).GT.zero )
THEN
540 vtst = dlapy2( vr( jj, j ), vr( jj, j+1 ) )
543 IF( vr( jj, j+1 ).EQ.zero .AND. abs( vr( jj, j ) ).GT.
544 $ vrmx )vrmx = abs( vr( jj, j ) )
546 IF( vrmx / vmx.LT.one-two*ulp )
547 $ result( 3 ) = ulpinv
555 IF( wi( j ).EQ.zero )
THEN
556 tnrm = dnrm2( n, vl( 1, j ), 1 )
557 ELSE IF( wi( j ).GT.zero )
THEN
558 tnrm = dlapy2( dnrm2( n, vl( 1, j ), 1 ),
559 $ dnrm2( n, vl( 1, j+1 ), 1 ) )
561 result( 4 ) = max( result( 4 ),
562 $ min( ulpinv, abs( tnrm-one ) / ulp ) )
563 IF( wi( j ).GT.zero )
THEN
567 vtst = dlapy2( vl( jj, j ), vl( jj, j+1 ) )
570 IF( vl( jj, j+1 ).EQ.zero .AND. abs( vl( jj, j ) ).GT.
571 $ vrmx )vrmx = abs( vl( jj, j ) )
573 IF( vrmx / vmx.LT.one-two*ulp )
574 $ result( 4 ) = ulpinv
580 DO 200 isens = 1, isensm
582 sense = sens( isens )
586 CALL dlacpy(
'F', n, n, a, lda, h, lda )
587 CALL dgeevx( balanc,
'N',
'N', sense, n, h, lda, wr1, wi1, dum,
588 $ 1, dum, 1, ilo1, ihi1, scale1, abnrm1, rcnde1,
589 $ rcndv1, work, lwork, iwork, iinfo )
590 IF( iinfo.NE.0 )
THEN
592 IF( jtype.NE.22 )
THEN
593 WRITE( nounit, fmt = 9998 )
'DGEEVX2', iinfo, n, jtype,
596 WRITE( nounit, fmt = 9999 )
'DGEEVX2', iinfo, n,
606 IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
607 $ result( 5 ) = ulpinv
612 IF( .NOT.nobal )
THEN
614 IF( scale( j ).NE.scale1( j ) )
615 $ result( 8 ) = ulpinv
618 $ result( 8 ) = ulpinv
620 $ result( 8 ) = ulpinv
621 IF( abnrm.NE.abnrm1 )
622 $ result( 8 ) = ulpinv
627 IF( isens.EQ.2 .AND. n.GT.1 )
THEN
629 IF( rcondv( j ).NE.rcndv1( j ) )
630 $ result( 9 ) = ulpinv
636 CALL dlacpy(
'F', n, n, a, lda, h, lda )
637 CALL dgeevx( balanc,
'N',
'V', sense, n, h, lda, wr1, wi1, dum,
638 $ 1, lre, ldlre, ilo1, ihi1, scale1, abnrm1, rcnde1,
639 $ rcndv1, work, lwork, iwork, iinfo )
640 IF( iinfo.NE.0 )
THEN
642 IF( jtype.NE.22 )
THEN
643 WRITE( nounit, fmt = 9998 )
'DGEEVX3', iinfo, n, jtype,
646 WRITE( nounit, fmt = 9999 )
'DGEEVX3', iinfo, n,
656 IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
657 $ result( 5 ) = ulpinv
664 IF( vr( j, jj ).NE.lre( j, jj ) )
665 $ result( 6 ) = ulpinv
671 IF( .NOT.nobal )
THEN
673 IF( scale( j ).NE.scale1( j ) )
674 $ result( 8 ) = ulpinv
677 $ result( 8 ) = ulpinv
679 $ result( 8 ) = ulpinv
680 IF( abnrm.NE.abnrm1 )
681 $ result( 8 ) = ulpinv
686 IF( isens.EQ.2 .AND. n.GT.1 )
THEN
688 IF( rcondv( j ).NE.rcndv1( j ) )
689 $ result( 9 ) = ulpinv
695 CALL dlacpy(
'F', n, n, a, lda, h, lda )
696 CALL dgeevx( balanc,
'V',
'N', sense, n, h, lda, wr1, wi1, lre,
697 $ ldlre, dum, 1, ilo1, ihi1, scale1, abnrm1, rcnde1,
698 $ rcndv1, work, lwork, iwork, iinfo )
699 IF( iinfo.NE.0 )
THEN
701 IF( jtype.NE.22 )
THEN
702 WRITE( nounit, fmt = 9998 )
'DGEEVX4', iinfo, n, jtype,
705 WRITE( nounit, fmt = 9999 )
'DGEEVX4', iinfo, n,
715 IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
716 $ result( 5 ) = ulpinv
723 IF( vl( j, jj ).NE.lre( j, jj ) )
724 $ result( 7 ) = ulpinv
730 IF( .NOT.nobal )
THEN
732 IF( scale( j ).NE.scale1( j ) )
733 $ result( 8 ) = ulpinv
736 $ result( 8 ) = ulpinv
738 $ result( 8 ) = ulpinv
739 IF( abnrm.NE.abnrm1 )
740 $ result( 8 ) = ulpinv
745 IF( isens.EQ.2 .AND. n.GT.1 )
THEN
747 IF( rcondv( j ).NE.rcndv1( j ) )
748 $ result( 9 ) = ulpinv
759 CALL dlacpy(
'F', n, n, a, lda, h, lda )
760 CALL dgeevx(
'N',
'V',
'V',
'B', n, h, lda, wr, wi, vl, ldvl,
761 $ vr, ldvr, ilo, ihi, scale, abnrm, rconde, rcondv,
762 $ work, lwork, iwork, iinfo )
763 IF( iinfo.NE.0 )
THEN
765 WRITE( nounit, fmt = 9999 )
'DGEEVX5', iinfo, n, iseed( 1 )
778 IF( wr( j ).LT.vrmin )
THEN
788 vrmin = rconde( kmin )
789 rconde( kmin ) = rconde( i )
791 vrmin = rcondv( kmin )
792 rcondv( kmin ) = rcondv( i )
800 eps = max( epsin, ulp )
801 v = max( dble( n )*eps*abnrm, smlnum )
805 IF( v.GT.rcondv( i )*rconde( i ) )
THEN
808 tol = v / rconde( i )
810 IF( v.GT.rcdvin( i )*rcdein( i ) )
THEN
813 tolin = v / rcdein( i )
815 tol = max( tol, smlnum / eps )
816 tolin = max( tolin, smlnum / eps )
817 IF( eps*( rcdvin( i )-tolin ).GT.rcondv( i )+tol )
THEN
819 ELSE IF( rcdvin( i )-tolin.GT.rcondv( i )+tol )
THEN
820 vmax = ( rcdvin( i )-tolin ) / ( rcondv( i )+tol )
821 ELSE IF( rcdvin( i )+tolin.LT.eps*( rcondv( i )-tol ) )
THEN
823 ELSE IF( rcdvin( i )+tolin.LT.rcondv( i )-tol )
THEN
824 vmax = ( rcondv( i )-tol ) / ( rcdvin( i )+tolin )
828 result( 10 ) = max( result( 10 ), vmax )
836 IF( v.GT.rcondv( i ) )
THEN
839 tol = v / rcondv( i )
841 IF( v.GT.rcdvin( i ) )
THEN
844 tolin = v / rcdvin( i )
846 tol = max( tol, smlnum / eps )
847 tolin = max( tolin, smlnum / eps )
848 IF( eps*( rcdein( i )-tolin ).GT.rconde( i )+tol )
THEN
850 ELSE IF( rcdein( i )-tolin.GT.rconde( i )+tol )
THEN
851 vmax = ( rcdein( i )-tolin ) / ( rconde( i )+tol )
852 ELSE IF( rcdein( i )+tolin.LT.eps*( rconde( i )-tol ) )
THEN
854 ELSE IF( rcdein( i )+tolin.LT.rconde( i )-tol )
THEN
855 vmax = ( rconde( i )-tol ) / ( rcdein( i )+tolin )
859 result( 11 ) = max( result( 11 ), vmax )
865 9999
FORMAT(
' DGET23: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
866 $ i6,
', INPUT EXAMPLE NUMBER = ', i4 )
867 9998
FORMAT(
' DGET23: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
868 $ i6,
', JTYPE=', i6,
', BALANC = ', a,
', ISEED=(',
869 $ 3( i5,
',' ), i5,
')' )
subroutine dgeevx(BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV, WORK, LWORK, IWORK, INFO)
DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
subroutine dget22(TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR, WI, WORK, RESULT)
DGET22
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dget23(COMP, BALANC, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, VR, LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN, RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT, WORK, LWORK, IWORK, INFO)
DGET23