zdrvst.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 1,778 行 · 第 1/5 页
F
1,778 行
DATA KTYPE / 1, 2, 5*4, 5*5, 3*8, 3*9 /
DATA KMAGN / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
$ 2, 3, 1, 2, 3 /
DATA KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
$ 0, 0, 4, 4, 4 /
* ..
* .. Executable Statements ..
*
* 1) Check for errors
*
NTESTT = 0
INFO = 0
*
BADNN = .FALSE.
NMAX = 1
DO 10 J = 1, NSIZES
NMAX = MAX( NMAX, NN( J ) )
IF( NN( J ).LT.0 )
$ BADNN = .TRUE.
10 CONTINUE
*
* Check for errors
*
IF( NSIZES.LT.0 ) THEN
INFO = -1
ELSE IF( BADNN ) THEN
INFO = -2
ELSE IF( NTYPES.LT.0 ) THEN
INFO = -3
ELSE IF( LDA.LT.NMAX ) THEN
INFO = -9
ELSE IF( LDU.LT.NMAX ) THEN
INFO = -16
ELSE IF( 2*MAX( 2, NMAX )**2.GT.LWORK ) THEN
INFO = -22
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'ZDRVST', -INFO )
RETURN
END IF
*
* Quick return if nothing to do
*
IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
$ RETURN
*
* More Important constants
*
UNFL = DLAMCH( 'Safe minimum' )
OVFL = DLAMCH( 'Overflow' )
CALL DLABAD( UNFL, OVFL )
ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
ULPINV = ONE / ULP
RTUNFL = SQRT( UNFL )
RTOVFL = SQRT( OVFL )
*
* Loop over sizes, types
*
DO 20 I = 1, 4
ISEED2( I ) = ISEED( I )
ISEED3( I ) = ISEED( I )
20 CONTINUE
*
NERRS = 0
NMATS = 0
*
DO 1220 JSIZE = 1, NSIZES
N = NN( JSIZE )
IF( N.GT.0 ) THEN
LGN = INT( LOG( DBLE( N ) ) / LOG( TWO ) )
IF( 2**LGN.LT.N )
$ LGN = LGN + 1
IF( 2**LGN.LT.N )
$ LGN = LGN + 1
LWEDC = MAX( 2*N+N*N, 2*N*N )
LRWEDC = 1 + 4*N + 2*N*LGN + 3*N**2
LIWEDC = 3 + 5*N
ELSE
LWEDC = 2
LRWEDC = 8
LIWEDC = 8
END IF
ANINV = ONE / DBLE( MAX( 1, N ) )
*
IF( NSIZES.NE.1 ) THEN
MTYPES = MIN( MAXTYP, NTYPES )
ELSE
MTYPES = MIN( MAXTYP+1, NTYPES )
END IF
*
DO 1210 JTYPE = 1, MTYPES
IF( .NOT.DOTYPE( JTYPE ) )
$ GO TO 1210
NMATS = NMATS + 1
NTEST = 0
*
DO 30 J = 1, 4
IOLDSD( J ) = ISEED( J )
30 CONTINUE
*
* 2) Compute "A"
*
* Control parameters:
*
* KMAGN KMODE KTYPE
* =1 O(1) clustered 1 zero
* =2 large clustered 2 identity
* =3 small exponential (none)
* =4 arithmetic diagonal, (w/ eigenvalues)
* =5 random log Hermitian, w/ eigenvalues
* =6 random (none)
* =7 random diagonal
* =8 random Hermitian
* =9 band Hermitian, w/ eigenvalues
*
IF( MTYPES.GT.MAXTYP )
$ GO TO 110
*
ITYPE = KTYPE( JTYPE )
IMODE = KMODE( JTYPE )
*
* Compute norm
*
GO TO ( 40, 50, 60 )KMAGN( JTYPE )
*
40 CONTINUE
ANORM = ONE
GO TO 70
*
50 CONTINUE
ANORM = ( RTOVFL*ULP )*ANINV
GO TO 70
*
60 CONTINUE
ANORM = RTUNFL*N*ULPINV
GO TO 70
*
70 CONTINUE
*
CALL ZLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
IINFO = 0
COND = ULPINV
*
* Special Matrices -- Identity & Jordan block
*
* Zero
*
IF( ITYPE.EQ.1 ) THEN
IINFO = 0
*
ELSE IF( ITYPE.EQ.2 ) THEN
*
* Identity
*
DO 80 JCOL = 1, N
A( JCOL, JCOL ) = ANORM
80 CONTINUE
*
ELSE IF( ITYPE.EQ.4 ) THEN
*
* Diagonal Matrix, [Eigen]values Specified
*
CALL ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
$ ANORM, 0, 0, 'N', A, LDA, WORK, IINFO )
*
ELSE IF( ITYPE.EQ.5 ) THEN
*
* Hermitian, eigenvalues specified
*
CALL ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
$ ANORM, N, N, 'N', A, LDA, WORK, IINFO )
*
ELSE IF( ITYPE.EQ.7 ) THEN
*
* Diagonal, random eigenvalues
*
CALL ZLATMR( N, N, 'S', ISEED, 'H', WORK, 6, ONE, CONE,
$ 'T', 'N', WORK( N+1 ), 1, ONE,
$ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
$ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
*
ELSE IF( ITYPE.EQ.8 ) THEN
*
* Hermitian, random eigenvalues
*
CALL ZLATMR( N, N, 'S', ISEED, 'H', WORK, 6, ONE, CONE,
$ 'T', 'N', WORK( N+1 ), 1, ONE,
$ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
$ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
*
ELSE IF( ITYPE.EQ.9 ) THEN
*
* Hermitian banded, eigenvalues specified
*
IHBW = INT( ( N-1 )*DLARND( 1, ISEED3 ) )
CALL ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
$ ANORM, IHBW, IHBW, 'Z', U, LDU, WORK,
$ IINFO )
*
* Store as dense matrix for most routines.
*
CALL ZLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
DO 100 IDIAG = -IHBW, IHBW
IROW = IHBW - IDIAG + 1
J1 = MAX( 1, IDIAG+1 )
J2 = MIN( N, N+IDIAG )
DO 90 J = J1, J2
I = J - IDIAG
A( I, J ) = U( IROW, J )
90 CONTINUE
100 CONTINUE
ELSE
IINFO = 1
END IF
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
RETURN
END IF
*
110 CONTINUE
*
ABSTOL = UNFL + UNFL
IF( N.LE.1 ) THEN
IL = 1
IU = N
ELSE
IL = 1 + INT( ( N-1 )*DLARND( 1, ISEED2 ) )
IU = 1 + INT( ( N-1 )*DLARND( 1, ISEED2 ) )
IF( IL.GT.IU ) THEN
ITEMP = IL
IL = IU
IU = ITEMP
END IF
END IF
*
* Perform tests storing upper or lower triangular
* part of matrix.
*
DO 1200 IUPLO = 0, 1
IF( IUPLO.EQ.0 ) THEN
UPLO = 'L'
ELSE
UPLO = 'U'
END IF
*
* Call ZHEEVD and CHEEVX.
*
CALL ZLACPY( ' ', N, N, A, LDA, V, LDU )
*
NTEST = NTEST + 1
CALL ZHEEVD( 'V', UPLO, N, A, LDU, D1, WORK, LWEDC,
$ RWORK, LRWEDC, IWORK, LIWEDC, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZHEEVD(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 130
END IF
END IF
*
* Do tests 1 and 2.
*
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 ZHEEVD( 'N', UPLO, N, A, LDU, D3, WORK, LWEDC,
$ RWORK, LRWEDC, IWORK, LIWEDC, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZHEEVD(N,' // UPLO //
$ ')', IINFO, N, JTYPE, IOLDSD
INFO = ABS( IINFO )
IF( IINFO.LT.0 ) THEN
RETURN
ELSE
RESULT( NTEST ) = ULPINV
GO TO 130
END IF
END IF
*
* Do test 3.
*
TEMP1 = ZERO
TEMP2 = ZERO
DO 120 J = 1, N
TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D3( J ) ) )
TEMP2 = MAX( TEMP2, ABS( D1( J )-D3( J ) ) )
120 CONTINUE
RESULT( NTEST ) = TEMP2 / MAX( UNFL,
$ ULP*MAX( TEMP1, TEMP2 ) )
*
130 CONTINUE
CALL ZLACPY( ' ', N, N, V, LDU, A, LDA )
*
NTEST = NTEST + 1
*
IF( N.GT.0 ) THEN
TEMP3 = MAX( ABS( D1( 1 ) ), ABS( D1( N ) ) )
IF( IL.NE.1 ) THEN
VL = D1( IL ) - MAX( HALF*( D1( IL )-D1( IL-1 ) ),
$ TEN*ULP*TEMP3, TEN*RTUNFL )
ELSE IF( N.GT.0 ) THEN
VL = D1( 1 ) - MAX( HALF*( D1( N )-D1( 1 ) ),
$ TEN*ULP*TEMP3, TEN*RTUNFL )
END IF
IF( IU.NE.N ) THEN
VU = D1( IU ) + MAX( HALF*( D1( IU+1 )-D1( IU ) ),
$ TEN*ULP*TEMP3, TEN*RTUNFL )
ELSE IF( N.GT.0 ) THEN
VU = D1( N ) + MAX( HALF*( D1( N )-D1( 1 ) ),
$ TEN*ULP*TEMP3, TEN*RTUNFL )
END IF
ELSE
TEMP3 = ZERO
VL = ZERO
VU = ONE
END IF
*
CALL ZHEEVX( 'V', 'A', UPLO, N, A, LDU, VL, VU, IL, IU,
$ ABSTOL, M, WA1, Z, LDU, WORK, LWORK, RWORK,
$ IWORK, IWORK( 5*N+1 ), IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZHEEVX(V,A,' // 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 150
END IF
END IF
*
* Do tests 4 and 5.
*
CALL ZLACPY( ' ', N, N, V, LDU, A, LDA )
*
CALL ZHET21( 1, UPLO, N, 0, A, LDU, WA1, D2, Z, LDU, V,
$ LDU, TAU, WORK, RWORK, RESULT( NTEST ) )
*
NTEST = NTEST + 2
CALL ZHEEVX( 'N', 'A', UPLO, N, A, LDU, VL, VU, IL, IU,
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?