zdrvst.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 1,778 行 · 第 1/5 页
F
1,778 行
*
750 CONTINUE
NTEST = NTEST + 1
IF( IUPLO.EQ.1 ) THEN
DO 770 J = 1, N
DO 760 I = MAX( 1, J-KD ), J
V( KD+1+I-J, J ) = A( I, J )
760 CONTINUE
770 CONTINUE
ELSE
DO 790 J = 1, N
DO 780 I = J, MIN( N, J+KD )
V( 1+I-J, J ) = A( I, J )
780 CONTINUE
790 CONTINUE
END IF
*
CALL ZHBEVX( 'V', 'I', UPLO, N, KD, V, LDU, U, LDU, VL,
$ VU, IL, IU, ABSTOL, M2, WA2, Z, LDU, WORK,
$ RWORK, IWORK, IWORK( 5*N+1 ), IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9998 )'ZHBEVX(V,I,' // UPLO //
$ ')', IINFO, N, KD, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
RESULT( NTEST+1 ) = ULPINV
RESULT( NTEST+2 ) = ULPINV
GO TO 840
END IF
END IF
*
* Do tests 31 and 32.
*
CALL ZHET22( 1, UPLO, N, M2, 0, A, LDU, WA2, D2, Z, LDU,
$ V, LDU, TAU, WORK, RWORK, RESULT( NTEST ) )
*
NTEST = NTEST + 2
*
IF( IUPLO.EQ.1 ) THEN
DO 810 J = 1, N
DO 800 I = MAX( 1, J-KD ), J
V( KD+1+I-J, J ) = A( I, J )
800 CONTINUE
810 CONTINUE
ELSE
DO 830 J = 1, N
DO 820 I = J, MIN( N, J+KD )
V( 1+I-J, J ) = A( I, J )
820 CONTINUE
830 CONTINUE
END IF
CALL ZHBEVX( 'N', 'I', UPLO, N, KD, V, LDU, U, LDU, VL,
$ VU, IL, IU, ABSTOL, M3, WA3, Z, LDU, WORK,
$ RWORK, IWORK, IWORK( 5*N+1 ), IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9998 )'ZHBEVX(N,I,' // UPLO //
$ ')', IINFO, N, KD, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
GO TO 840
END IF
END IF
*
* Do test 33.
*
TEMP1 = DSXT1( 1, WA2, M2, WA3, M3, ABSTOL, ULP, UNFL )
TEMP2 = DSXT1( 1, WA3, M3, WA2, M2, ABSTOL, ULP, UNFL )
IF( N.GT.0 ) THEN
TEMP3 = MAX( ABS( WA1( 1 ) ), ABS( WA1( N ) ) )
ELSE
TEMP3 = ZERO
END IF
RESULT( NTEST ) = ( TEMP1+TEMP2 ) /
$ MAX( UNFL, TEMP3*ULP )
*
* Load array V with the upper or lower triangular part
* of the matrix in band form.
*
840 CONTINUE
NTEST = NTEST + 1
IF( IUPLO.EQ.1 ) THEN
DO 860 J = 1, N
DO 850 I = MAX( 1, J-KD ), J
V( KD+1+I-J, J ) = A( I, J )
850 CONTINUE
860 CONTINUE
ELSE
DO 880 J = 1, N
DO 870 I = J, MIN( N, J+KD )
V( 1+I-J, J ) = A( I, J )
870 CONTINUE
880 CONTINUE
END IF
CALL ZHBEVX( 'V', 'V', UPLO, N, KD, V, LDU, U, LDU, VL,
$ VU, IL, IU, ABSTOL, M2, WA2, Z, LDU, WORK,
$ RWORK, IWORK, IWORK( 5*N+1 ), IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9998 )'ZHBEVX(V,V,' // UPLO //
$ ')', IINFO, N, KD, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
RESULT( NTEST+1 ) = ULPINV
RESULT( NTEST+2 ) = ULPINV
GO TO 930
END IF
END IF
*
* Do tests 34 and 35.
*
CALL ZHET22( 1, UPLO, N, M2, 0, A, LDU, WA2, D2, Z, LDU,
$ V, LDU, TAU, WORK, RWORK, RESULT( NTEST ) )
*
NTEST = NTEST + 2
*
IF( IUPLO.EQ.1 ) THEN
DO 900 J = 1, N
DO 890 I = MAX( 1, J-KD ), J
V( KD+1+I-J, J ) = A( I, J )
890 CONTINUE
900 CONTINUE
ELSE
DO 920 J = 1, N
DO 910 I = J, MIN( N, J+KD )
V( 1+I-J, J ) = A( I, J )
910 CONTINUE
920 CONTINUE
END IF
CALL ZHBEVX( 'N', 'V', UPLO, N, KD, V, LDU, U, LDU, VL,
$ VU, IL, IU, ABSTOL, M3, WA3, Z, LDU, WORK,
$ RWORK, IWORK, IWORK( 5*N+1 ), IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9998 )'ZHBEVX(N,V,' // UPLO //
$ ')', IINFO, N, KD, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
GO TO 930
END IF
END IF
*
IF( M3.EQ.0 .AND. N.GT.0 ) THEN
RESULT( NTEST ) = ULPINV
GO TO 930
END IF
*
* Do test 36.
*
TEMP1 = DSXT1( 1, WA2, M2, WA3, M3, ABSTOL, ULP, UNFL )
TEMP2 = DSXT1( 1, WA3, M3, WA2, M2, ABSTOL, ULP, UNFL )
IF( N.GT.0 ) THEN
TEMP3 = MAX( ABS( WA1( 1 ) ), ABS( WA1( N ) ) )
ELSE
TEMP3 = ZERO
END IF
RESULT( NTEST ) = ( TEMP1+TEMP2 ) /
$ MAX( UNFL, TEMP3*ULP )
*
930 CONTINUE
*
* Call ZHEEV
*
CALL ZLACPY( ' ', N, N, A, LDA, V, LDU )
*
NTEST = NTEST + 1
CALL ZHEEV( 'V', UPLO, N, A, LDU, D1, WORK, LWORK, RWORK,
$ IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZHEEV(V,' // UPLO // ')',
$ IINFO, N, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
RESULT( NTEST+1 ) = ULPINV
RESULT( NTEST+2 ) = ULPINV
GO TO 950
END IF
END IF
*
* Do tests 37 and 38
*
CALL ZHET21( 1, UPLO, N, 0, V, LDU, D1, D2, A, LDU, Z,
$ LDU, TAU, WORK, RWORK, RESULT( NTEST ) )
*
CALL ZLACPY( ' ', N, N, V, LDU, A, LDA )
*
NTEST = NTEST + 2
CALL ZHEEV( 'N', UPLO, N, A, LDU, D3, WORK, LWORK, RWORK,
$ IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZHEEV(N,' // UPLO // ')',
$ IINFO, N, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
GO TO 950
END IF
END IF
*
* Do test 39
*
TEMP1 = ZERO
TEMP2 = ZERO
DO 940 J = 1, N
TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D3( J ) ) )
TEMP2 = MAX( TEMP2, ABS( D1( J )-D3( J ) ) )
940 CONTINUE
RESULT( NTEST ) = TEMP2 / MAX( UNFL,
$ ULP*MAX( TEMP1, TEMP2 ) )
*
950 CONTINUE
*
CALL ZLACPY( ' ', N, N, V, LDU, A, LDA )
*
* Call ZHPEV
*
* Load array WORK with the upper or lower triangular
* part of the matrix in packed form.
*
IF( IUPLO.EQ.1 ) THEN
INDX = 1
DO 970 J = 1, N
DO 960 I = 1, J
WORK( INDX ) = A( I, J )
INDX = INDX + 1
960 CONTINUE
970 CONTINUE
ELSE
INDX = 1
DO 990 J = 1, N
DO 980 I = J, N
WORK( INDX ) = A( I, J )
INDX = INDX + 1
980 CONTINUE
990 CONTINUE
END IF
*
NTEST = NTEST + 1
INDWRK = N*( N+1 ) / 2 + 1
CALL ZHPEV( 'V', UPLO, N, WORK, D1, Z, LDU,
$ WORK( INDWRK ), RWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZHPEV(V,' // UPLO // ')',
$ IINFO, N, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
RESULT( NTEST+1 ) = ULPINV
RESULT( NTEST+2 ) = ULPINV
GO TO 1050
END IF
END IF
*
* Do tests 40 and 41.
*
CALL ZHET21( 1, UPLO, N, 0, A, LDA, D1, D2, Z, LDU, V,
$ LDU, TAU, WORK, RWORK, RESULT( NTEST ) )
*
IF( IUPLO.EQ.1 ) THEN
INDX = 1
DO 1010 J = 1, N
DO 1000 I = 1, J
WORK( INDX ) = A( I, J )
INDX = INDX + 1
1000 CONTINUE
1010 CONTINUE
ELSE
INDX = 1
DO 1030 J = 1, N
DO 1020 I = J, N
WORK( INDX ) = A( I, J )
INDX = INDX + 1
1020 CONTINUE
1030 CONTINUE
END IF
*
NTEST = NTEST + 2
INDWRK = N*( N+1 ) / 2 + 1
CALL ZHPEV( 'N', UPLO, N, WORK, D3, Z, LDU,
$ WORK( INDWRK ), RWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZHPEV(N,' // UPLO // ')',
$ IINFO, N, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
GO TO 1050
END IF
END IF
*
* Do test 42
*
TEMP1 = ZERO
TEMP2 = ZERO
DO 1040 J = 1, N
TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D3( J ) ) )
TEMP2 = MAX( TEMP2, ABS( D1( J )-D3( J ) ) )
1040 CONTINUE
RESULT( NTEST ) = TEMP2 / MAX( UNFL,
$ ULP*MAX( TEMP1, TEMP2 ) )
*
1050 CONTINUE
*
* Call ZHBEV
*
IF( JTYPE.LE.7 ) THEN
KD = 0
ELSE IF( JTYPE.GE.8 .AND. JTYPE.LE.15 ) THEN
KD = MAX( N-1, 0 )
ELSE
KD = IHBW
END IF
*
* Load array V with the upper or lower triangular part
* of the matrix in band form.
*
IF( IUPLO.EQ.1 ) THEN
DO 1070 J = 1, N
DO 1060 I = MAX( 1, J-KD ), J
V( KD+1+I-J, J ) = A( I, J )
1060 CONTINUE
1070 CONTINUE
ELSE
DO 1090 J = 1, N
DO 1080 I = J, MIN( N, J+KD )
V( 1+I-J, J ) = A( I, J )
1080 CONTINUE
1090 CONTINUE
END IF
*
NTEST = NTEST + 1
CALL ZHBEV( 'V', UPLO, N, KD, V, LDU, D1, Z, LDU, WORK,
$ RWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9998 )'ZHBEV(V,' // UPLO // ')',
$ ⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?