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,
')' )
subroutine xerbla(srname, info)
subroutine dget22(transa, transe, transw, n, a, lda, e, lde, wr, wi, work, result)
DGET22
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
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 dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.