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 )
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' )
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
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
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 xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
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.
real function slamch(CMACH)
SLAMCH
real function snrm2(N, X, INCX)
SNRM2
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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 ...
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).