schkst.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 1,804 行 · 第 1/5 页
F
1,804 行
SUBROUTINE SCHKST( 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, 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, LWORK, NOUNIT, NSIZES,
$ NTYPES
REAL THRESH
* ..
* .. Array Arguments ..
LOGICAL DOTYPE( * )
INTEGER ISEED( 4 ), IWORK( * ), NN( * )
REAL A( LDA, * ), AP( * ), D1( * ), D2( * ),
$ D3( * ), D4( * ), D5( * ), RESULT( * ),
$ SD( * ), SE( * ), TAU( * ), U( LDU, * ),
$ V( LDU, * ), VP( * ), WA1( * ), WA2( * ),
$ WA3( * ), WORK( * ), WR( * ), Z( LDU, * )
* ..
*
* Purpose
* =======
*
* SCHKST checks the symmetric eigenvalue problem routines.
*
* SSYTRD factors A as U S U' , where ' means transpose,
* S is symmetric tridiagonal, and U is orthogonal.
* SSYTRD can use either just the lower or just the upper triangle
* of A; SCHKST 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.
*
* SSPTRD does the same as SSYTRD, except that A and V are stored
* in "packed" format.
*
* SORGTR constructs the matrix U from the contents of V and TAU.
*
* SOPGTR constructs the matrix U from the contents of VP and TAU.
*
* SSTEQR factors S as Z D1 Z' , where Z is the orthogonal
* 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.
*
* SSTERF computes D3, the matrix of eigenvalues, by the
* PWK method, which does not yield eigenvectors.
*
* SPTEQR factors S as Z4 D4 Z4' , for a
* symmetric positive definite tridiagonal matrix.
* D5 is the matrix of eigenvalues computed when Z is not
* computed.
*
* SSTEBZ 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.
*
* SSTEIN computes Y, the eigenvectors of S, given the
* eigenvalues.
*
* SSTEDC factors S as Z D1 Z' , where Z is the orthogonal
* matrix of eigenvectors and D1 is a diagonal matrix with
* the eigenvalues on the diagonal ('I' option). It may also
* update an input orthogonal matrix, usually the output
* from SSYTRD/SORGTR or SSPTRD/SOPGTR ('V' option). It may
* also just compute eigenvalues ('N' option).
*
* SSTEMR factors S as Z D1 Z' , where Z is the orthogonal
* matrix of eigenvectors and D1 is a diagonal matrix with
* the eigenvalues on the diagonal ('I' option). SSTEMR
* uses the Relatively Robust Representation whenever possible.
*
* When SCHKST 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 symmetric eigenroutines. For each matrix, a number
* of tests will be performed:
*
* (1) | A - V S V' | / ( |A| n ulp ) SSYTRD( UPLO='U', ... )
*
* (2) | I - UV' | / ( n ulp ) SORGTR( UPLO='U', ... )
*
* (3) | A - V S V' | / ( |A| n ulp ) SSYTRD( UPLO='L', ... )
*
* (4) | I - UV' | / ( n ulp ) SORGTR( UPLO='L', ... )
*
* (5-8) Same as 1-4, but for SSPTRD and SOPGTR.
*
* (9) | S - Z D Z' | / ( |S| n ulp ) SSTEQR('V',...)
*
* (10) | I - ZZ' | / ( n ulp ) SSTEQR('V',...)
*
* (11) | D1 - D2 | / ( |D1| ulp ) SSTEQR('N',...)
*
* (12) | D1 - D3 | / ( |D1| ulp ) SSTERF
*
* (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
* SSTECH)
*
* For S positive definite,
*
* (14) | S - Z4 D4 Z4' | / ( |S| n ulp ) SPTEQR('V',...)
*
* (15) | I - Z4 Z4' | / ( n ulp ) SPTEQR('V',...)
*
* (16) | D4 - D5 | / ( 100 |D4| ulp ) SPTEQR('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
* SSTEBZ( 'A', 'E', ...)
*
* (18) | WA1 - D3 | / ( |D3| ulp ) SSTEBZ( 'A', 'E', ...)
*
* (19) ( max { min | WA2(i)-WA3(j) | } +
* i j
* max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
* i j
* SSTEBZ( 'I', 'E', ...)
*
* (20) | S - Y WA1 Y' | / ( |S| n ulp ) SSTEBZ, SSTEIN
*
* (21) | I - Y Y' | / ( n ulp ) SSTEBZ, SSTEIN
*
* (22) | S - Z D Z' | / ( |S| n ulp ) SSTEDC('I')
*
* (23) | I - ZZ' | / ( n ulp ) SSTEDC('I')
*
* (24) | S - Z D Z' | / ( |S| n ulp ) SSTEDC('V')
*
* (25) | I - ZZ' | / ( n ulp ) SSTEDC('V')
*
* (26) | D1 - D2 | / ( |D1| ulp ) SSTEDC('V') and
* SSTEDC('N')
*
* Test 27 is disabled at the moment because SSTEMR 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
* SSTEMR('V', 'A')
*
* (28) max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
* i
* omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
* SSTEMR('V', 'I')
*
* Tests 29 through 34 are disable at present because SSTEMR
* does not handle partial specturm requests.
*
* (29) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'I')
*
* (30) | I - ZZ' | / ( n ulp ) SSTEMR('V', 'I')
*
* (31) ( max { min | WA2(i)-WA3(j) | } +
* i j
* max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
* i j
* SSTEMR('N', 'I') vs. SSTEMR('V', 'I')
*
* (32) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'V')
*
* (33) | I - ZZ' | / ( n ulp ) SSTEMR('V', 'V')
*
* (34) ( max { min | WA2(i)-WA3(j) | } +
* i j
* max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
* i j
* SSTEMR('N', 'V') vs. SSTEMR('V', 'V')
*
* (35) | S - Z D Z' | / ( |S| n ulp ) SSTEMR('V', 'A')
*
* (36) | I - ZZ' | / ( n ulp ) SSTEMR('V', 'A')
*
* (37) ( max { min | WA2(i)-WA3(j) | } +
* i j
* max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
* i j
* SSTEMR('N', 'A') vs. SSTEMR('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 orthogonal 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 orthogonal 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 orthogonal 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) Symmetric 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,
* SCHKST 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, SCHKST
* 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 SCHKST to continue the same random number
* sequence.
*
* THRESH (input) REAL
* 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) REAL 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) REAL array of
* dimension( max(NN)*max(NN+1)/2 )
* The matrix A stored in packed format.
*
* SD (workspace/output) REAL array of
* dimension( max(NN) )
* The diagonal of the tridiagonal matrix computed by SSYTRD.
* On exit, SD and SE contain the tridiagonal form of the
* matrix in A.
*
* SE (workspace/output) REAL array of
* dimension( max(NN) )
* The off-diagonal of the tridiagonal matrix computed by
* SSYTRD. On exit, SD and SE contain the tridiagonal form of
* the matrix in A.
*
* D1 (workspace/output) REAL array of
* dimension( max(NN) )
* The eigenvalues of A, as computed by SSTEQR simlutaneously
* with Z. On exit, the eigenvalues in D1 correspond with the
* matrix in A.
*
* D2 (workspace/output) REAL array of
* dimension( max(NN) )
* The eigenvalues of A, as computed by SSTEQR if Z is not
* computed. On exit, the eigenvalues in D2 correspond with
* the matrix in A.
*
* D3 (workspace/output) REAL array of
* dimension( max(NN) )
* The eigenvalues of A, as computed by SSTERF. On exit, the
* eigenvalues in D3 correspond with the matrix in A.
*
* U (workspace/output) REAL array of
* dimension( LDU, max(NN) ).
* The orthogonal matrix computed by SSYTRD + SORGTR.
*
* 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) REAL array of
* dimension( LDU, max(NN) ).
* The Housholder vectors computed by SSYTRD 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 SSYTRD, 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 SORGTR, set those entries to
* 1 before using them, and then restore them later.)
*
* VP (workspace) REAL array of
* dimension( max(NN)*max(NN+1)/2 )
* The matrix V stored in packed format.
*
* TAU (workspace/output) REAL array of
* dimension( max(NN) )
* The Householder factors computed by SSYTRD in reducing A
* to tridiagonal form.
*
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?