375 SUBROUTINE sget23( 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,
394 INTEGER ISEED( 4 ), IWORK( * )
395 REAL 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( * )
408 parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
410 parameter ( epsin = 5.9605e-8 )
415 INTEGER I, IHI, IHI1, IINFO, ILO, ILO1, ISENS, ISENSM,
417 REAL ABNRM, ABNRM1, EPS, SMLNUM, TNRM, TOL, TOLIN,
418 $ ulp, ulpinv, v, vimin, vmax, vmx, vrmin, vrmx,
423 REAL DUM( 1 ), RES( 2 )
427 REAL SLAMCH, SLAPY2, SNRM2
428 EXTERNAL lsame, slamch, slapy2, snrm2
434 INTRINSIC abs, max, min, real
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(
'SGET23', -info )
483 ulp = slamch(
'Precision' )
484 smlnum = slamch(
'S' )
489 IF( lwork.GE.6*n+n*n )
THEN
496 CALL slacpy(
'F', n, n, a, lda, h, lda )
497 CALL sgeevx( 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 )
'SGEEVX1', iinfo, n, jtype,
506 WRITE( nounit, fmt = 9999 )
'SGEEVX1', iinfo, n, iseed( 1 )
514 CALL sget22(
'N',
'N',
'N', n, a, lda, vr, ldvr, wr, wi, work,
516 result( 1 ) = res( 1 )
520 CALL sget22(
'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 = snrm2( n, vr( 1, j ), 1 )
530 ELSE IF( wi( j ).GT.zero )
THEN
531 tnrm = slapy2( snrm2( n, vr( 1, j ), 1 ),
532 $ snrm2( 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 = slapy2( 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 = snrm2( n, vl( 1, j ), 1 )
557 ELSE IF( wi( j ).GT.zero )
THEN
558 tnrm = slapy2( snrm2( n, vl( 1, j ), 1 ),
559 $ snrm2( 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 = slapy2( 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 slacpy(
'F', n, n, a, lda, h, lda )
587 CALL sgeevx( 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 )
'SGEEVX2', iinfo, n, jtype,
596 WRITE( nounit, fmt = 9999 )
'SGEEVX2', 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 slacpy(
'F', n, n, a, lda, h, lda )
637 CALL sgeevx( 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 )
'SGEEVX3', iinfo, n, jtype,
646 WRITE( nounit, fmt = 9999 )
'SGEEVX3', 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 slacpy(
'F', n, n, a, lda, h, lda )
696 CALL sgeevx( 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 )
'SGEEVX4', iinfo, n, jtype,
705 WRITE( nounit, fmt = 9999 )
'SGEEVX4', 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 slacpy(
'F', n, n, a, lda, h, lda )
760 CALL sgeevx(
'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 )
'SGEEVX5', 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(
REAL( 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(
' SGET23: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
866 $ i6,
', INPUT EXAMPLE NUMBER = ', i4 )
867 9998
FORMAT(
' SGET23: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
868 $ i6,
', JTYPE=', i6,
', BALANC = ', a,
', ISEED=(',
869 $ 3( i5,
',' ), i5,
')' )
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
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 xerbla(SRNAME, INFO)
XERBLA
subroutine sget22(TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR, WI, WORK, RESULT)
SGET22
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.