ddrvst.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 1,693 行 · 第 1/5 页
F
1,693 行
* 19= | A - U S U' | / ( |A| n ulp ) DSTEVR('V','I', ... )
* 20= | I - U U' | / ( n ulp ) DSTEVR('V','I', ... )
* 21= |D(with Z) - D(w/o Z)| / (|D| ulp) DSTEVR('N','I', ... )
* 22= | A - U S U' | / ( |A| n ulp ) DSTEVR('V','V', ... )
* 23= | I - U U' | / ( n ulp ) DSTEVR('V','V', ... )
* 24= |D(with Z) - D(w/o Z)| / (|D| ulp) DSTEVR('N','V', ... )
*
* 25= | A - U S U' | / ( |A| n ulp ) DSYEV('L','V', ... )
* 26= | I - U U' | / ( n ulp ) DSYEV('L','V', ... )
* 27= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEV('L','N', ... )
* 28= | A - U S U' | / ( |A| n ulp ) DSYEVX('L','V','A', ... )
* 29= | I - U U' | / ( n ulp ) DSYEVX('L','V','A', ... )
* 30= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVX('L','N','A', ... )
* 31= | A - U S U' | / ( |A| n ulp ) DSYEVX('L','V','I', ... )
* 32= | I - U U' | / ( n ulp ) DSYEVX('L','V','I', ... )
* 33= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVX('L','N','I', ... )
* 34= | A - U S U' | / ( |A| n ulp ) DSYEVX('L','V','V', ... )
* 35= | I - U U' | / ( n ulp ) DSYEVX('L','V','V', ... )
* 36= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVX('L','N','V', ... )
* 37= | A - U S U' | / ( |A| n ulp ) DSPEV('L','V', ... )
* 38= | I - U U' | / ( n ulp ) DSPEV('L','V', ... )
* 39= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEV('L','N', ... )
* 40= | A - U S U' | / ( |A| n ulp ) DSPEVX('L','V','A', ... )
* 41= | I - U U' | / ( n ulp ) DSPEVX('L','V','A', ... )
* 42= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVX('L','N','A', ... )
* 43= | A - U S U' | / ( |A| n ulp ) DSPEVX('L','V','I', ... )
* 44= | I - U U' | / ( n ulp ) DSPEVX('L','V','I', ... )
* 45= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVX('L','N','I', ... )
* 46= | A - U S U' | / ( |A| n ulp ) DSPEVX('L','V','V', ... )
* 47= | I - U U' | / ( n ulp ) DSPEVX('L','V','V', ... )
* 48= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVX('L','N','V', ... )
* 49= | A - U S U' | / ( |A| n ulp ) DSBEV('L','V', ... )
* 50= | I - U U' | / ( n ulp ) DSBEV('L','V', ... )
* 51= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEV('L','N', ... )
* 52= | A - U S U' | / ( |A| n ulp ) DSBEVX('L','V','A', ... )
* 53= | I - U U' | / ( n ulp ) DSBEVX('L','V','A', ... )
* 54= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVX('L','N','A', ... )
* 55= | A - U S U' | / ( |A| n ulp ) DSBEVX('L','V','I', ... )
* 56= | I - U U' | / ( n ulp ) DSBEVX('L','V','I', ... )
* 57= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVX('L','N','I', ... )
* 58= | A - U S U' | / ( |A| n ulp ) DSBEVX('L','V','V', ... )
* 59= | I - U U' | / ( n ulp ) DSBEVX('L','V','V', ... )
* 60= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVX('L','N','V', ... )
* 61= | A - U S U' | / ( |A| n ulp ) DSYEVD('L','V', ... )
* 62= | I - U U' | / ( n ulp ) DSYEVD('L','V', ... )
* 63= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVD('L','N', ... )
* 64= | A - U S U' | / ( |A| n ulp ) DSPEVD('L','V', ... )
* 65= | I - U U' | / ( n ulp ) DSPEVD('L','V', ... )
* 66= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVD('L','N', ... )
* 67= | A - U S U' | / ( |A| n ulp ) DSBEVD('L','V', ... )
* 68= | I - U U' | / ( n ulp ) DSBEVD('L','V', ... )
* 69= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVD('L','N', ... )
* 70= | A - U S U' | / ( |A| n ulp ) DSYEVR('L','V','A', ... )
* 71= | I - U U' | / ( n ulp ) DSYEVR('L','V','A', ... )
* 72= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVR('L','N','A', ... )
* 73= | A - U S U' | / ( |A| n ulp ) DSYEVR('L','V','I', ... )
* 74= | I - U U' | / ( n ulp ) DSYEVR('L','V','I', ... )
* 75= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVR('L','N','I', ... )
* 76= | A - U S U' | / ( |A| n ulp ) DSYEVR('L','V','V', ... )
* 77= | I - U U' | / ( n ulp ) DSYEVR('L','V','V', ... )
* 78= |D(with Z) - D(w/o Z)| / (|D| ulp) DSYEVR('L','N','V', ... )
*
* Tests 25 through 78 are repeated (as tests 79 through 132)
* with UPLO='U'
*
* To be added in 1999
*
* 79= | A - U S U' | / ( |A| n ulp ) DSPEVR('L','V','A', ... )
* 80= | I - U U' | / ( n ulp ) DSPEVR('L','V','A', ... )
* 81= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVR('L','N','A', ... )
* 82= | A - U S U' | / ( |A| n ulp ) DSPEVR('L','V','I', ... )
* 83= | I - U U' | / ( n ulp ) DSPEVR('L','V','I', ... )
* 84= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVR('L','N','I', ... )
* 85= | A - U S U' | / ( |A| n ulp ) DSPEVR('L','V','V', ... )
* 86= | I - U U' | / ( n ulp ) DSPEVR('L','V','V', ... )
* 87= |D(with Z) - D(w/o Z)| / (|D| ulp) DSPEVR('L','N','V', ... )
* 88= | A - U S U' | / ( |A| n ulp ) DSBEVR('L','V','A', ... )
* 89= | I - U U' | / ( n ulp ) DSBEVR('L','V','A', ... )
* 90= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVR('L','N','A', ... )
* 91= | A - U S U' | / ( |A| n ulp ) DSBEVR('L','V','I', ... )
* 92= | I - U U' | / ( n ulp ) DSBEVR('L','V','I', ... )
* 93= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVR('L','N','I', ... )
* 94= | A - U S U' | / ( |A| n ulp ) DSBEVR('L','V','V', ... )
* 95= | I - U U' | / ( n ulp ) DSBEVR('L','V','V', ... )
* 96= |D(with Z) - D(w/o Z)| / (|D| ulp) DSBEVR('L','N','V', ... )
*
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE, TWO, TEN
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
$ TEN = 10.0D0 )
DOUBLE PRECISION HALF
PARAMETER ( HALF = 0.5D0 )
INTEGER MAXTYP
PARAMETER ( MAXTYP = 18 )
* ..
* .. Local Scalars ..
LOGICAL BADNN
CHARACTER UPLO
INTEGER I, IDIAG, IHBW, IINFO, IL, IMODE, INDX, IROW,
$ ITEMP, ITYPE, IU, IUPLO, J, J1, J2, JCOL,
$ JSIZE, JTYPE, KD, LGN, LIWEDC, LWEDC, M, M2,
$ M3, MTYPES, N, NERRS, NMATS, NMAX, NTEST,
$ NTESTT
DOUBLE PRECISION ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL,
$ RTUNFL, TEMP1, TEMP2, TEMP3, ULP, ULPINV, UNFL,
$ VL, VU
* ..
* .. Local Arrays ..
INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
$ ISEED3( 4 ), KMAGN( MAXTYP ), KMODE( MAXTYP ),
$ KTYPE( MAXTYP )
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH, DLARND, DSXT1
EXTERNAL DLAMCH, DLARND, DSXT1
* ..
* .. External Subroutines ..
EXTERNAL ALASVM, DLABAD, DLACPY, DLAFTS, DLASET, DLATMR,
$ DLATMS, DSBEV, DSBEVD, DSBEVX, DSPEV, DSPEVD,
$ DSPEVX, DSTEV, DSTEVD, DSTEVR, DSTEVX, DSTT21,
$ DSTT22, DSYEV, DSYEVD, DSYEVR, DSYEVX, DSYT21,
$ DSYT22, XERBLA
* ..
* .. Scalars in Common ..
CHARACTER*6 SRNAMT
* ..
* .. Common blocks ..
COMMON / SRNAMC / SRNAMT
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, INT, LOG, MAX, MIN, SQRT
* ..
* .. Data statements ..
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 ..
*
* Keep ftrnchek happy
*
VL = ZERO
VU = ZERO
*
* 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 = -21
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DDRVST', -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 1740 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 + 4*N**2
c LIWEDC = 6 + 6*N + 5*N*LGN
LIWEDC = 3 + 5*N
ELSE
LWEDC = 9
c LIWEDC = 12
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 1730 JTYPE = 1, MTYPES
*
IF( .NOT.DOTYPE( JTYPE ) )
$ GO TO 1730
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 symmetric, w/ eigenvalues
* =6 random (none)
* =7 random diagonal
* =8 random symmetric
* =9 band symmetric, 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 DLASET( 'Full', LDA, N, ZERO, ZERO, 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 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
*
IDUMMA( 1 ) = 1
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
*
IDUMMA( 1 ) = 1
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
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?