dchkhs.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 1,085 行 · 第 1/3 页
F
1,085 行
* OVFL, UNFL Overflow and underflow thresholds.
* ULP, ULPINV Finest relative precision and its inverse.
* RTOVFL, RTUNFL,
* RTULP, RTULPI Square roots of the previous 4 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) )
* KCONDS(j) Selects whether CONDS is to be 1 or
* 1/sqrt(ulp). (0 means irrelevant.)
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
INTEGER MAXTYP
PARAMETER ( MAXTYP = 21 )
* ..
* .. Local Scalars ..
LOGICAL BADNN, MATCH
INTEGER I, IHI, IINFO, ILO, IMODE, IN, ITYPE, J, JCOL,
$ JJ, JSIZE, JTYPE, K, MTYPES, N, N1, NERRS,
$ NMATS, NMAX, NSELC, NSELR, NTEST, NTESTT
DOUBLE PRECISION ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP,
$ RTULPI, RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL
* ..
* .. Local Arrays ..
CHARACTER ADUMMA( 1 )
INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
$ KMAGN( MAXTYP ), KMODE( MAXTYP ),
$ KTYPE( MAXTYP )
DOUBLE PRECISION DUMMA( 6 )
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH
EXTERNAL DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL DCOPY, DGEHRD, DGEMM, DGET10, DGET22, DHSEIN,
$ DHSEQR, DHST01, DLABAD, DLACPY, DLAFTS, DLASET,
$ DLASUM, DLATME, DLATMR, DLATMS, DORGHR, DORMHR,
$ DTREVC, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX, MIN, SQRT
* ..
* .. Data statements ..
DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
$ 3, 1, 2, 3 /
DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
$ 1, 5, 5, 5, 4, 3, 1 /
DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
* ..
* .. Executable Statements ..
*
* Check for errors
*
NTESTT = 0
INFO = 0
*
BADNN = .FALSE.
NMAX = 0
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( THRESH.LT.ZERO ) THEN
INFO = -6
ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
INFO = -9
ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN
INFO = -14
ELSE IF( 4*NMAX*NMAX+2.GT.NWORK ) THEN
INFO = -28
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DCHKHS', -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 = DLAMCH( 'Overflow' )
CALL DLABAD( UNFL, OVFL )
ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
ULPINV = ONE / ULP
RTUNFL = SQRT( UNFL )
RTOVFL = SQRT( OVFL )
RTULP = SQRT( ULP )
RTULPI = ONE / RTULP
*
* Loop over sizes, types
*
NERRS = 0
NMATS = 0
*
DO 270 JSIZE = 1, NSIZES
N = NN( JSIZE )
IF( N.EQ.0 )
$ GO TO 270
N1 = MAX( 1, N )
ANINV = ONE / DBLE( N1 )
*
IF( NSIZES.NE.1 ) THEN
MTYPES = MIN( MAXTYP, NTYPES )
ELSE
MTYPES = MIN( MAXTYP+1, NTYPES )
END IF
*
DO 260 JTYPE = 1, MTYPES
IF( .NOT.DOTYPE( JTYPE ) )
$ GO TO 260
NMATS = NMATS + 1
NTEST = 0
*
* Save ISEED in case of an error.
*
DO 20 J = 1, 4
IOLDSD( J ) = ISEED( J )
20 CONTINUE
*
* Initialize RESULT
*
DO 30 J = 1, 14
RESULT( J ) = ZERO
30 CONTINUE
*
* Compute "A"
*
* Control parameters:
*
* KMAGN KCONDS KMODE KTYPE
* =1 O(1) 1 clustered 1 zero
* =2 large large clustered 2 identity
* =3 small exponential Jordan
* =4 arithmetic diagonal, (w/ eigenvalues)
* =5 random log symmetric, w/ eigenvalues
* =6 random general, w/ eigenvalues
* =7 random diagonal
* =8 random symmetric
* =9 random general
* =10 random triangular
*
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
COND = ULPINV
*
* Special Matrices
*
IF( ITYPE.EQ.1 ) THEN
*
* Zero
*
IINFO = 0
*
ELSE IF( ITYPE.EQ.2 ) THEN
*
* Identity
*
DO 80 JCOL = 1, N
A( JCOL, JCOL ) = ANORM
80 CONTINUE
*
ELSE IF( ITYPE.EQ.3 ) THEN
*
* Jordan Block
*
DO 90 JCOL = 1, N
A( JCOL, JCOL ) = ANORM
IF( JCOL.GT.1 )
$ A( JCOL, JCOL-1 ) = ONE
90 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.6 ) THEN
*
* General, eigenvalues specified
*
IF( KCONDS( JTYPE ).EQ.1 ) THEN
CONDS = ONE
ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
CONDS = RTULPI
ELSE
CONDS = ZERO
END IF
*
ADUMMA( 1 ) = ' '
CALL DLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE,
$ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4,
$ CONDS, N, N, ANORM, A, LDA, WORK( 2*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
*
* General, random eigenvalues
*
CALL DLATMR( N, N, 'S', ISEED, 'N', 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.10 ) THEN
*
* Triangular, random eigenvalues
*
CALL DLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE,
$ 'T', 'N', WORK( N+1 ), 1, ONE,
$ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
$ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
*
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 DGEHRD to compute H and U, do tests.
*
CALL DLACPY( ' ', N, N, A, LDA, H, LDA )
*
NTEST = 1
*
ILO = 1
IHI = N
*
CALL DGEHRD( N, ILO, IHI, H, LDA, WORK, WORK( N+1 ),
$ NWORK-N, IINFO )
*
IF( IINFO.NE.0 ) THEN
RESULT( 1 ) = ULPINV
WRITE( NOUNIT, FMT = 9999 )'DGEHRD', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
GO TO 250
END IF
*
DO 120 J = 1, N - 1
UU( J+1, J ) = ZERO
DO 110 I = J + 2, N
U( I, J ) = H( I, J )
UU( I, J ) = H( I, J )
H( I, J ) = ZERO
110 CONTINUE
120 CONTINUE
CALL DCOPY( N-1, WORK, 1, TAU, 1 )
CALL DORGHR( N, ILO, IHI, U, LDU, WORK, WORK( N+1 ),
$ NWORK-N, IINFO )
NTEST = 2
*
CALL DHST01( N, ILO, IHI, A, LDA, H, LDA, U, LDU, WORK,
$ NWORK, RESULT( 1 ) )
*
* Call DHSEQR to compute T1, T2 and Z, do tests.
*
* Eigenvalues only (WR3,WI3)
*
CALL DLACPY( ' ', N, N, H, LDA, T2, LDA )
NTEST = 3
RESULT( 3 ) = ULPINV
*
CALL DHSEQR( 'E', 'N', N, ILO, IHI, T2, LDA, WR3, WI3, UZ,
$ LDU, WORK, NWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'DHSEQR(E)', IINFO, N, JTYPE,
$ IOLDSD
IF( IINFO.LE.N+2 ) THEN
INFO = ABS( IINFO )
GO TO 250
END IF
END IF
*
* Eigenvalues (WR1,WI1) and Full Schur Form (T2)
*
CALL DLACPY( ' ', N, N, H, LDA, T2, LDA )
*
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?