cdrves.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 803 行 · 第 1/2 页
F
803 行
ULPINV = ONE / ULP
RTULP = SQRT( ULP )
RTULPI = ONE / RTULP
*
* Loop over sizes, types
*
NERRS = 0
*
DO 240 JSIZE = 1, NSIZES
N = NN( JSIZE )
IF( NSIZES.NE.1 ) THEN
MTYPES = MIN( MAXTYP, NTYPES )
ELSE
MTYPES = MIN( MAXTYP+1, NTYPES )
END IF
*
DO 230 JTYPE = 1, MTYPES
IF( .NOT.DOTYPE( JTYPE ) )
$ GO TO 230
*
* 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
*
IF( ITYPE.EQ.1 ) THEN
*
* Zero
*
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
*
* Symmetric, 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 )
CALL CLASET( 'Full', N-3, 1, CZERO, CZERO, A( 3, 1 ),
$ LDA )
CALL CLASET( 'Full', N-3, 2, CZERO, CZERO,
$ A( 3, N-1 ), LDA )
CALL CLASET( 'Full', 1, N, CZERO, CZERO, A( N, 1 ),
$ LDA )
END IF
*
ELSE IF( ITYPE.EQ.10 ) THEN
*
* Triangular, 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, 0,
$ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
*
ELSE
*
IINFO = 1
END IF
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9992 )'Generator', IINFO, N, JTYPE,
$ IOLDSD
INFO = ABS( IINFO )
RETURN
END IF
*
90 CONTINUE
*
* Test for minimal and generous workspace
*
DO 220 IWK = 1, 2
IF( IWK.EQ.1 ) THEN
NNWORK = 3*N
ELSE
NNWORK = 5*N + 2*N**2
END IF
NNWORK = MAX( NNWORK, 1 )
*
* Initialize RESULT
*
DO 100 J = 1, 13
RESULT( J ) = -ONE
100 CONTINUE
*
* Test with and without sorting of eigenvalues
*
DO 180 ISORT = 0, 1
IF( ISORT.EQ.0 ) THEN
SORT = 'N'
RSUB = 0
ELSE
SORT = 'S'
RSUB = 6
END IF
*
* Compute Schur form and Schur vectors, and test them
*
CALL CLACPY( 'F', N, N, A, LDA, H, LDA )
CALL CGEES( 'V', SORT, CSLECT, N, H, LDA, SDIM, W, VS,
$ LDVS, WORK, NNWORK, RWORK, BWORK, IINFO )
IF( IINFO.NE.0 ) THEN
RESULT( 1+RSUB ) = ULPINV
WRITE( NOUNIT, FMT = 9992 )'CGEES1', IINFO, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
GO TO 190
END IF
*
* Do Test (1) or Test (7)
*
RESULT( 1+RSUB ) = ZERO
DO 120 J = 1, N - 1
DO 110 I = J + 1, N
IF( H( I, J ).NE.ZERO )
$ RESULT( 1+RSUB ) = ULPINV
110 CONTINUE
120 CONTINUE
*
* Do Tests (2) and (3) or Tests (8) and (9)
*
LWORK = MAX( 1, 2*N*N )
CALL CHST01( N, 1, N, A, LDA, H, LDA, VS, LDVS, WORK,
$ LWORK, RWORK, RES )
RESULT( 2+RSUB ) = RES( 1 )
RESULT( 3+RSUB ) = RES( 2 )
*
* Do Test (4) or Test (10)
*
RESULT( 4+RSUB ) = ZERO
DO 130 I = 1, N
IF( H( I, I ).NE.W( I ) )
$ RESULT( 4+RSUB ) = ULPINV
130 CONTINUE
*
* Do Test (5) or Test (11)
*
CALL CLACPY( 'F', N, N, A, LDA, HT, LDA )
CALL CGEES( 'N', SORT, CSLECT, N, HT, LDA, SDIM, WT,
$ VS, LDVS, WORK, NNWORK, RWORK, BWORK,
$ IINFO )
IF( IINFO.NE.0 ) THEN
RESULT( 5+RSUB ) = ULPINV
WRITE( NOUNIT, FMT = 9992 )'CGEES2', IINFO, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
GO TO 190
END IF
*
RESULT( 5+RSUB ) = ZERO
DO 150 J = 1, N
DO 140 I = 1, N
IF( H( I, J ).NE.HT( I, J ) )
$ RESULT( 5+RSUB ) = ULPINV
140 CONTINUE
150 CONTINUE
*
* Do Test (6) or Test (12)
*
RESULT( 6+RSUB ) = ZERO
DO 160 I = 1, N
IF( W( I ).NE.WT( I ) )
$ RESULT( 6+RSUB ) = ULPINV
160 CONTINUE
*
* Do Test (13)
*
IF( ISORT.EQ.1 ) THEN
RESULT( 13 ) = ZERO
KNTEIG = 0
DO 170 I = 1, N
IF( CSLECT( W( I ) ) )
$ KNTEIG = KNTEIG + 1
IF( I.LT.N ) THEN
IF( CSLECT( W( I+1 ) ) .AND.
$ ( .NOT.CSLECT( W( I ) ) ) )RESULT( 13 )
$ = ULPINV
END IF
170 CONTINUE
IF( SDIM.NE.KNTEIG )
$ RESULT( 13 ) = ULPINV
END IF
*
180 CONTINUE
*
* End of Loop -- Check for RESULT(j) > THRESH
*
190 CONTINUE
*
NTEST = 0
NFAIL = 0
DO 200 J = 1, 13
IF( RESULT( J ).GE.ZERO )
$ NTEST = NTEST + 1
IF( RESULT( J ).GE.THRESH )
$ NFAIL = NFAIL + 1
200 CONTINUE
*
IF( NFAIL.GT.0 )
$ NTESTF = NTESTF + 1
IF( NTESTF.EQ.1 ) THEN
WRITE( NOUNIT, FMT = 9999 )PATH
WRITE( NOUNIT, FMT = 9998 )
WRITE( NOUNIT, FMT = 9997 )
WRITE( NOUNIT, FMT = 9996 )
WRITE( NOUNIT, FMT = 9995 )THRESH
WRITE( NOUNIT, FMT = 9994 )
NTESTF = 2
END IF
*
DO 210 J = 1, 13
IF( RESULT( J ).GE.THRESH ) THEN
WRITE( NOUNIT, FMT = 9993 )N, IWK, IOLDSD, JTYPE,
$ J, RESULT( J )
END IF
210 CONTINUE
*
NERRS = NERRS + NFAIL
NTESTT = NTESTT + NTEST
*
220 CONTINUE
230 CONTINUE
240 CONTINUE
*
* Summary
*
CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT )
*
9999 FORMAT( / 1X, A3, ' -- Complex Schur Form Decomposition Driver',
$ / ' Matrix types (see CDRVES for details): ' )
*
9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
$ ' ', ' 5=Diagonal: geometr. spaced entries.',
$ / ' 2=Identity matrix. ', ' 6=Diagona',
$ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
$ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
$ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
$ 'mall, evenly spaced.' )
9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
$ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
$ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
$ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
$ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
$ 'lex ', A6, / ' 12=Well-cond., random complex ', A6, ' ',
$ ' 17=Ill-cond., large rand. complx ', A4, / ' 13=Ill-condi',
$ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
$ ' complx ', A4 )
9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
$ 'with small random entries.', / ' 20=Matrix with large ran',
$ 'dom entries. ', / )
9995 FORMAT( ' Tests performed with test threshold =', F8.2,
$ / ' ( A denotes A on input and T denotes A on output)',
$ / / ' 1 = 0 if T in Schur form (no sort), ',
$ ' 1/ulp otherwise', /
$ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
$ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ',
$ / ' 4 = 0 if W are eigenvalues of T (no sort),',
$ ' 1/ulp otherwise', /
$ ' 5 = 0 if T same no matter if VS computed (no sort),',
$ ' 1/ulp otherwise', /
$ ' 6 = 0 if W same no matter if VS computed (no sort)',
$ ', 1/ulp otherwise' )
9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
$ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
$ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
$ / ' 10 = 0 if W are eigenvalues of T (sort),',
$ ' 1/ulp otherwise', /
$ ' 11 = 0 if T same no matter if VS computed (sort),',
$ ' 1/ulp otherwise', /
$ ' 12 = 0 if W same no matter if VS computed (sort),',
$ ' 1/ulp otherwise', /
$ ' 13 = 0 if sorting succesful, 1/ulp otherwise', / )
9993 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ),
$ ' type ', I2, ', test(', I2, ')=', G10.3 )
9992 FORMAT( ' CDRVES: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
$ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
*
RETURN
*
* End of CDRVES
*
END
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?