373 SUBROUTINE sget23( 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,
391 INTEGER ISEED( 4 ), IWORK( * )
392 REAL 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( * )
405 PARAMETER ( ZERO = 0.0e0, one = 1.0e0, two = 2.0e0 )
407 PARAMETER ( EPSIN = 5.9605e-8 )
412 INTEGER I, IHI, IHI1, IINFO, ILO, ILO1, ISENS, ISENSM,
414 REAL ABNRM, ABNRM1, EPS, SMLNUM, TNRM, TOL, TOLIN,
415 $ ulp, ulpinv, v, vimin, vmax, vmx, vrmin, vrmx,
420 REAL DUM( 1 ), RES( 2 )
424 REAL SLAMCH, SLAPY2, SNRM2
425 EXTERNAL LSAME, SLAMCH, SLAPY2, SNRM2
431 INTRINSIC abs, max, min, real
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(
'SGET23', -info )
480 ulp = slamch(
'Precision' )
481 smlnum = slamch(
'S' )
486 IF( lwork.GE.6*n+n*n )
THEN
493 CALL slacpy(
'F', n, n, a, lda, h, lda )
494 CALL sgeevx( 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 )
'SGEEVX1', iinfo, n, jtype,
503 WRITE( nounit, fmt = 9999 )
'SGEEVX1', iinfo, n, iseed( 1 )
511 CALL sget22(
'N',
'N',
'N', n, a, lda, vr, ldvr, wr, wi, work,
513 result( 1 ) = res( 1 )
517 CALL sget22(
'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 = snrm2( n, vr( 1, j ), 1 )
527 ELSE IF( wi( j ).GT.zero )
THEN
528 tnrm = slapy2( snrm2( n, vr( 1, j ), 1 ),
529 $ snrm2( 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 = slapy2( 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 = snrm2( n, vl( 1, j ), 1 )
554 ELSE IF( wi( j ).GT.zero )
THEN
555 tnrm = slapy2( snrm2( n, vl( 1, j ), 1 ),
556 $ snrm2( 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 = slapy2( 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 slacpy(
'F', n, n, a, lda, h, lda )
584 CALL sgeevx( 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 )
'SGEEVX2', iinfo, n, jtype,
593 WRITE( nounit, fmt = 9999 )
'SGEEVX2', 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 slacpy(
'F', n, n, a, lda, h, lda )
634 CALL sgeevx( 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 )
'SGEEVX3', iinfo, n, jtype,
643 WRITE( nounit, fmt = 9999 )
'SGEEVX3', 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 slacpy(
'F', n, n, a, lda, h, lda )
693 CALL sgeevx( 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 )
'SGEEVX4', iinfo, n, jtype,
702 WRITE( nounit, fmt = 9999 )
'SGEEVX4', 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 slacpy(
'F', n, n, a, lda, h, lda )
757 CALL sgeevx(
'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 )
'SGEEVX5', 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( real( 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(
' SGET23: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
863 $ i6,
', INPUT EXAMPLE NUMBER = ', i4 )
864 9998
FORMAT(
' SGET23: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
865 $ i6,
', JTYPE=', i6,
', BALANC = ', a,
', ISEED=(',
866 $ 3( i5,
',' ), i5,
')' )
subroutine xerbla(srname, info)
subroutine sgeevx(balanc, jobvl, jobvr, sense, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, ilo, ihi, scale, abnrm, rconde, rcondv, work, lwork, iwork, info)
SGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sget22(transa, transe, transw, n, a, lda, e, lde, wr, wi, work, result)
SGET22
subroutine sget23(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)
SGET23