dchkst.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 1,801 行 · 第 1/5 页
F
1,801 行
* Z (workspace/output) DOUBLE PRECISION array of
* dimension( LDU, max(NN) ).
* The orthogonal matrix of eigenvectors computed by DSTEQR,
* DPTEQR, and DSTEIN.
*
* WORK (workspace/output) DOUBLE PRECISION array of
* dimension( LWORK )
*
* LWORK (input) INTEGER
* The number of entries in WORK. This must be at least
* 1 + 4 * Nmax + 2 * Nmax * lg Nmax + 3 * Nmax**2
* where Nmax = max( NN(j), 2 ) and lg = log base 2.
*
* IWORK (workspace/output) INTEGER array,
* dimension (6 + 6*Nmax + 5 * Nmax * lg Nmax )
* where Nmax = max( NN(j), 2 ) and lg = log base 2.
* Workspace.
*
* RESULT (output) DOUBLE PRECISION array, dimension (26)
* The values computed by the tests described above.
* The values are currently limited to 1/ulp, to avoid
* overflow.
*
* INFO (output) INTEGER
* If 0, then everything ran OK.
* -1: NSIZES < 0
* -2: Some NN(j) < 0
* -3: NTYPES < 0
* -5: THRESH < 0
* -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
* -23: LDU < 1 or LDU < NMAX.
* -29: LWORK too small.
* If DLATMR, SLATMS, DSYTRD, DORGTR, DSTEQR, SSTERF,
* or DORMC2 returns an error code, the
* absolute value of it is returned.
*
*-----------------------------------------------------------------------
*
* Some Local Variables and Parameters:
* ---- ----- --------- --- ----------
* ZERO, ONE Real 0 and 1.
* MAXTYP The number of types defined.
* NTEST The number of tests performed, or which can
* be performed so far, for the current matrix.
* NTESTT The total number of tests performed so far.
* NBLOCK Blocksize as returned by ENVIR.
* NMAX Largest value in NN.
* NMATS The number of matrices generated so far.
* NERRS The number of tests which have exceeded THRESH
* so far.
* COND, IMODE Values to be passed to the matrix generators.
* ANORM Norm of A; passed to matrix generators.
*
* OVFL, UNFL Overflow and underflow thresholds.
* ULP, ULPINV Finest relative precision and its inverse.
* RTOVFL, RTUNFL Square roots of the previous 2 values.
* The following four arrays decode JTYPE:
* KTYPE(j) The general type (1-10) for type "j".
* KMODE(j) The MODE value to be passed to the matrix
* generator for type "j".
* KMAGN(j) The order of magnitude ( O(1),
* O(overflow^(1/2) ), O(underflow^(1/2) )
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE, TWO, EIGHT, TEN, HUN
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
$ EIGHT = 8.0D0, TEN = 10.0D0, HUN = 100.0D0 )
DOUBLE PRECISION HALF
PARAMETER ( HALF = ONE / TWO )
INTEGER MAXTYP
PARAMETER ( MAXTYP = 21 )
LOGICAL SRANGE
PARAMETER ( SRANGE = .FALSE. )
LOGICAL SREL
PARAMETER ( SREL = .FALSE. )
* ..
* .. Local Scalars ..
LOGICAL BADNN, TRYRAC
INTEGER I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, J, JC,
$ JR, JSIZE, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC,
$ M, M2, M3, MTYPES, N, NAP, NBLOCK, NERRS,
$ NMATS, NMAX, NSPLIT, NTEST, NTESTT
DOUBLE PRECISION ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL,
$ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP,
$ ULPINV, UNFL, VL, VU
* ..
* .. Local Arrays ..
INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
$ KMAGN( MAXTYP ), KMODE( MAXTYP ),
$ KTYPE( MAXTYP )
DOUBLE PRECISION DUMMA( 1 )
* ..
* .. External Functions ..
INTEGER ILAENV
DOUBLE PRECISION DLAMCH, DLARND, DSXT1
EXTERNAL ILAENV, DLAMCH, DLARND, DSXT1
* ..
* .. External Subroutines ..
EXTERNAL DCOPY, DLABAD, DLACPY, DLASET, DLASUM, DLATMR,
$ DLATMS, DOPGTR, DORGTR, DPTEQR, DSPT21, DSPTRD,
$ DSTEBZ, DSTECH, DSTEDC, DSTEMR, DSTEIN, DSTEQR,
$ DSTERF, DSTT21, DSTT22, DSYT21, DSYTRD, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, INT, LOG, MAX, MIN, SQRT
* ..
* .. Data statements ..
DATA KTYPE / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8,
$ 8, 8, 9, 9, 9, 9, 9, 10 /
DATA KMAGN / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
$ 2, 3, 1, 1, 1, 2, 3, 1 /
DATA KMODE / 0, 0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
$ 0, 0, 4, 3, 1, 4, 4, 3 /
* ..
* .. Executable Statements ..
*
* Keep ftnchek happy
IDUMMA( 1 ) = 1
*
* Check for errors
*
NTESTT = 0
INFO = 0
*
* Important constants
*
BADNN = .FALSE.
TRYRAC = .TRUE.
NMAX = 1
DO 10 J = 1, NSIZES
NMAX = MAX( NMAX, NN( J ) )
IF( NN( J ).LT.0 )
$ BADNN = .TRUE.
10 CONTINUE
*
NBLOCK = ILAENV( 1, 'DSYTRD', 'L', NMAX, -1, -1, -1 )
NBLOCK = MIN( NMAX, MAX( 1, NBLOCK ) )
*
* 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 = -23
ELSE IF( 2*MAX( 2, NMAX )**2.GT.LWORK ) THEN
INFO = -29
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DCHKST', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
$ RETURN
*
* More Important constants
*
UNFL = DLAMCH( 'Safe minimum' )
OVFL = ONE / UNFL
CALL DLABAD( UNFL, OVFL )
ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
ULPINV = ONE / ULP
LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) )
RTUNFL = SQRT( UNFL )
RTOVFL = SQRT( OVFL )
*
* Loop over sizes, types
*
DO 20 I = 1, 4
ISEED2( I ) = ISEED( I )
20 CONTINUE
NERRS = 0
NMATS = 0
*
DO 310 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 = 1 + 4*N + 2*N*LGN + 3*N**2
LIWEDC = 6 + 6*N + 5*N*LGN
ELSE
LWEDC = 8
LIWEDC = 12
END IF
NAP = ( N*( N+1 ) ) / 2
ANINV = ONE / DBLE( MAX( 1, N ) )
*
IF( NSIZES.NE.1 ) THEN
MTYPES = MIN( MAXTYP, NTYPES )
ELSE
MTYPES = MIN( MAXTYP+1, NTYPES )
END IF
*
DO 300 JTYPE = 1, MTYPES
IF( .NOT.DOTYPE( JTYPE ) )
$ GO TO 300
NMATS = NMATS + 1
NTEST = 0
*
DO 30 J = 1, 4
IOLDSD( J ) = ISEED( J )
30 CONTINUE
*
* 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 symmetric, w/ eigenvalues
* =6 random (none)
* =7 random diagonal
* =8 random symmetric
* =9 positive definite
* =10 diagonally dominant tridiagonal
*
IF( MTYPES.GT.MAXTYP )
$ GO TO 100
*
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 DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
IINFO = 0
IF( JTYPE.LE.15 ) THEN
COND = ULPINV
ELSE
COND = ULPINV*ANINV / TEN
END IF
*
* Special Matrices -- Identity & Jordan block
*
* Zero
*
IF( ITYPE.EQ.1 ) THEN
IINFO = 0
*
ELSE IF( ITYPE.EQ.2 ) THEN
*
* Identity
*
DO 80 JC = 1, N
A( JC, JC ) = ANORM
80 CONTINUE
*
ELSE IF( ITYPE.EQ.4 ) THEN
*
* Diagonal Matrix, [Eigen]values Specified
*
CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
$ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
$ IINFO )
*
*
ELSE IF( ITYPE.EQ.5 ) THEN
*
* Symmetric, eigenvalues specified
*
CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
$ ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
$ IINFO )
*
ELSE IF( ITYPE.EQ.7 ) THEN
*
* Diagonal, random eigenvalues
*
CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
$ '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
*
* Symmetric, random eigenvalues
*
CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
$ '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
*
* Positive definite, eigenvalues specified.
*
CALL DLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND,
$ ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
$ IINFO )
*
ELSE IF( ITYPE.EQ.10 ) THEN
*
* Positive definite tridiagonal, eigenvalues specified.
*
CALL DLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND,
$ ANORM, 1, 1, 'N', A, LDA, WORK( N+1 ),
$ IINFO )
DO 90 I = 2, N
TEMP1 = ABS( A( I-1, I ) ) /
$ SQRT( ABS( A( I-1, I-1 )*A( I, I ) ) )
IF( TEMP1.GT.HALF ) THEN
A( I-1, I ) = HALF*SQRT( ABS( A( I-1, I-1 )*A( I,
$ I ) ) )
A( I, I-1 ) = A( I-1, I )
END IF
90 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
*
100 CONTINUE
*
* Call DSYTRD and DORGTR to compute S and U from
* upper triangle.
*
CALL DLACPY( 'U', N, N, A, LDA, V, LDU )
*
NTEST = 1
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?