sget23.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 722 行 · 第 1/2 页
F
722 行
RESULT( 1 ) = RES( 1 )
*
* Do Test (2)
*
CALL SGET22( 'T', 'N', 'T', N, A, LDA, VL, LDVL, WR, WI, WORK,
$ RES )
RESULT( 2 ) = RES( 1 )
*
* Do Test (3)
*
DO 30 J = 1, N
TNRM = ONE
IF( WI( J ).EQ.ZERO ) THEN
TNRM = SNRM2( N, VR( 1, J ), 1 )
ELSE IF( WI( J ).GT.ZERO ) THEN
TNRM = SLAPY2( SNRM2( N, VR( 1, J ), 1 ),
$ SNRM2( N, VR( 1, J+1 ), 1 ) )
END IF
RESULT( 3 ) = MAX( RESULT( 3 ),
$ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) )
IF( WI( J ).GT.ZERO ) THEN
VMX = ZERO
VRMX = ZERO
DO 20 JJ = 1, N
VTST = SLAPY2( VR( JJ, J ), VR( JJ, J+1 ) )
IF( VTST.GT.VMX )
$ VMX = VTST
IF( VR( JJ, J+1 ).EQ.ZERO .AND. ABS( VR( JJ, J ) ).GT.
$ VRMX )VRMX = ABS( VR( JJ, J ) )
20 CONTINUE
IF( VRMX / VMX.LT.ONE-TWO*ULP )
$ RESULT( 3 ) = ULPINV
END IF
30 CONTINUE
*
* Do Test (4)
*
DO 50 J = 1, N
TNRM = ONE
IF( WI( J ).EQ.ZERO ) THEN
TNRM = SNRM2( N, VL( 1, J ), 1 )
ELSE IF( WI( J ).GT.ZERO ) THEN
TNRM = SLAPY2( SNRM2( N, VL( 1, J ), 1 ),
$ SNRM2( N, VL( 1, J+1 ), 1 ) )
END IF
RESULT( 4 ) = MAX( RESULT( 4 ),
$ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) )
IF( WI( J ).GT.ZERO ) THEN
VMX = ZERO
VRMX = ZERO
DO 40 JJ = 1, N
VTST = SLAPY2( VL( JJ, J ), VL( JJ, J+1 ) )
IF( VTST.GT.VMX )
$ VMX = VTST
IF( VL( JJ, J+1 ).EQ.ZERO .AND. ABS( VL( JJ, J ) ).GT.
$ VRMX )VRMX = ABS( VL( JJ, J ) )
40 CONTINUE
IF( VRMX / VMX.LT.ONE-TWO*ULP )
$ RESULT( 4 ) = ULPINV
END IF
50 CONTINUE
*
* Test for all options of computing condition numbers
*
DO 200 ISENS = 1, ISENSM
*
SENSE = SENS( ISENS )
*
* Compute eigenvalues only, and test them
*
CALL SLACPY( 'F', N, N, A, LDA, H, LDA )
CALL SGEEVX( BALANC, 'N', 'N', SENSE, N, H, LDA, WR1, WI1, DUM,
$ 1, DUM, 1, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1,
$ RCNDV1, WORK, LWORK, IWORK, IINFO )
IF( IINFO.NE.0 ) THEN
RESULT( 1 ) = ULPINV
IF( JTYPE.NE.22 ) THEN
WRITE( NOUNIT, FMT = 9998 )'SGEEVX2', IINFO, N, JTYPE,
$ BALANC, ISEED
ELSE
WRITE( NOUNIT, FMT = 9999 )'SGEEVX2', IINFO, N,
$ ISEED( 1 )
END IF
INFO = ABS( IINFO )
GO TO 190
END IF
*
* Do Test (5)
*
DO 60 J = 1, N
IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) )
$ RESULT( 5 ) = ULPINV
60 CONTINUE
*
* Do Test (8)
*
IF( .NOT.NOBAL ) THEN
DO 70 J = 1, N
IF( SCALE( J ).NE.SCALE1( J ) )
$ RESULT( 8 ) = ULPINV
70 CONTINUE
IF( ILO.NE.ILO1 )
$ RESULT( 8 ) = ULPINV
IF( IHI.NE.IHI1 )
$ RESULT( 8 ) = ULPINV
IF( ABNRM.NE.ABNRM1 )
$ RESULT( 8 ) = ULPINV
END IF
*
* Do Test (9)
*
IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN
DO 80 J = 1, N
IF( RCONDV( J ).NE.RCNDV1( J ) )
$ RESULT( 9 ) = ULPINV
80 CONTINUE
END IF
*
* Compute eigenvalues and right eigenvectors, and test them
*
CALL SLACPY( 'F', N, N, A, LDA, H, LDA )
CALL SGEEVX( BALANC, 'N', 'V', SENSE, N, H, LDA, WR1, WI1, DUM,
$ 1, LRE, LDLRE, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1,
$ RCNDV1, WORK, LWORK, IWORK, IINFO )
IF( IINFO.NE.0 ) THEN
RESULT( 1 ) = ULPINV
IF( JTYPE.NE.22 ) THEN
WRITE( NOUNIT, FMT = 9998 )'SGEEVX3', IINFO, N, JTYPE,
$ BALANC, ISEED
ELSE
WRITE( NOUNIT, FMT = 9999 )'SGEEVX3', IINFO, N,
$ ISEED( 1 )
END IF
INFO = ABS( IINFO )
GO TO 190
END IF
*
* Do Test (5) again
*
DO 90 J = 1, N
IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) )
$ RESULT( 5 ) = ULPINV
90 CONTINUE
*
* Do Test (6)
*
DO 110 J = 1, N
DO 100 JJ = 1, N
IF( VR( J, JJ ).NE.LRE( J, JJ ) )
$ RESULT( 6 ) = ULPINV
100 CONTINUE
110 CONTINUE
*
* Do Test (8) again
*
IF( .NOT.NOBAL ) THEN
DO 120 J = 1, N
IF( SCALE( J ).NE.SCALE1( J ) )
$ RESULT( 8 ) = ULPINV
120 CONTINUE
IF( ILO.NE.ILO1 )
$ RESULT( 8 ) = ULPINV
IF( IHI.NE.IHI1 )
$ RESULT( 8 ) = ULPINV
IF( ABNRM.NE.ABNRM1 )
$ RESULT( 8 ) = ULPINV
END IF
*
* Do Test (9) again
*
IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN
DO 130 J = 1, N
IF( RCONDV( J ).NE.RCNDV1( J ) )
$ RESULT( 9 ) = ULPINV
130 CONTINUE
END IF
*
* Compute eigenvalues and left eigenvectors, and test them
*
CALL SLACPY( 'F', N, N, A, LDA, H, LDA )
CALL SGEEVX( BALANC, 'V', 'N', SENSE, N, H, LDA, WR1, WI1, LRE,
$ LDLRE, DUM, 1, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1,
$ RCNDV1, WORK, LWORK, IWORK, IINFO )
IF( IINFO.NE.0 ) THEN
RESULT( 1 ) = ULPINV
IF( JTYPE.NE.22 ) THEN
WRITE( NOUNIT, FMT = 9998 )'SGEEVX4', IINFO, N, JTYPE,
$ BALANC, ISEED
ELSE
WRITE( NOUNIT, FMT = 9999 )'SGEEVX4', IINFO, N,
$ ISEED( 1 )
END IF
INFO = ABS( IINFO )
GO TO 190
END IF
*
* Do Test (5) again
*
DO 140 J = 1, N
IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) )
$ RESULT( 5 ) = ULPINV
140 CONTINUE
*
* Do Test (7)
*
DO 160 J = 1, N
DO 150 JJ = 1, N
IF( VL( J, JJ ).NE.LRE( J, JJ ) )
$ RESULT( 7 ) = ULPINV
150 CONTINUE
160 CONTINUE
*
* Do Test (8) again
*
IF( .NOT.NOBAL ) THEN
DO 170 J = 1, N
IF( SCALE( J ).NE.SCALE1( J ) )
$ RESULT( 8 ) = ULPINV
170 CONTINUE
IF( ILO.NE.ILO1 )
$ RESULT( 8 ) = ULPINV
IF( IHI.NE.IHI1 )
$ RESULT( 8 ) = ULPINV
IF( ABNRM.NE.ABNRM1 )
$ RESULT( 8 ) = ULPINV
END IF
*
* Do Test (9) again
*
IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN
DO 180 J = 1, N
IF( RCONDV( J ).NE.RCNDV1( J ) )
$ RESULT( 9 ) = ULPINV
180 CONTINUE
END IF
*
190 CONTINUE
*
200 CONTINUE
*
* If COMP, compare condition numbers to precomputed ones
*
IF( COMP ) THEN
CALL SLACPY( 'F', N, N, A, LDA, H, LDA )
CALL SGEEVX( 'N', 'V', 'V', 'B', N, H, LDA, WR, WI, VL, LDVL,
$ VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV,
$ WORK, LWORK, IWORK, IINFO )
IF( IINFO.NE.0 ) THEN
RESULT( 1 ) = ULPINV
WRITE( NOUNIT, FMT = 9999 )'SGEEVX5', IINFO, N, ISEED( 1 )
INFO = ABS( IINFO )
GO TO 250
END IF
*
* Sort eigenvalues and condition numbers lexicographically
* to compare with inputs
*
DO 220 I = 1, N - 1
KMIN = I
VRMIN = WR( I )
VIMIN = WI( I )
DO 210 J = I + 1, N
IF( WR( J ).LT.VRMIN ) THEN
KMIN = J
VRMIN = WR( J )
VIMIN = WI( J )
END IF
210 CONTINUE
WR( KMIN ) = WR( I )
WI( KMIN ) = WI( I )
WR( I ) = VRMIN
WI( I ) = VIMIN
VRMIN = RCONDE( KMIN )
RCONDE( KMIN ) = RCONDE( I )
RCONDE( I ) = VRMIN
VRMIN = RCONDV( KMIN )
RCONDV( KMIN ) = RCONDV( I )
RCONDV( I ) = VRMIN
220 CONTINUE
*
* Compare condition numbers for eigenvectors
* taking their condition numbers into account
*
RESULT( 10 ) = ZERO
EPS = MAX( EPSIN, ULP )
V = MAX( REAL( N )*EPS*ABNRM, SMLNUM )
IF( ABNRM.EQ.ZERO )
$ V = ONE
DO 230 I = 1, N
IF( V.GT.RCONDV( I )*RCONDE( I ) ) THEN
TOL = RCONDV( I )
ELSE
TOL = V / RCONDE( I )
END IF
IF( V.GT.RCDVIN( I )*RCDEIN( I ) ) THEN
TOLIN = RCDVIN( I )
ELSE
TOLIN = V / RCDEIN( I )
END IF
TOL = MAX( TOL, SMLNUM / EPS )
TOLIN = MAX( TOLIN, SMLNUM / EPS )
IF( EPS*( RCDVIN( I )-TOLIN ).GT.RCONDV( I )+TOL ) THEN
VMAX = ONE / EPS
ELSE IF( RCDVIN( I )-TOLIN.GT.RCONDV( I )+TOL ) THEN
VMAX = ( RCDVIN( I )-TOLIN ) / ( RCONDV( I )+TOL )
ELSE IF( RCDVIN( I )+TOLIN.LT.EPS*( RCONDV( I )-TOL ) ) THEN
VMAX = ONE / EPS
ELSE IF( RCDVIN( I )+TOLIN.LT.RCONDV( I )-TOL ) THEN
VMAX = ( RCONDV( I )-TOL ) / ( RCDVIN( I )+TOLIN )
ELSE
VMAX = ONE
END IF
RESULT( 10 ) = MAX( RESULT( 10 ), VMAX )
230 CONTINUE
*
* Compare condition numbers for eigenvalues
* taking their condition numbers into account
*
RESULT( 11 ) = ZERO
DO 240 I = 1, N
IF( V.GT.RCONDV( I ) ) THEN
TOL = ONE
ELSE
TOL = V / RCONDV( I )
END IF
IF( V.GT.RCDVIN( I ) ) THEN
TOLIN = ONE
ELSE
TOLIN = V / RCDVIN( I )
END IF
TOL = MAX( TOL, SMLNUM / EPS )
TOLIN = MAX( TOLIN, SMLNUM / EPS )
IF( EPS*( RCDEIN( I )-TOLIN ).GT.RCONDE( I )+TOL ) THEN
VMAX = ONE / EPS
ELSE IF( RCDEIN( I )-TOLIN.GT.RCONDE( I )+TOL ) THEN
VMAX = ( RCDEIN( I )-TOLIN ) / ( RCONDE( I )+TOL )
ELSE IF( RCDEIN( I )+TOLIN.LT.EPS*( RCONDE( I )-TOL ) ) THEN
VMAX = ONE / EPS
ELSE IF( RCDEIN( I )+TOLIN.LT.RCONDE( I )-TOL ) THEN
VMAX = ( RCONDE( I )-TOL ) / ( RCDEIN( I )+TOLIN )
ELSE
VMAX = ONE
END IF
RESULT( 11 ) = MAX( RESULT( 11 ), VMAX )
240 CONTINUE
250 CONTINUE
*
END IF
*
9999 FORMAT( ' SGET23: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
$ I6, ', INPUT EXAMPLE NUMBER = ', I4 )
9998 FORMAT( ' SGET23: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
$ I6, ', JTYPE=', I6, ', BALANC = ', A, ', ISEED=(',
$ 3( I5, ',' ), I5, ')' )
*
RETURN
*
* End of SGET23
*
END
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?