zdrvst.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 1,778 行 · 第 1/5 页
F
1,778 行
INDX = 1
DO 390 J = 1, N
DO 380 I = 1, J
WORK( INDX ) = A( I, J )
INDX = INDX + 1
380 CONTINUE
390 CONTINUE
ELSE
INDX = 1
DO 410 J = 1, N
DO 400 I = J, N
WORK( INDX ) = A( I, J )
INDX = INDX + 1
400 CONTINUE
410 CONTINUE
END IF
*
CALL ZHPEVX( 'V', 'I', UPLO, N, WORK, VL, VU, IL, IU,
$ ABSTOL, M2, WA2, Z, LDU, V, RWORK, IWORK,
$ IWORK( 5*N+1 ), IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZHPEVX(V,I,' // 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 460
END IF
END IF
*
* Do tests 19 and 20.
*
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
INDX = 1
DO 430 J = 1, N
DO 420 I = 1, J
WORK( INDX ) = A( I, J )
INDX = INDX + 1
420 CONTINUE
430 CONTINUE
ELSE
INDX = 1
DO 450 J = 1, N
DO 440 I = J, N
WORK( INDX ) = A( I, J )
INDX = INDX + 1
440 CONTINUE
450 CONTINUE
END IF
*
CALL ZHPEVX( 'N', 'I', UPLO, N, WORK, VL, VU, IL, IU,
$ ABSTOL, M3, WA3, Z, LDU, V, RWORK, IWORK,
$ IWORK( 5*N+1 ), IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZHPEVX(N,I,' // UPLO //
$ ')', IINFO, N, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
GO TO 460
END IF
END IF
*
* Do test 21.
*
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 )
*
460 CONTINUE
NTEST = NTEST + 1
IF( IUPLO.EQ.1 ) THEN
INDX = 1
DO 480 J = 1, N
DO 470 I = 1, J
WORK( INDX ) = A( I, J )
INDX = INDX + 1
470 CONTINUE
480 CONTINUE
ELSE
INDX = 1
DO 500 J = 1, N
DO 490 I = J, N
WORK( INDX ) = A( I, J )
INDX = INDX + 1
490 CONTINUE
500 CONTINUE
END IF
*
CALL ZHPEVX( 'V', 'V', UPLO, N, WORK, VL, VU, IL, IU,
$ ABSTOL, M2, WA2, Z, LDU, V, RWORK, IWORK,
$ IWORK( 5*N+1 ), IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZHPEVX(V,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 550
END IF
END IF
*
* Do tests 22 and 23.
*
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
INDX = 1
DO 520 J = 1, N
DO 510 I = 1, J
WORK( INDX ) = A( I, J )
INDX = INDX + 1
510 CONTINUE
520 CONTINUE
ELSE
INDX = 1
DO 540 J = 1, N
DO 530 I = J, N
WORK( INDX ) = A( I, J )
INDX = INDX + 1
530 CONTINUE
540 CONTINUE
END IF
*
CALL ZHPEVX( 'N', 'V', UPLO, N, WORK, VL, VU, IL, IU,
$ ABSTOL, M3, WA3, Z, LDU, V, RWORK, IWORK,
$ IWORK( 5*N+1 ), IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZHPEVX(N,V,' // UPLO //
$ ')', IINFO, N, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
GO TO 550
END IF
END IF
*
IF( M3.EQ.0 .AND. N.GT.0 ) THEN
RESULT( NTEST ) = ULPINV
GO TO 550
END IF
*
* Do test 24.
*
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 )
*
550 CONTINUE
*
* Call ZHBEVD and CHBEVX.
*
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 570 J = 1, N
DO 560 I = MAX( 1, J-KD ), J
V( KD+1+I-J, J ) = A( I, J )
560 CONTINUE
570 CONTINUE
ELSE
DO 590 J = 1, N
DO 580 I = J, MIN( N, J+KD )
V( 1+I-J, J ) = A( I, J )
580 CONTINUE
590 CONTINUE
END IF
*
NTEST = NTEST + 1
CALL ZHBEVD( 'V', UPLO, N, KD, V, LDU, D1, Z, LDU, WORK,
$ LWEDC, RWORK, LRWEDC, IWORK, LIWEDC, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9998 )'ZHBEVD(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 650
END IF
END IF
*
* Do tests 25 and 26.
*
CALL ZHET21( 1, UPLO, N, 0, A, LDA, D1, D2, Z, LDU, V,
$ LDU, TAU, WORK, RWORK, RESULT( NTEST ) )
*
IF( IUPLO.EQ.1 ) THEN
DO 610 J = 1, N
DO 600 I = MAX( 1, J-KD ), J
V( KD+1+I-J, J ) = A( I, J )
600 CONTINUE
610 CONTINUE
ELSE
DO 630 J = 1, N
DO 620 I = J, MIN( N, J+KD )
V( 1+I-J, J ) = A( I, J )
620 CONTINUE
630 CONTINUE
END IF
*
NTEST = NTEST + 2
CALL ZHBEVD( 'N', UPLO, N, KD, V, LDU, D3, Z, LDU, WORK,
$ LWEDC, RWORK, LRWEDC, IWORK, LIWEDC, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9998 )'ZHBEVD(N,' // UPLO //
$ ')', IINFO, N, KD, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
GO TO 650
END IF
END IF
*
* Do test 27.
*
TEMP1 = ZERO
TEMP2 = ZERO
DO 640 J = 1, N
TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D3( J ) ) )
TEMP2 = MAX( TEMP2, ABS( D1( J )-D3( J ) ) )
640 CONTINUE
RESULT( NTEST ) = TEMP2 / MAX( UNFL,
$ ULP*MAX( TEMP1, TEMP2 ) )
*
* Load array V with the upper or lower triangular part
* of the matrix in band form.
*
650 CONTINUE
IF( IUPLO.EQ.1 ) THEN
DO 670 J = 1, N
DO 660 I = MAX( 1, J-KD ), J
V( KD+1+I-J, J ) = A( I, J )
660 CONTINUE
670 CONTINUE
ELSE
DO 690 J = 1, N
DO 680 I = J, MIN( N, J+KD )
V( 1+I-J, J ) = A( I, J )
680 CONTINUE
690 CONTINUE
END IF
*
NTEST = NTEST + 1
CALL ZHBEVX( 'V', 'A', UPLO, N, KD, V, LDU, U, LDU, VL,
$ VU, IL, IU, ABSTOL, M, WA1, Z, LDU, WORK,
$ RWORK, IWORK, IWORK( 5*N+1 ), IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZHBEVX(V,A,' // 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 750
END IF
END IF
*
* Do tests 28 and 29.
*
CALL ZHET21( 1, UPLO, N, 0, A, LDU, WA1, D2, Z, LDU, V,
$ LDU, TAU, WORK, RWORK, RESULT( NTEST ) )
*
NTEST = NTEST + 2
*
IF( IUPLO.EQ.1 ) THEN
DO 710 J = 1, N
DO 700 I = MAX( 1, J-KD ), J
V( KD+1+I-J, J ) = A( I, J )
700 CONTINUE
710 CONTINUE
ELSE
DO 730 J = 1, N
DO 720 I = J, MIN( N, J+KD )
V( 1+I-J, J ) = A( I, J )
720 CONTINUE
730 CONTINUE
END IF
*
CALL ZHBEVX( 'N', 'A', 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(N,A,' // UPLO //
$ ')', IINFO, N, KD, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
GO TO 750
END IF
END IF
*
* Do test 30.
*
TEMP1 = ZERO
TEMP2 = ZERO
DO 740 J = 1, N
TEMP1 = MAX( TEMP1, ABS( WA1( J ) ), ABS( WA2( J ) ) )
TEMP2 = MAX( TEMP2, ABS( WA1( J )-WA2( J ) ) )
740 CONTINUE
RESULT( NTEST ) = TEMP2 / MAX( UNFL,
$ ULP*MAX( TEMP1, TEMP2 ) )
*
* Load array V with the upper or lower triangular part
* of the matrix in band form.
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?