schkst.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 1,804 行 · 第 1/5 页
F
1,804 行
CALL SSYTRD( 'U', N, V, LDU, SD, SE, TAU, WORK, LWORK,
$ IINFO )
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SSYTRD(U)', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 1 ) = ULPINV
GO TO 280
END IF
END IF
*
CALL SLACPY( 'U', N, N, V, LDU, U, LDU )
*
NTEST = 2
CALL SORGTR( 'U', N, U, LDU, TAU, WORK, LWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SORGTR(U)', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 2 ) = ULPINV
GO TO 280
END IF
END IF
*
* Do tests 1 and 2
*
CALL SSYT21( 2, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V,
$ LDU, TAU, WORK, RESULT( 1 ) )
CALL SSYT21( 3, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V,
$ LDU, TAU, WORK, RESULT( 2 ) )
*
* Call SSYTRD and SORGTR to compute S and U from
* lower triangle, do tests.
*
CALL SLACPY( 'L', N, N, A, LDA, V, LDU )
*
NTEST = 3
CALL SSYTRD( 'L', N, V, LDU, SD, SE, TAU, WORK, LWORK,
$ IINFO )
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SSYTRD(L)', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 3 ) = ULPINV
GO TO 280
END IF
END IF
*
CALL SLACPY( 'L', N, N, V, LDU, U, LDU )
*
NTEST = 4
CALL SORGTR( 'L', N, U, LDU, TAU, WORK, LWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SORGTR(L)', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 4 ) = ULPINV
GO TO 280
END IF
END IF
*
CALL SSYT21( 2, 'Lower', N, 1, A, LDA, SD, SE, U, LDU, V,
$ LDU, TAU, WORK, RESULT( 3 ) )
CALL SSYT21( 3, 'Lower', N, 1, A, LDA, SD, SE, U, LDU, V,
$ LDU, TAU, WORK, RESULT( 4 ) )
*
* Store the upper triangle of A in AP
*
I = 0
DO 120 JC = 1, N
DO 110 JR = 1, JC
I = I + 1
AP( I ) = A( JR, JC )
110 CONTINUE
120 CONTINUE
*
* Call SSPTRD and SOPGTR to compute S and U from AP
*
CALL SCOPY( NAP, AP, 1, VP, 1 )
*
NTEST = 5
CALL SSPTRD( 'U', N, VP, SD, SE, TAU, IINFO )
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SSPTRD(U)', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 5 ) = ULPINV
GO TO 280
END IF
END IF
*
NTEST = 6
CALL SOPGTR( 'U', N, VP, TAU, U, LDU, WORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SOPGTR(U)', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 6 ) = ULPINV
GO TO 280
END IF
END IF
*
* Do tests 5 and 6
*
CALL SSPT21( 2, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU,
$ WORK, RESULT( 5 ) )
CALL SSPT21( 3, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU,
$ WORK, RESULT( 6 ) )
*
* Store the lower triangle of A in AP
*
I = 0
DO 140 JC = 1, N
DO 130 JR = JC, N
I = I + 1
AP( I ) = A( JR, JC )
130 CONTINUE
140 CONTINUE
*
* Call SSPTRD and SOPGTR to compute S and U from AP
*
CALL SCOPY( NAP, AP, 1, VP, 1 )
*
NTEST = 7
CALL SSPTRD( 'L', N, VP, SD, SE, TAU, IINFO )
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SSPTRD(L)', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 7 ) = ULPINV
GO TO 280
END IF
END IF
*
NTEST = 8
CALL SOPGTR( 'L', N, VP, TAU, U, LDU, WORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SOPGTR(L)', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 8 ) = ULPINV
GO TO 280
END IF
END IF
*
CALL SSPT21( 2, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU,
$ WORK, RESULT( 7 ) )
CALL SSPT21( 3, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU,
$ WORK, RESULT( 8 ) )
*
* Call SSTEQR to compute D1, D2, and Z, do tests.
*
* Compute D1 and Z
*
CALL SCOPY( N, SD, 1, D1, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, WORK, 1 )
CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
*
NTEST = 9
CALL SSTEQR( 'V', N, D1, WORK, Z, LDU, WORK( N+1 ), IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SSTEQR(V)', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 9 ) = ULPINV
GO TO 280
END IF
END IF
*
* Compute D2
*
CALL SCOPY( N, SD, 1, D2, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, WORK, 1 )
*
NTEST = 11
CALL SSTEQR( 'N', N, D2, WORK, WORK( N+1 ), LDU,
$ WORK( N+1 ), IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SSTEQR(N)', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 11 ) = ULPINV
GO TO 280
END IF
END IF
*
* Compute D3 (using PWK method)
*
CALL SCOPY( N, SD, 1, D3, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, WORK, 1 )
*
NTEST = 12
CALL SSTERF( N, D3, WORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SSTERF', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 12 ) = ULPINV
GO TO 280
END IF
END IF
*
* Do Tests 9 and 10
*
CALL SSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
$ RESULT( 9 ) )
*
* Do Tests 11 and 12
*
TEMP1 = ZERO
TEMP2 = ZERO
TEMP3 = ZERO
TEMP4 = ZERO
*
DO 150 J = 1, N
TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
TEMP3 = MAX( TEMP3, ABS( D1( J ) ), ABS( D3( J ) ) )
TEMP4 = MAX( TEMP4, ABS( D1( J )-D3( J ) ) )
150 CONTINUE
*
RESULT( 11 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
RESULT( 12 ) = TEMP4 / MAX( UNFL, ULP*MAX( TEMP3, TEMP4 ) )
*
* Do Test 13 -- Sturm Sequence Test of Eigenvalues
* Go up by factors of two until it succeeds
*
NTEST = 13
TEMP1 = THRESH*( HALF-ULP )
*
DO 160 J = 0, LOG2UI
CALL SSTECH( N, SD, SE, D1, TEMP1, WORK, IINFO )
IF( IINFO.EQ.0 )
$ GO TO 170
TEMP1 = TEMP1*TWO
160 CONTINUE
*
170 CONTINUE
RESULT( 13 ) = TEMP1
*
* For positive definite matrices ( JTYPE.GT.15 ) call SPTEQR
* and do tests 14, 15, and 16 .
*
IF( JTYPE.GT.15 ) THEN
*
* Compute D4 and Z4
*
CALL SCOPY( N, SD, 1, D4, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, WORK, 1 )
CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
*
NTEST = 14
CALL SPTEQR( 'V', N, D4, WORK, Z, LDU, WORK( N+1 ),
$ IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SPTEQR(V)', IINFO, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 14 ) = ULPINV
GO TO 280
END IF
END IF
*
* Do Tests 14 and 15
*
CALL SSTT21( N, 0, SD, SE, D4, DUMMA, Z, LDU, WORK,
$ RESULT( 14 ) )
*
* Compute D5
*
CALL SCOPY( N, SD, 1, D5, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, WORK, 1 )
*
NTEST = 16
CALL SPTEQR( 'N', N, D5, WORK, Z, LDU, WORK( N+1 ),
$ IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SPTEQR(N)', IINFO, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( 16 ) = ULPINV
GO TO 280
END IF
END IF
*
* Do Test 16
*
TEMP1 = ZERO
TEMP2 = ZERO
DO 180 J = 1, N
TEMP1 = MAX( TEMP1, ABS( D4( J ) ), ABS( D5( J ) ) )
TEMP2 = MAX( TEMP2, ABS( D4( J )-D5( J ) ) )
180 CONTINUE
*
RESULT( 16 ) = TEMP2 / MAX( UNFL,
$ HUN*ULP*MAX( TEMP1, TEMP2 ) )
ELSE
RESULT( 14 ) = ZERO
RESULT( 15 ) = ZERO
RESULT( 16 ) = ZERO
END IF
*
* Call SSTEBZ with different options and do tests 17-18.
*
* If S is positive definite and diagonally dominant,
* ask for all eigenvalues with high relative accuracy.
*
VL = ZERO
VU = ZERO
IL = 0
IU = 0
IF( JTYPE.EQ.21 ) THEN
NTEST = 17
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?