cdrvev.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 835 行 · 第 1/3 页
F
835 行
*
* OVFL, UNFL Overflow and underflow thresholds.
* ULP, ULPINV Finest relative precision and its inverse.
* 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) Selectw whether CONDS is to be 1 or
* 1/sqrt(ulp). (0 means irrelevant.)
*
* =====================================================================
*
* .. Parameters ..
COMPLEX CZERO
PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )
COMPLEX CONE
PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
REAL TWO
PARAMETER ( TWO = 2.0E+0 )
INTEGER MAXTYP
PARAMETER ( MAXTYP = 21 )
* ..
* .. Local Scalars ..
LOGICAL BADNN
CHARACTER*3 PATH
INTEGER IINFO, IMODE, ITYPE, IWK, J, JCOL, JJ, JSIZE,
$ JTYPE, MTYPES, N, NERRS, NFAIL, NMAX,
$ NNWORK, NTEST, NTESTF, NTESTT
REAL ANORM, COND, CONDS, OVFL, RTULP, RTULPI, TNRM,
$ ULP, ULPINV, UNFL, VMX, VRMX, VTST
* ..
* .. Local Arrays ..
INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
$ KMAGN( MAXTYP ), KMODE( MAXTYP ),
$ KTYPE( MAXTYP )
REAL RES( 2 )
COMPLEX DUM( 1 )
* ..
* .. External Functions ..
REAL SCNRM2, SLAMCH
EXTERNAL SCNRM2, SLAMCH
* ..
* .. External Subroutines ..
EXTERNAL CGEEV, CGET22, CLACPY, CLATME, CLATMR, CLATMS,
$ CLASET, SLABAD, SLASUM, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, AIMAG, CMPLX, MAX, MIN, REAL, 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 ..
*
PATH( 1: 1 ) = 'Complex precision'
PATH( 2: 3 ) = 'EV'
*
* Check for errors
*
NTESTT = 0
NTESTF = 0
INFO = 0
*
* Important constants
*
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( NOUNIT.LE.0 ) THEN
INFO = -7
ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
INFO = -9
ELSE IF( LDVL.LT.1 .OR. LDVL.LT.NMAX ) THEN
INFO = -14
ELSE IF( LDVR.LT.1 .OR. LDVR.LT.NMAX ) THEN
INFO = -16
ELSE IF( LDLRE.LT.1 .OR. LDLRE.LT.NMAX ) THEN
INFO = -28
ELSE IF( 5*NMAX+2*NMAX**2.GT.NWORK ) THEN
INFO = -21
END IF
*
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'CDRVEV', -INFO )
RETURN
END IF
*
* Quick return if nothing to do
*
IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
$ RETURN
*
* More Important constants
*
UNFL = SLAMCH( 'Safe minimum' )
OVFL = ONE / UNFL
CALL SLABAD( UNFL, OVFL )
ULP = SLAMCH( 'Precision' )
ULPINV = ONE / ULP
RTULP = SQRT( ULP )
RTULPI = ONE / RTULP
*
* Loop over sizes, types
*
NERRS = 0
*
DO 270 JSIZE = 1, NSIZES
N = NN( JSIZE )
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
*
* Save ISEED in case of an error.
*
DO 20 J = 1, 4
IOLDSD( J ) = ISEED( J )
20 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 90
*
ITYPE = KTYPE( JTYPE )
IMODE = KMODE( JTYPE )
*
* Compute norm
*
GO TO ( 30, 40, 50 )KMAGN( JTYPE )
*
30 CONTINUE
ANORM = ONE
GO TO 60
*
40 CONTINUE
ANORM = OVFL*ULP
GO TO 60
*
50 CONTINUE
ANORM = UNFL*ULPINV
GO TO 60
*
60 CONTINUE
*
CALL CLASET( '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 70 JCOL = 1, N
A( JCOL, JCOL ) = CMPLX( ANORM )
70 CONTINUE
*
ELSE IF( ITYPE.EQ.3 ) THEN
*
* Jordan Block
*
DO 80 JCOL = 1, N
A( JCOL, JCOL ) = CMPLX( ANORM )
IF( JCOL.GT.1 )
$ A( JCOL, JCOL-1 ) = CONE
80 CONTINUE
*
ELSE IF( ITYPE.EQ.4 ) THEN
*
* Diagonal Matrix, [Eigen]values Specified
*
CALL CLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
$ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
$ IINFO )
*
ELSE IF( ITYPE.EQ.5 ) THEN
*
* Hermitian, eigenvalues specified
*
CALL CLATMS( N, N, 'S', ISEED, 'H', RWORK, 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
*
CALL CLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE,
$ ' ', 'T', 'T', 'T', RWORK, 4, CONDS, N, N,
$ ANORM, A, LDA, WORK( 2*N+1 ), IINFO )
*
ELSE IF( ITYPE.EQ.7 ) THEN
*
* Diagonal, random eigenvalues
*
CALL CLATMR( N, N, 'D', ISEED, 'N', 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
*
* Symmetric, random eigenvalues
*
CALL CLATMR( N, N, 'D', 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
*
* General, random eigenvalues
*
CALL CLATMR( N, N, 'D', ISEED, 'N', 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 )
IF( N.GE.4 ) THEN
CALL CLASET( 'Full', 2, N, CZERO, CZERO, A, LDA )
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?