373 SUBROUTINE dget23( COMP, BALANC, JTYPE, THRESH, ISEED, NOUNIT, N,
374 $ A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, VR,
375 $ LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN,
376 $ RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT,
377 $ WORK, LWORK, IWORK, INFO )
386 INTEGER INFO, JTYPE, LDA, LDLRE, LDVL, LDVR, LWORK, N,
388 DOUBLE PRECISION THRESH
391 INTEGER ISEED( 4 ), IWORK( * )
392 DOUBLE PRECISION A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
393 $ RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
394 $ RCNDV1( * ), RCONDE( * ), RCONDV( * ),
395 $ result( 11 ), scale( * ), scale1( * ),
396 $ vl( ldvl, * ), vr( ldvr, * ), wi( * ),
397 $ wi1( * ), work( * ), wr( * ), wr1( * )
404 DOUBLE PRECISION ZERO, ONE, TWO
405 PARAMETER ( ZERO = 0.0d0, one = 1.0d0, two = 2.0d0 )
406 DOUBLE PRECISION EPSIN
407 PARAMETER ( EPSIN = 5.9605d-8 )
412 INTEGER I, IHI, IHI1, IINFO, ILO, ILO1, ISENS, ISENSM,
414 DOUBLE PRECISION ABNRM, ABNRM1, EPS, SMLNUM, TNRM, TOL, TOLIN,
415 $ ulp, ulpinv, v, vimin, vmax, vmx, vrmin, vrmx,
420 DOUBLE PRECISION DUM( 1 ), RES( 2 )
424 DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2
425 EXTERNAL LSAME, DLAMCH, DLAPY2, DNRM2
431 INTRINSIC abs, dble, max, min
434 DATA sens /
'N',
'V' /
440 nobal = lsame( balanc,
'N' )
441 balok = nobal .OR. lsame( balanc,
'P' ) .OR.
442 $ lsame( balanc,
'S' ) .OR. lsame( balanc,
'B' )
444 IF( .NOT.balok )
THEN
446 ELSE IF( thresh.LT.zero )
THEN
448 ELSE IF( nounit.LE.0 )
THEN
450 ELSE IF( n.LT.0 )
THEN
452 ELSE IF( lda.LT.1 .OR. lda.LT.n )
THEN
454 ELSE IF( ldvl.LT.1 .OR. ldvl.LT.n )
THEN
456 ELSE IF( ldvr.LT.1 .OR. ldvr.LT.n )
THEN
458 ELSE IF( ldlre.LT.1 .OR. ldlre.LT.n )
THEN
460 ELSE IF( lwork.LT.3*n .OR. ( comp .AND. lwork.LT.6*n+n*n ) )
THEN
465 CALL xerbla(
'DGET23', -info )
480 ulp = dlamch(
'Precision' )
481 smlnum = dlamch(
'S' )
486 IF( lwork.GE.6*n+n*n )
THEN
493 CALL dlacpy(
'F', n, n, a, lda, h, lda )
494 CALL dgeevx( balanc,
'V',
'V', sense, n, h, lda, wr, wi, vl, ldvl,
495 $ vr, ldvr, ilo, ihi, scale, abnrm, rconde, rcondv,
496 $ work, lwork, iwork, iinfo )
497 IF( iinfo.NE.0 )
THEN
499 IF( jtype.NE.22 )
THEN
500 WRITE( nounit, fmt = 9998 )
'DGEEVX1', iinfo, n, jtype,
503 WRITE( nounit, fmt = 9999 )
'DGEEVX1', iinfo, n, iseed( 1 )
511 CALL dget22(
'N',
'N',
'N', n, a, lda, vr, ldvr, wr, wi, work,
513 result( 1 ) = res( 1 )
517 CALL dget22(
'T',
'N',
'T', n, a, lda, vl, ldvl, wr, wi, work,
519 result( 2 ) = res( 1 )
525 IF( wi( j ).EQ.zero )
THEN
526 tnrm = dnrm2( n, vr( 1, j ), 1 )
527 ELSE IF( wi( j ).GT.zero )
THEN
528 tnrm = dlapy2( dnrm2( n, vr( 1, j ), 1 ),
529 $ dnrm2( n, vr( 1, j+1 ), 1 ) )
531 result( 3 ) = max( result( 3 ),
532 $ min( ulpinv, abs( tnrm-one ) / ulp ) )
533 IF( wi( j ).GT.zero )
THEN
537 vtst = dlapy2( vr( jj, j ), vr( jj, j+1 ) )
540 IF( vr( jj, j+1 ).EQ.zero .AND. abs( vr( jj, j ) ).GT.
541 $ vrmx )vrmx = abs( vr( jj, j ) )
543 IF( vrmx / vmx.LT.one-two*ulp )
544 $ result( 3 ) = ulpinv
552 IF( wi( j ).EQ.zero )
THEN
553 tnrm = dnrm2( n, vl( 1, j ), 1 )
554 ELSE IF( wi( j ).GT.zero )
THEN
555 tnrm = dlapy2( dnrm2( n, vl( 1, j ), 1 ),
556 $ dnrm2( n, vl( 1, j+1 ), 1 ) )
558 result( 4 ) = max( result( 4 ),
559 $ min( ulpinv, abs( tnrm-one ) / ulp ) )
560 IF( wi( j ).GT.zero )
THEN
564 vtst = dlapy2( vl( jj, j ), vl( jj, j+1 ) )
567 IF( vl( jj, j+1 ).EQ.zero .AND. abs( vl( jj, j ) ).GT.
568 $ vrmx )vrmx = abs( vl( jj, j ) )
570 IF( vrmx / vmx.LT.one-two*ulp )
571 $ result( 4 ) = ulpinv
577 DO 200 isens = 1, isensm
579 sense = sens( isens )
583 CALL dlacpy(
'F', n, n, a, lda, h, lda )
584 CALL dgeevx( balanc,
'N',
'N', sense, n, h, lda, wr1, wi1, dum,
585 $ 1, dum, 1, ilo1, ihi1, scale1, abnrm1, rcnde1,
586 $ rcndv1, work, lwork, iwork, iinfo )
587 IF( iinfo.NE.0 )
THEN
589 IF( jtype.NE.22 )
THEN
590 WRITE( nounit, fmt = 9998 )
'DGEEVX2', iinfo, n, jtype,
593 WRITE( nounit, fmt = 9999 )
'DGEEVX2', iinfo, n,
603 IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
604 $ result( 5 ) = ulpinv
609 IF( .NOT.nobal )
THEN
611 IF( scale( j ).NE.scale1( j ) )
612 $ result( 8 ) = ulpinv
615 $ result( 8 ) = ulpinv
617 $ result( 8 ) = ulpinv
618 IF( abnrm.NE.abnrm1 )
619 $ result( 8 ) = ulpinv
624 IF( isens.EQ.2 .AND. n.GT.1 )
THEN
626 IF( rcondv( j ).NE.rcndv1( j ) )
627 $ result( 9 ) = ulpinv
633 CALL dlacpy(
'F', n, n, a, lda, h, lda )
634 CALL dgeevx( balanc,
'N',
'V', sense, n, h, lda, wr1, wi1, dum,
635 $ 1, lre, ldlre, ilo1, ihi1, scale1, abnrm1, rcnde1,
636 $ rcndv1, work, lwork, iwork, iinfo )
637 IF( iinfo.NE.0 )
THEN
639 IF( jtype.NE.22 )
THEN
640 WRITE( nounit, fmt = 9998 )
'DGEEVX3', iinfo, n, jtype,
643 WRITE( nounit, fmt = 9999 )
'DGEEVX3', iinfo, n,
653 IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
654 $ result( 5 ) = ulpinv
661 IF( vr( j, jj ).NE.lre( j, jj ) )
662 $ result( 6 ) = ulpinv
668 IF( .NOT.nobal )
THEN
670 IF( scale( j ).NE.scale1( j ) )
671 $ result( 8 ) = ulpinv
674 $ result( 8 ) = ulpinv
676 $ result( 8 ) = ulpinv
677 IF( abnrm.NE.abnrm1 )
678 $ result( 8 ) = ulpinv
683 IF( isens.EQ.2 .AND. n.GT.1 )
THEN
685 IF( rcondv( j ).NE.rcndv1( j ) )
686 $ result( 9 ) = ulpinv
692 CALL dlacpy(
'F', n, n, a, lda, h, lda )
693 CALL dgeevx( balanc,
'V',
'N', sense, n, h, lda, wr1, wi1, lre,
694 $ ldlre, dum, 1, ilo1, ihi1, scale1, abnrm1, rcnde1,
695 $ rcndv1, work, lwork, iwork, iinfo )
696 IF( iinfo.NE.0 )
THEN
698 IF( jtype.NE.22 )
THEN
699 WRITE( nounit, fmt = 9998 )
'DGEEVX4', iinfo, n, jtype,
702 WRITE( nounit, fmt = 9999 )
'DGEEVX4', iinfo, n,
712 IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
713 $ result( 5 ) = ulpinv
720 IF( vl( j, jj ).NE.lre( j, jj ) )
721 $ result( 7 ) = ulpinv
727 IF( .NOT.nobal )
THEN
729 IF( scale( j ).NE.scale1( j ) )
730 $ result( 8 ) = ulpinv
733 $ result( 8 ) = ulpinv
735 $ result( 8 ) = ulpinv
736 IF( abnrm.NE.abnrm1 )
737 $ result( 8 ) = ulpinv
742 IF( isens.EQ.2 .AND. n.GT.1 )
THEN
744 IF( rcondv( j ).NE.rcndv1( j ) )
745 $ result( 9 ) = ulpinv
756 CALL dlacpy(
'F', n, n, a, lda, h, lda )
757 CALL dgeevx(
'N',
'V',
'V',
'B', n, h, lda, wr, wi, vl, ldvl,
758 $ vr, ldvr, ilo, ihi, scale, abnrm, rconde, rcondv,
759 $ work, lwork, iwork, iinfo )
760 IF( iinfo.NE.0 )
THEN
762 WRITE( nounit, fmt = 9999 )
'DGEEVX5', iinfo, n, iseed( 1 )
775 IF( wr( j ).LT.vrmin )
THEN
785 vrmin = rconde( kmin )
786 rconde( kmin ) = rconde( i )
788 vrmin = rcondv( kmin )
789 rcondv( kmin ) = rcondv( i )
797 eps = max( epsin, ulp )
798 v = max( dble( n )*eps*abnrm, smlnum )
802 IF( v.GT.rcondv( i )*rconde( i ) )
THEN
805 tol = v / rconde( i )
807 IF( v.GT.rcdvin( i )*rcdein( i ) )
THEN
810 tolin = v / rcdein( i )
812 tol = max( tol, smlnum / eps )
813 tolin = max( tolin, smlnum / eps )
814 IF( eps*( rcdvin( i )-tolin ).GT.rcondv( i )+tol )
THEN
816 ELSE IF( rcdvin( i )-tolin.GT.rcondv( i )+tol )
THEN
817 vmax = ( rcdvin( i )-tolin ) / ( rcondv( i )+tol )
818 ELSE IF( rcdvin( i )+tolin.LT.eps*( rcondv( i )-tol ) )
THEN
820 ELSE IF( rcdvin( i )+tolin.LT.rcondv( i )-tol )
THEN
821 vmax = ( rcondv( i )-tol ) / ( rcdvin( i )+tolin )
825 result( 10 ) = max( result( 10 ), vmax )
833 IF( v.GT.rcondv( i ) )
THEN
836 tol = v / rcondv( i )
838 IF( v.GT.rcdvin( i ) )
THEN
841 tolin = v / rcdvin( i )
843 tol = max( tol, smlnum / eps )
844 tolin = max( tolin, smlnum / eps )
845 IF( eps*( rcdein( i )-tolin ).GT.rconde( i )+tol )
THEN
847 ELSE IF( rcdein( i )-tolin.GT.rconde( i )+tol )
THEN
848 vmax = ( rcdein( i )-tolin ) / ( rconde( i )+tol )
849 ELSE IF( rcdein( i )+tolin.LT.eps*( rconde( i )-tol ) )
THEN
851 ELSE IF( rcdein( i )+tolin.LT.rconde( i )-tol )
THEN
852 vmax = ( rconde( i )-tol ) / ( rcdein( i )+tolin )
856 result( 11 ) = max( result( 11 ), vmax )
862 9999
FORMAT(
' DGET23: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
863 $ i6,
', INPUT EXAMPLE NUMBER = ', i4 )
864 9998
FORMAT(
' DGET23: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
865 $ i6,
', JTYPE=', i6,
', BALANC = ', a,
', ISEED=(',
866 $ 3( i5,
',' ), i5,
')' )