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 )
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' )
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
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
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,
')' )