📄 zchkst.f
字号:
SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
$ NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5,
$ WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK,
$ LWORK, RWORK, LRWORK, IWORK, LIWORK, RESULT,
$ INFO )
IMPLICIT NONE
*
* -- LAPACK test routine (version 3.1) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* November 2006
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LDU, LIWORK, LRWORK, LWORK, NOUNIT,
$ NSIZES, NTYPES
DOUBLE PRECISION THRESH
* ..
* .. Array Arguments ..
LOGICAL DOTYPE( * )
INTEGER ISEED( 4 ), IWORK( * ), NN( * )
DOUBLE PRECISION D1( * ), D2( * ), D3( * ), D4( * ), D5( * ),
$ RESULT( * ), RWORK( * ), SD( * ), SE( * ),
$ WA1( * ), WA2( * ), WA3( * ), WR( * )
COMPLEX*16 A( LDA, * ), AP( * ), TAU( * ), U( LDU, * ),
$ V( LDU, * ), VP( * ), WORK( * ), Z( LDU, * )
* ..
*
* Purpose
* =======
*
* ZCHKST checks the Hermitian eigenvalue problem routines.
*
* ZHETRD factors A as U S U* , where * means conjugate transpose,
* S is real symmetric tridiagonal, and U is unitary.
* ZHETRD can use either just the lower or just the upper triangle
* of A; ZCHKST checks both cases.
* U is represented as a product of Householder
* transformations, whose vectors are stored in the first
* n-1 columns of V, and whose scale factors are in TAU.
*
* ZHPTRD does the same as ZHETRD, except that A and V are stored
* in "packed" format.
*
* ZUNGTR constructs the matrix U from the contents of V and TAU.
*
* ZUPGTR constructs the matrix U from the contents of VP and TAU.
*
* ZSTEQR factors S as Z D1 Z* , where Z is the unitary
* matrix of eigenvectors and D1 is a diagonal matrix with
* the eigenvalues on the diagonal. D2 is the matrix of
* eigenvalues computed when Z is not computed.
*
* DSTERF computes D3, the matrix of eigenvalues, by the
* PWK method, which does not yield eigenvectors.
*
* ZPTEQR factors S as Z4 D4 Z4* , for a
* Hermitian positive definite tridiagonal matrix.
* D5 is the matrix of eigenvalues computed when Z is not
* computed.
*
* DSTEBZ computes selected eigenvalues. WA1, WA2, and
* WA3 will denote eigenvalues computed to high
* absolute accuracy, with different range options.
* WR will denote eigenvalues computed to high relative
* accuracy.
*
* ZSTEIN computes Y, the eigenvectors of S, given the
* eigenvalues.
*
* ZSTEDC factors S as Z D1 Z* , where Z is the unitary
* matrix of eigenvectors and D1 is a diagonal matrix with
* the eigenvalues on the diagonal ('I' option). It may also
* update an input unitary matrix, usually the output
* from ZHETRD/ZUNGTR or ZHPTRD/ZUPGTR ('V' option). It may
* also just compute eigenvalues ('N' option).
*
* ZSTEMR factors S as Z D1 Z* , where Z is the unitary
* matrix of eigenvectors and D1 is a diagonal matrix with
* the eigenvalues on the diagonal ('I' option). ZSTEMR
* uses the Relatively Robust Representation whenever possible.
*
* When ZCHKST is called, a number of matrix "sizes" ("n's") and a
* number of matrix "types" are specified. For each size ("n")
* and each type of matrix, one matrix will be generated and used
* to test the Hermitian eigenroutines. For each matrix, a number
* of tests will be performed:
*
* (1) | A - V S V* | / ( |A| n ulp ) ZHETRD( UPLO='U', ... )
*
* (2) | I - UV* | / ( n ulp ) ZUNGTR( UPLO='U', ... )
*
* (3) | A - V S V* | / ( |A| n ulp ) ZHETRD( UPLO='L', ... )
*
* (4) | I - UV* | / ( n ulp ) ZUNGTR( UPLO='L', ... )
*
* (5-8) Same as 1-4, but for ZHPTRD and ZUPGTR.
*
* (9) | S - Z D Z* | / ( |S| n ulp ) ZSTEQR('V',...)
*
* (10) | I - ZZ* | / ( n ulp ) ZSTEQR('V',...)
*
* (11) | D1 - D2 | / ( |D1| ulp ) ZSTEQR('N',...)
*
* (12) | D1 - D3 | / ( |D1| ulp ) DSTERF
*
* (13) 0 if the true eigenvalues (computed by sturm count)
* of S are within THRESH of
* those in D1. 2*THRESH if they are not. (Tested using
* DSTECH)
*
* For S positive definite,
*
* (14) | S - Z4 D4 Z4* | / ( |S| n ulp ) ZPTEQR('V',...)
*
* (15) | I - Z4 Z4* | / ( n ulp ) ZPTEQR('V',...)
*
* (16) | D4 - D5 | / ( 100 |D4| ulp ) ZPTEQR('N',...)
*
* When S is also diagonally dominant by the factor gamma < 1,
*
* (17) max | D4(i) - WR(i) | / ( |D4(i)| omega ) ,
* i
* omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
* DSTEBZ( 'A', 'E', ...)
*
* (18) | WA1 - D3 | / ( |D3| ulp ) DSTEBZ( 'A', 'E', ...)
*
* (19) ( max { min | WA2(i)-WA3(j) | } +
* i j
* max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
* i j
* DSTEBZ( 'I', 'E', ...)
*
* (20) | S - Y WA1 Y* | / ( |S| n ulp ) DSTEBZ, ZSTEIN
*
* (21) | I - Y Y* | / ( n ulp ) DSTEBZ, ZSTEIN
*
* (22) | S - Z D Z* | / ( |S| n ulp ) ZSTEDC('I')
*
* (23) | I - ZZ* | / ( n ulp ) ZSTEDC('I')
*
* (24) | S - Z D Z* | / ( |S| n ulp ) ZSTEDC('V')
*
* (25) | I - ZZ* | / ( n ulp ) ZSTEDC('V')
*
* (26) | D1 - D2 | / ( |D1| ulp ) ZSTEDC('V') and
* ZSTEDC('N')
*
* Test 27 is disabled at the moment because ZSTEMR does not
* guarantee high relatvie accuracy.
*
* (27) max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
* i
* omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
* ZSTEMR('V', 'A')
*
* (28) max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
* i
* omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
* ZSTEMR('V', 'I')
*
* Tests 29 through 34 are disable at present because ZSTEMR
* does not handle partial specturm requests.
*
* (29) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'I')
*
* (30) | I - ZZ* | / ( n ulp ) ZSTEMR('V', 'I')
*
* (31) ( max { min | WA2(i)-WA3(j) | } +
* i j
* max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
* i j
* ZSTEMR('N', 'I') vs. CSTEMR('V', 'I')
*
* (32) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'V')
*
* (33) | I - ZZ* | / ( n ulp ) ZSTEMR('V', 'V')
*
* (34) ( max { min | WA2(i)-WA3(j) | } +
* i j
* max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
* i j
* ZSTEMR('N', 'V') vs. CSTEMR('V', 'V')
*
* (35) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'A')
*
* (36) | I - ZZ* | / ( n ulp ) ZSTEMR('V', 'A')
*
* (37) ( max { min | WA2(i)-WA3(j) | } +
* i j
* max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
* i j
* ZSTEMR('N', 'A') vs. CSTEMR('V', 'A')
*
* The "sizes" are specified by an array NN(1:NSIZES); the value of
* each element NN(j) specifies one size.
* The "types" are specified by a logical array DOTYPE( 1:NTYPES );
* if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
* Currently, the list of possible types is:
*
* (1) The zero matrix.
* (2) The identity matrix.
*
* (3) A diagonal matrix with evenly spaced entries
* 1, ..., ULP and random signs.
* (ULP = (first number larger than 1) - 1 )
* (4) A diagonal matrix with geometrically spaced entries
* 1, ..., ULP and random signs.
* (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
* and random signs.
*
* (6) Same as (4), but multiplied by SQRT( overflow threshold )
* (7) Same as (4), but multiplied by SQRT( underflow threshold )
*
* (8) A matrix of the form U* D U, where U is unitary and
* D has evenly spaced entries 1, ..., ULP with random signs
* on the diagonal.
*
* (9) A matrix of the form U* D U, where U is unitary and
* D has geometrically spaced entries 1, ..., ULP with random
* signs on the diagonal.
*
* (10) A matrix of the form U* D U, where U is unitary and
* D has "clustered" entries 1, ULP,..., ULP with random
* signs on the diagonal.
*
* (11) Same as (8), but multiplied by SQRT( overflow threshold )
* (12) Same as (8), but multiplied by SQRT( underflow threshold )
*
* (13) Hermitian matrix with random entries chosen from (-1,1).
* (14) Same as (13), but multiplied by SQRT( overflow threshold )
* (15) Same as (13), but multiplied by SQRT( underflow threshold )
* (16) Same as (8), but diagonal elements are all positive.
* (17) Same as (9), but diagonal elements are all positive.
* (18) Same as (10), but diagonal elements are all positive.
* (19) Same as (16), but multiplied by SQRT( overflow threshold )
* (20) Same as (16), but multiplied by SQRT( underflow threshold )
* (21) A diagonally dominant tridiagonal matrix with geometrically
* spaced diagonal entries 1, ..., ULP.
*
* Arguments
* =========
*
* NSIZES (input) INTEGER
* The number of sizes of matrices to use. If it is zero,
* ZCHKST does nothing. It must be at least zero.
*
* NN (input) INTEGER array, dimension (NSIZES)
* An array containing the sizes to be used for the matrices.
* Zero values will be skipped. The values must be at least
* zero.
*
* NTYPES (input) INTEGER
* The number of elements in DOTYPE. If it is zero, ZCHKST
* does nothing. It must be at least zero. If it is MAXTYP+1
* and NSIZES is 1, then an additional type, MAXTYP+1 is
* defined, which is to use whatever matrix is in A. This
* is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
* DOTYPE(MAXTYP+1) is .TRUE. .
*
* DOTYPE (input) LOGICAL array, dimension (NTYPES)
* If DOTYPE(j) is .TRUE., then for each size in NN a
* matrix of that size and of type j will be generated.
* If NTYPES is smaller than the maximum number of types
* defined (PARAMETER MAXTYP), then types NTYPES+1 through
* MAXTYP will not be generated. If NTYPES is larger
* than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
* will be ignored.
*
* ISEED (input/output) INTEGER array, dimension (4)
* On entry ISEED specifies the seed of the random number
* generator. The array elements should be between 0 and 4095;
* if not they will be reduced mod 4096. Also, ISEED(4) must
* be odd. The random number generator uses a linear
* congruential sequence limited to small integers, and so
* should produce machine independent random numbers. The
* values of ISEED are changed on exit, and can be used in the
* next call to ZCHKST to continue the same random number
* sequence.
*
* THRESH (input) DOUBLE PRECISION
* A test will count as "failed" if the "error", computed as
* described above, exceeds THRESH. Note that the error
* is scaled to be O(1), so THRESH should be a reasonably
* small multiple of 1, e.g., 10 or 100. In particular,
* it should not depend on the precision (single vs. double)
* or the size of the matrix. It must be at least zero.
*
* NOUNIT (input) INTEGER
* The FORTRAN unit number for printing out error messages
* (e.g., if a routine returns IINFO not equal to 0.)
*
* A (input/workspace/output) COMPLEX*16 array of
* dimension ( LDA , max(NN) )
* Used to hold the matrix whose eigenvalues are to be
* computed. On exit, A contains the last matrix actually
* used.
*
* LDA (input) INTEGER
* The leading dimension of A. It must be at
* least 1 and at least max( NN ).
*
* AP (workspace) COMPLEX*16 array of
* dimension( max(NN)*max(NN+1)/2 )
* The matrix A stored in packed format.
*
* SD (workspace/output) DOUBLE PRECISION array of
* dimension( max(NN) )
* The diagonal of the tridiagonal matrix computed by ZHETRD.
* On exit, SD and SE contain the tridiagonal form of the
* matrix in A.
*
* SE (workspace/output) DOUBLE PRECISION array of
* dimension( max(NN) )
* The off-diagonal of the tridiagonal matrix computed by
* ZHETRD. On exit, SD and SE contain the tridiagonal form of
* the matrix in A.
*
* D1 (workspace/output) DOUBLE PRECISION array of
* dimension( max(NN) )
* The eigenvalues of A, as computed by ZSTEQR simlutaneously
* with Z. On exit, the eigenvalues in D1 correspond with the
* matrix in A.
*
* D2 (workspace/output) DOUBLE PRECISION array of
* dimension( max(NN) )
* The eigenvalues of A, as computed by ZSTEQR if Z is not
* computed. On exit, the eigenvalues in D2 correspond with
* the matrix in A.
*
* D3 (workspace/output) DOUBLE PRECISION array of
* dimension( max(NN) )
* The eigenvalues of A, as computed by DSTERF. On exit, the
* eigenvalues in D3 correspond with the matrix in A.
*
* U (workspace/output) COMPLEX*16 array of
* dimension( LDU, max(NN) ).
* The unitary matrix computed by ZHETRD + ZUNGTR.
*
* LDU (input) INTEGER
* The leading dimension of U, Z, and V. It must be at least 1
* and at least max( NN ).
*
* V (workspace/output) COMPLEX*16 array of
* dimension( LDU, max(NN) ).
* The Housholder vectors computed by ZHETRD in reducing A to
* tridiagonal form. The vectors computed with UPLO='U' are
* in the upper triangle, and the vectors computed with UPLO='L'
* are in the lower triangle. (As described in ZHETRD, the
* sub- and superdiagonal are not set to 1, although the
* true Householder vector has a 1 in that position. The
* routines that use V, such as ZUNGTR, set those entries to
* 1 before using them, and then restore them later.)
*
* VP (workspace) COMPLEX*16 array of
* dimension( max(NN)*max(NN+1)/2 )
* The matrix V stored in packed format.
*
* TAU (workspace/output) COMPLEX*16 array of
* dimension( max(NN) )
* The Householder factors computed by ZHETRD in reducing A
* to tridiagonal form.
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -