1980 lines
66 KiB
Fortran
1980 lines
66 KiB
Fortran
*> \brief \b DCHKST
|
|
*
|
|
* =========== DOCUMENTATION ===========
|
|
*
|
|
* Online html documentation available at
|
|
* http://www.netlib.org/lapack/explore-html/
|
|
*
|
|
* Definition:
|
|
* ===========
|
|
*
|
|
* SUBROUTINE DCHKST( 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 )
|
|
*
|
|
* .. Scalar Arguments ..
|
|
* INTEGER INFO, LDA, LDU, LIWORK, LWORK, NOUNIT, NSIZES,
|
|
* $ NTYPES
|
|
* DOUBLE PRECISION THRESH
|
|
* ..
|
|
* .. Array Arguments ..
|
|
* LOGICAL DOTYPE( * )
|
|
* INTEGER ISEED( 4 ), IWORK( * ), NN( * )
|
|
* DOUBLE PRECISION A( LDA, * ), AP( * ), D1( * ), D2( * ),
|
|
* $ D3( * ), D4( * ), D5( * ), RESULT( * ),
|
|
* $ SD( * ), SE( * ), TAU( * ), U( LDU, * ),
|
|
* $ V( LDU, * ), VP( * ), WA1( * ), WA2( * ),
|
|
* $ WA3( * ), WORK( * ), WR( * ), Z( LDU, * )
|
|
* ..
|
|
*
|
|
*
|
|
*> \par Purpose:
|
|
* =============
|
|
*>
|
|
*> \verbatim
|
|
*>
|
|
*> DCHKST checks the symmetric eigenvalue problem routines.
|
|
*>
|
|
*> DSYTRD factors A as U S U' , where ' means transpose,
|
|
*> S is symmetric tridiagonal, and U is orthogonal.
|
|
*> DSYTRD can use either just the lower or just the upper triangle
|
|
*> of A; DCHKST 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.
|
|
*>
|
|
*> DSPTRD does the same as DSYTRD, except that A and V are stored
|
|
*> in "packed" format.
|
|
*>
|
|
*> DORGTR constructs the matrix U from the contents of V and TAU.
|
|
*>
|
|
*> DOPGTR constructs the matrix U from the contents of VP and TAU.
|
|
*>
|
|
*> DSTEQR 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.
|
|
*>
|
|
*> DSTERF computes D3, the matrix of eigenvalues, by the
|
|
*> PWK method, which does not yield eigenvectors.
|
|
*>
|
|
*> DPTEQR 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.
|
|
*>
|
|
*> 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.
|
|
*>
|
|
*> DSTEIN computes Y, the eigenvectors of S, given the
|
|
*> eigenvalues.
|
|
*>
|
|
*> DSTEDC 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 DSYTRD/DORGTR or DSPTRD/DOPGTR ('V' option). It may
|
|
*> also just compute eigenvalues ('N' option).
|
|
*>
|
|
*> DSTEMR 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). DSTEMR
|
|
*> uses the Relatively Robust Representation whenever possible.
|
|
*>
|
|
*> When DCHKST 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 ) DSYTRD( UPLO='U', ... )
|
|
*>
|
|
*> (2) | I - UV' | / ( n ulp ) DORGTR( UPLO='U', ... )
|
|
*>
|
|
*> (3) | A - V S V' | / ( |A| n ulp ) DSYTRD( UPLO='L', ... )
|
|
*>
|
|
*> (4) | I - UV' | / ( n ulp ) DORGTR( UPLO='L', ... )
|
|
*>
|
|
*> (5-8) Same as 1-4, but for DSPTRD and DOPGTR.
|
|
*>
|
|
*> (9) | S - Z D Z' | / ( |S| n ulp ) DSTEQR('V',...)
|
|
*>
|
|
*> (10) | I - ZZ' | / ( n ulp ) DSTEQR('V',...)
|
|
*>
|
|
*> (11) | D1 - D2 | / ( |D1| ulp ) DSTEQR('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 ) DPTEQR('V',...)
|
|
*>
|
|
*> (15) | I - Z4 Z4' | / ( n ulp ) DPTEQR('V',...)
|
|
*>
|
|
*> (16) | D4 - D5 | / ( 100 |D4| ulp ) DPTEQR('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, SSTEIN
|
|
*>
|
|
*> (21) | I - Y Y' | / ( n ulp ) DSTEBZ, SSTEIN
|
|
*>
|
|
*> (22) | S - Z D Z' | / ( |S| n ulp ) DSTEDC('I')
|
|
*>
|
|
*> (23) | I - ZZ' | / ( n ulp ) DSTEDC('I')
|
|
*>
|
|
*> (24) | S - Z D Z' | / ( |S| n ulp ) DSTEDC('V')
|
|
*>
|
|
*> (25) | I - ZZ' | / ( n ulp ) DSTEDC('V')
|
|
*>
|
|
*> (26) | D1 - D2 | / ( |D1| ulp ) DSTEDC('V') and
|
|
*> DSTEDC('N')
|
|
*>
|
|
*> Test 27 is disabled at the moment because DSTEMR 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
|
|
*> DSTEMR('V', 'A')
|
|
*>
|
|
*> (28) max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
|
|
*> i
|
|
*> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
|
|
*> DSTEMR('V', 'I')
|
|
*>
|
|
*> Tests 29 through 34 are disable at present because DSTEMR
|
|
*> does not handle partial specturm requests.
|
|
*>
|
|
*> (29) | S - Z D Z' | / ( |S| n ulp ) DSTEMR('V', 'I')
|
|
*>
|
|
*> (30) | I - ZZ' | / ( n ulp ) DSTEMR('V', 'I')
|
|
*>
|
|
*> (31) ( max { min | WA2(i)-WA3(j) | } +
|
|
*> i j
|
|
*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
|
|
*> i j
|
|
*> DSTEMR('N', 'I') vs. SSTEMR('V', 'I')
|
|
*>
|
|
*> (32) | S - Z D Z' | / ( |S| n ulp ) DSTEMR('V', 'V')
|
|
*>
|
|
*> (33) | I - ZZ' | / ( n ulp ) DSTEMR('V', 'V')
|
|
*>
|
|
*> (34) ( max { min | WA2(i)-WA3(j) | } +
|
|
*> i j
|
|
*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
|
|
*> i j
|
|
*> DSTEMR('N', 'V') vs. SSTEMR('V', 'V')
|
|
*>
|
|
*> (35) | S - Z D Z' | / ( |S| n ulp ) DSTEMR('V', 'A')
|
|
*>
|
|
*> (36) | I - ZZ' | / ( n ulp ) DSTEMR('V', 'A')
|
|
*>
|
|
*> (37) ( max { min | WA2(i)-WA3(j) | } +
|
|
*> i j
|
|
*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
|
|
*> i j
|
|
*> DSTEMR('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.
|
|
*> \endverbatim
|
|
*
|
|
* Arguments:
|
|
* ==========
|
|
*
|
|
*> \param[in] NSIZES
|
|
*> \verbatim
|
|
*> NSIZES is INTEGER
|
|
*> The number of sizes of matrices to use. If it is zero,
|
|
*> DCHKST does nothing. It must be at least zero.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] NN
|
|
*> \verbatim
|
|
*> NN is 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.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] NTYPES
|
|
*> \verbatim
|
|
*> NTYPES is INTEGER
|
|
*> The number of elements in DOTYPE. If it is zero, DCHKST
|
|
*> 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. .
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] DOTYPE
|
|
*> \verbatim
|
|
*> DOTYPE is 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.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in,out] ISEED
|
|
*> \verbatim
|
|
*> ISEED is 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 DCHKST to continue the same random number
|
|
*> sequence.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] THRESH
|
|
*> \verbatim
|
|
*> THRESH is 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.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] NOUNIT
|
|
*> \verbatim
|
|
*> NOUNIT is INTEGER
|
|
*> The FORTRAN unit number for printing out error messages
|
|
*> (e.g., if a routine returns IINFO not equal to 0.)
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in,out] A
|
|
*> \verbatim
|
|
*> A is DOUBLE PRECISION 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.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] LDA
|
|
*> \verbatim
|
|
*> LDA is INTEGER
|
|
*> The leading dimension of A. It must be at
|
|
*> least 1 and at least max( NN ).
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] AP
|
|
*> \verbatim
|
|
*> AP is DOUBLE PRECISION array of
|
|
*> dimension( max(NN)*max(NN+1)/2 )
|
|
*> The matrix A stored in packed format.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] SD
|
|
*> \verbatim
|
|
*> SD is DOUBLE PRECISION array of
|
|
*> dimension( max(NN) )
|
|
*> The diagonal of the tridiagonal matrix computed by DSYTRD.
|
|
*> On exit, SD and SE contain the tridiagonal form of the
|
|
*> matrix in A.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] SE
|
|
*> \verbatim
|
|
*> SE is DOUBLE PRECISION array of
|
|
*> dimension( max(NN) )
|
|
*> The off-diagonal of the tridiagonal matrix computed by
|
|
*> DSYTRD. On exit, SD and SE contain the tridiagonal form of
|
|
*> the matrix in A.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] D1
|
|
*> \verbatim
|
|
*> D1 is DOUBLE PRECISION array of
|
|
*> dimension( max(NN) )
|
|
*> The eigenvalues of A, as computed by DSTEQR simlutaneously
|
|
*> with Z. On exit, the eigenvalues in D1 correspond with the
|
|
*> matrix in A.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] D2
|
|
*> \verbatim
|
|
*> D2 is DOUBLE PRECISION array of
|
|
*> dimension( max(NN) )
|
|
*> The eigenvalues of A, as computed by DSTEQR if Z is not
|
|
*> computed. On exit, the eigenvalues in D2 correspond with
|
|
*> the matrix in A.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] D3
|
|
*> \verbatim
|
|
*> D3 is 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.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] D4
|
|
*> \verbatim
|
|
*> D4 is DOUBLE PRECISION array of
|
|
*> dimension( max(NN) )
|
|
*> The eigenvalues of A, as computed by DPTEQR(V).
|
|
*> DPTEQR factors S as Z4 D4 Z4*
|
|
*> On exit, the eigenvalues in D4 correspond with the matrix in A.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] D5
|
|
*> \verbatim
|
|
*> D5 is DOUBLE PRECISION array of
|
|
*> dimension( max(NN) )
|
|
*> The eigenvalues of A, as computed by DPTEQR(N)
|
|
*> when Z is not computed. On exit, the
|
|
*> eigenvalues in D4 correspond with the matrix in A.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] WA1
|
|
*> \verbatim
|
|
*> WA1 is DOUBLE PRECISION array of
|
|
*> dimension( max(NN) )
|
|
*> All eigenvalues of A, computed to high
|
|
*> absolute accuracy, with different range options.
|
|
*> as computed by DSTEBZ.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] WA2
|
|
*> \verbatim
|
|
*> WA2 is DOUBLE PRECISION array of
|
|
*> dimension( max(NN) )
|
|
*> Selected eigenvalues of A, computed to high
|
|
*> absolute accuracy, with different range options.
|
|
*> as computed by DSTEBZ.
|
|
*> Choose random values for IL and IU, and ask for the
|
|
*> IL-th through IU-th eigenvalues.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] WA3
|
|
*> \verbatim
|
|
*> WA3 is DOUBLE PRECISION array of
|
|
*> dimension( max(NN) )
|
|
*> Selected eigenvalues of A, computed to high
|
|
*> absolute accuracy, with different range options.
|
|
*> as computed by DSTEBZ.
|
|
*> Determine the values VL and VU of the IL-th and IU-th
|
|
*> eigenvalues and ask for all eigenvalues in this range.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] WR
|
|
*> \verbatim
|
|
*> WR is DOUBLE PRECISION array of
|
|
*> dimension( max(NN) )
|
|
*> All eigenvalues of A, computed to high
|
|
*> absolute accuracy, with different options.
|
|
*> as computed by DSTEBZ.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] U
|
|
*> \verbatim
|
|
*> U is DOUBLE PRECISION array of
|
|
*> dimension( LDU, max(NN) ).
|
|
*> The orthogonal matrix computed by DSYTRD + DORGTR.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] LDU
|
|
*> \verbatim
|
|
*> LDU is INTEGER
|
|
*> The leading dimension of U, Z, and V. It must be at least 1
|
|
*> and at least max( NN ).
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] V
|
|
*> \verbatim
|
|
*> V is DOUBLE PRECISION array of
|
|
*> dimension( LDU, max(NN) ).
|
|
*> The Housholder vectors computed by DSYTRD 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 DSYTRD, 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 DORGTR, set those entries to
|
|
*> 1 before using them, and then restore them later.)
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] VP
|
|
*> \verbatim
|
|
*> VP is DOUBLE PRECISION array of
|
|
*> dimension( max(NN)*max(NN+1)/2 )
|
|
*> The matrix V stored in packed format.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] TAU
|
|
*> \verbatim
|
|
*> TAU is DOUBLE PRECISION array of
|
|
*> dimension( max(NN) )
|
|
*> The Householder factors computed by DSYTRD in reducing A
|
|
*> to tridiagonal form.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] Z
|
|
*> \verbatim
|
|
*> Z is DOUBLE PRECISION array of
|
|
*> dimension( LDU, max(NN) ).
|
|
*> The orthogonal matrix of eigenvectors computed by DSTEQR,
|
|
*> DPTEQR, and DSTEIN.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] WORK
|
|
*> \verbatim
|
|
*> WORK is DOUBLE PRECISION array of
|
|
*> dimension( LWORK )
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] LWORK
|
|
*> \verbatim
|
|
*> LWORK is INTEGER
|
|
*> The number of entries in WORK. This must be at least
|
|
*> 1 + 4 * Nmax + 2 * Nmax * lg Nmax + 3 * Nmax**2
|
|
*> where Nmax = max( NN(j), 2 ) and lg = log base 2.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] IWORK
|
|
*> \verbatim
|
|
*> IWORK is INTEGER array,
|
|
*> Workspace.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] LIWORK
|
|
*> \verbatim
|
|
*> LIWORK is INTEGER
|
|
*> The number of entries in IWORK. This must be at least
|
|
*> 6 + 6*Nmax + 5 * Nmax * lg Nmax
|
|
*> where Nmax = max( NN(j), 2 ) and lg = log base 2.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] RESULT
|
|
*> \verbatim
|
|
*> RESULT is DOUBLE PRECISION array, dimension (26)
|
|
*> The values computed by the tests described above.
|
|
*> The values are currently limited to 1/ulp, to avoid
|
|
*> overflow.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] INFO
|
|
*> \verbatim
|
|
*> INFO is INTEGER
|
|
*> If 0, then everything ran OK.
|
|
*> -1: NSIZES < 0
|
|
*> -2: Some NN(j) < 0
|
|
*> -3: NTYPES < 0
|
|
*> -5: THRESH < 0
|
|
*> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
|
|
*> -23: LDU < 1 or LDU < NMAX.
|
|
*> -29: LWORK too small.
|
|
*> If DLATMR, SLATMS, DSYTRD, DORGTR, DSTEQR, SSTERF,
|
|
*> or DORMC2 returns an error code, the
|
|
*> absolute value of it is returned.
|
|
*>
|
|
*>-----------------------------------------------------------------------
|
|
*>
|
|
*> Some Local Variables and Parameters:
|
|
*> ---- ----- --------- --- ----------
|
|
*> ZERO, ONE Real 0 and 1.
|
|
*> MAXTYP The number of types defined.
|
|
*> NTEST The number of tests performed, or which can
|
|
*> be performed so far, for the current matrix.
|
|
*> NTESTT The total number of tests performed so far.
|
|
*> NBLOCK Blocksize as returned by ENVIR.
|
|
*> NMAX Largest value in NN.
|
|
*> NMATS The number of matrices generated so far.
|
|
*> NERRS The number of tests which have exceeded THRESH
|
|
*> so far.
|
|
*> COND, IMODE Values to be passed to the matrix generators.
|
|
*> ANORM Norm of A; passed to matrix generators.
|
|
*>
|
|
*> OVFL, UNFL Overflow and underflow thresholds.
|
|
*> ULP, ULPINV Finest relative precision and its inverse.
|
|
*> RTOVFL, RTUNFL Square roots of the previous 2 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) )
|
|
*> \endverbatim
|
|
*
|
|
* Authors:
|
|
* ========
|
|
*
|
|
*> \author Univ. of Tennessee
|
|
*> \author Univ. of California Berkeley
|
|
*> \author Univ. of Colorado Denver
|
|
*> \author NAG Ltd.
|
|
*
|
|
*> \date December 2016
|
|
*
|
|
*> \ingroup double_eig
|
|
*
|
|
* =====================================================================
|
|
SUBROUTINE DCHKST( 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 )
|
|
*
|
|
* -- LAPACK test routine (version 3.7.0) --
|
|
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
|
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
|
* December 2016
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INFO, LDA, LDU, LIWORK, LWORK, NOUNIT, NSIZES,
|
|
$ NTYPES
|
|
DOUBLE PRECISION THRESH
|
|
* ..
|
|
* .. Array Arguments ..
|
|
LOGICAL DOTYPE( * )
|
|
INTEGER ISEED( 4 ), IWORK( * ), NN( * )
|
|
DOUBLE PRECISION A( LDA, * ), AP( * ), D1( * ), D2( * ),
|
|
$ D3( * ), D4( * ), D5( * ), RESULT( * ),
|
|
$ SD( * ), SE( * ), TAU( * ), U( LDU, * ),
|
|
$ V( LDU, * ), VP( * ), WA1( * ), WA2( * ),
|
|
$ WA3( * ), WORK( * ), WR( * ), Z( LDU, * )
|
|
* ..
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ZERO, ONE, TWO, EIGHT, TEN, HUN
|
|
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
|
|
$ EIGHT = 8.0D0, TEN = 10.0D0, HUN = 100.0D0 )
|
|
DOUBLE PRECISION HALF
|
|
PARAMETER ( HALF = ONE / TWO )
|
|
INTEGER MAXTYP
|
|
PARAMETER ( MAXTYP = 21 )
|
|
LOGICAL SRANGE
|
|
PARAMETER ( SRANGE = .FALSE. )
|
|
LOGICAL SREL
|
|
PARAMETER ( SREL = .FALSE. )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL BADNN, TRYRAC
|
|
INTEGER I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, J, JC,
|
|
$ JR, JSIZE, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC,
|
|
$ M, M2, M3, MTYPES, N, NAP, NBLOCK, NERRS,
|
|
$ NMATS, NMAX, NSPLIT, NTEST, NTESTT
|
|
DOUBLE PRECISION ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL,
|
|
$ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP,
|
|
$ ULPINV, UNFL, VL, VU
|
|
* ..
|
|
* .. Local Arrays ..
|
|
INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
|
|
$ KMAGN( MAXTYP ), KMODE( MAXTYP ),
|
|
$ KTYPE( MAXTYP )
|
|
DOUBLE PRECISION DUMMA( 1 )
|
|
* ..
|
|
* .. External Functions ..
|
|
INTEGER ILAENV
|
|
DOUBLE PRECISION DLAMCH, DLARND, DSXT1
|
|
EXTERNAL ILAENV, DLAMCH, DLARND, DSXT1
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DCOPY, DLABAD, DLACPY, DLASET, DLASUM, DLATMR,
|
|
$ DLATMS, DOPGTR, DORGTR, DPTEQR, DSPT21, DSPTRD,
|
|
$ DSTEBZ, DSTECH, DSTEDC, DSTEMR, DSTEIN, DSTEQR,
|
|
$ DSTERF, DSTT21, DSTT22, DSYT21, DSYTRD, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC ABS, DBLE, INT, LOG, MAX, MIN, SQRT
|
|
* ..
|
|
* .. Data statements ..
|
|
DATA KTYPE / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8,
|
|
$ 8, 8, 9, 9, 9, 9, 9, 10 /
|
|
DATA KMAGN / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
|
|
$ 2, 3, 1, 1, 1, 2, 3, 1 /
|
|
DATA KMODE / 0, 0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
|
|
$ 0, 0, 4, 3, 1, 4, 4, 3 /
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Keep ftnchek happy
|
|
IDUMMA( 1 ) = 1
|
|
*
|
|
* Check for errors
|
|
*
|
|
NTESTT = 0
|
|
INFO = 0
|
|
*
|
|
* Important constants
|
|
*
|
|
BADNN = .FALSE.
|
|
TRYRAC = .TRUE.
|
|
NMAX = 1
|
|
DO 10 J = 1, NSIZES
|
|
NMAX = MAX( NMAX, NN( J ) )
|
|
IF( NN( J ).LT.0 )
|
|
$ BADNN = .TRUE.
|
|
10 CONTINUE
|
|
*
|
|
NBLOCK = ILAENV( 1, 'DSYTRD', 'L', NMAX, -1, -1, -1 )
|
|
NBLOCK = MIN( NMAX, MAX( 1, NBLOCK ) )
|
|
*
|
|
* 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 = -23
|
|
ELSE IF( 2*MAX( 2, NMAX )**2.GT.LWORK ) THEN
|
|
INFO = -29
|
|
END IF
|
|
*
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DCHKST', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
* More Important constants
|
|
*
|
|
UNFL = DLAMCH( 'Safe minimum' )
|
|
OVFL = ONE / UNFL
|
|
CALL DLABAD( UNFL, OVFL )
|
|
ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
|
|
ULPINV = ONE / ULP
|
|
LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) )
|
|
RTUNFL = SQRT( UNFL )
|
|
RTOVFL = SQRT( OVFL )
|
|
*
|
|
* Loop over sizes, types
|
|
*
|
|
DO 20 I = 1, 4
|
|
ISEED2( I ) = ISEED( I )
|
|
20 CONTINUE
|
|
NERRS = 0
|
|
NMATS = 0
|
|
*
|
|
DO 310 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
|
|
LIWEDC = 6 + 6*N + 5*N*LGN
|
|
ELSE
|
|
LWEDC = 8
|
|
LIWEDC = 12
|
|
END IF
|
|
NAP = ( N*( N+1 ) ) / 2
|
|
ANINV = ONE / DBLE( MAX( 1, N ) )
|
|
*
|
|
IF( NSIZES.NE.1 ) THEN
|
|
MTYPES = MIN( MAXTYP, NTYPES )
|
|
ELSE
|
|
MTYPES = MIN( MAXTYP+1, NTYPES )
|
|
END IF
|
|
*
|
|
DO 300 JTYPE = 1, MTYPES
|
|
IF( .NOT.DOTYPE( JTYPE ) )
|
|
$ GO TO 300
|
|
NMATS = NMATS + 1
|
|
NTEST = 0
|
|
*
|
|
DO 30 J = 1, 4
|
|
IOLDSD( J ) = ISEED( J )
|
|
30 CONTINUE
|
|
*
|
|
* 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 positive definite
|
|
* =10 diagonally dominant tridiagonal
|
|
*
|
|
IF( MTYPES.GT.MAXTYP )
|
|
$ GO TO 100
|
|
*
|
|
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
|
|
IF( JTYPE.LE.15 ) THEN
|
|
COND = ULPINV
|
|
ELSE
|
|
COND = ULPINV*ANINV / TEN
|
|
END IF
|
|
*
|
|
* Special Matrices -- Identity & Jordan block
|
|
*
|
|
* Zero
|
|
*
|
|
IF( ITYPE.EQ.1 ) THEN
|
|
IINFO = 0
|
|
*
|
|
ELSE IF( ITYPE.EQ.2 ) THEN
|
|
*
|
|
* Identity
|
|
*
|
|
DO 80 JC = 1, N
|
|
A( JC, JC ) = 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
|
|
*
|
|
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
|
|
*
|
|
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
|
|
*
|
|
* Positive definite, eigenvalues specified.
|
|
*
|
|
CALL DLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND,
|
|
$ ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
|
|
$ IINFO )
|
|
*
|
|
ELSE IF( ITYPE.EQ.10 ) THEN
|
|
*
|
|
* Positive definite tridiagonal, eigenvalues specified.
|
|
*
|
|
CALL DLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND,
|
|
$ ANORM, 1, 1, 'N', A, LDA, WORK( N+1 ),
|
|
$ IINFO )
|
|
DO 90 I = 2, N
|
|
TEMP1 = ABS( A( I-1, I ) ) /
|
|
$ SQRT( ABS( A( I-1, I-1 )*A( I, I ) ) )
|
|
IF( TEMP1.GT.HALF ) THEN
|
|
A( I-1, I ) = HALF*SQRT( ABS( A( I-1, I-1 )*A( I,
|
|
$ I ) ) )
|
|
A( I, I-1 ) = A( I-1, I )
|
|
END IF
|
|
90 CONTINUE
|
|
*
|
|
ELSE
|
|
*
|
|
IINFO = 1
|
|
END IF
|
|
*
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
100 CONTINUE
|
|
*
|
|
* Call DSYTRD and DORGTR to compute S and U from
|
|
* upper triangle.
|
|
*
|
|
CALL DLACPY( 'U', N, N, A, LDA, V, LDU )
|
|
*
|
|
NTEST = 1
|
|
CALL DSYTRD( 'U', N, V, LDU, SD, SE, TAU, WORK, LWORK,
|
|
$ IINFO )
|
|
*
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSYTRD(U)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 1 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
CALL DLACPY( 'U', N, N, V, LDU, U, LDU )
|
|
*
|
|
NTEST = 2
|
|
CALL DORGTR( 'U', N, U, LDU, TAU, WORK, LWORK, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DORGTR(U)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 2 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do tests 1 and 2
|
|
*
|
|
CALL DSYT21( 2, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V,
|
|
$ LDU, TAU, WORK, RESULT( 1 ) )
|
|
CALL DSYT21( 3, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V,
|
|
$ LDU, TAU, WORK, RESULT( 2 ) )
|
|
*
|
|
* Call DSYTRD and DORGTR to compute S and U from
|
|
* lower triangle, do tests.
|
|
*
|
|
CALL DLACPY( 'L', N, N, A, LDA, V, LDU )
|
|
*
|
|
NTEST = 3
|
|
CALL DSYTRD( 'L', N, V, LDU, SD, SE, TAU, WORK, LWORK,
|
|
$ IINFO )
|
|
*
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSYTRD(L)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 3 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
CALL DLACPY( 'L', N, N, V, LDU, U, LDU )
|
|
*
|
|
NTEST = 4
|
|
CALL DORGTR( 'L', N, U, LDU, TAU, WORK, LWORK, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DORGTR(L)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 4 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
CALL DSYT21( 2, 'Lower', N, 1, A, LDA, SD, SE, U, LDU, V,
|
|
$ LDU, TAU, WORK, RESULT( 3 ) )
|
|
CALL DSYT21( 3, 'Lower', N, 1, A, LDA, SD, SE, U, LDU, V,
|
|
$ LDU, TAU, WORK, RESULT( 4 ) )
|
|
*
|
|
* Store the upper triangle of A in AP
|
|
*
|
|
I = 0
|
|
DO 120 JC = 1, N
|
|
DO 110 JR = 1, JC
|
|
I = I + 1
|
|
AP( I ) = A( JR, JC )
|
|
110 CONTINUE
|
|
120 CONTINUE
|
|
*
|
|
* Call DSPTRD and DOPGTR to compute S and U from AP
|
|
*
|
|
CALL DCOPY( NAP, AP, 1, VP, 1 )
|
|
*
|
|
NTEST = 5
|
|
CALL DSPTRD( 'U', N, VP, SD, SE, TAU, IINFO )
|
|
*
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSPTRD(U)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 5 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
NTEST = 6
|
|
CALL DOPGTR( 'U', N, VP, TAU, U, LDU, WORK, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DOPGTR(U)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 6 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do tests 5 and 6
|
|
*
|
|
CALL DSPT21( 2, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU,
|
|
$ WORK, RESULT( 5 ) )
|
|
CALL DSPT21( 3, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU,
|
|
$ WORK, RESULT( 6 ) )
|
|
*
|
|
* Store the lower triangle of A in AP
|
|
*
|
|
I = 0
|
|
DO 140 JC = 1, N
|
|
DO 130 JR = JC, N
|
|
I = I + 1
|
|
AP( I ) = A( JR, JC )
|
|
130 CONTINUE
|
|
140 CONTINUE
|
|
*
|
|
* Call DSPTRD and DOPGTR to compute S and U from AP
|
|
*
|
|
CALL DCOPY( NAP, AP, 1, VP, 1 )
|
|
*
|
|
NTEST = 7
|
|
CALL DSPTRD( 'L', N, VP, SD, SE, TAU, IINFO )
|
|
*
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSPTRD(L)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 7 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
NTEST = 8
|
|
CALL DOPGTR( 'L', N, VP, TAU, U, LDU, WORK, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DOPGTR(L)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 8 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
CALL DSPT21( 2, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU,
|
|
$ WORK, RESULT( 7 ) )
|
|
CALL DSPT21( 3, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU,
|
|
$ WORK, RESULT( 8 ) )
|
|
*
|
|
* Call DSTEQR to compute D1, D2, and Z, do tests.
|
|
*
|
|
* Compute D1 and Z
|
|
*
|
|
CALL DCOPY( N, SD, 1, D1, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
|
|
*
|
|
NTEST = 9
|
|
CALL DSTEQR( 'V', N, D1, WORK, Z, LDU, WORK( N+1 ), IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEQR(V)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 9 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Compute D2
|
|
*
|
|
CALL DCOPY( N, SD, 1, D2, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
*
|
|
NTEST = 11
|
|
CALL DSTEQR( 'N', N, D2, WORK, WORK( N+1 ), LDU,
|
|
$ WORK( N+1 ), IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEQR(N)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 11 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Compute D3 (using PWK method)
|
|
*
|
|
CALL DCOPY( N, SD, 1, D3, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
*
|
|
NTEST = 12
|
|
CALL DSTERF( N, D3, WORK, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTERF', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 12 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do Tests 9 and 10
|
|
*
|
|
CALL DSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
|
|
$ RESULT( 9 ) )
|
|
*
|
|
* Do Tests 11 and 12
|
|
*
|
|
TEMP1 = ZERO
|
|
TEMP2 = ZERO
|
|
TEMP3 = ZERO
|
|
TEMP4 = ZERO
|
|
*
|
|
DO 150 J = 1, N
|
|
TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
|
|
TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
|
|
TEMP3 = MAX( TEMP3, ABS( D1( J ) ), ABS( D3( J ) ) )
|
|
TEMP4 = MAX( TEMP4, ABS( D1( J )-D3( J ) ) )
|
|
150 CONTINUE
|
|
*
|
|
RESULT( 11 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
|
|
RESULT( 12 ) = TEMP4 / MAX( UNFL, ULP*MAX( TEMP3, TEMP4 ) )
|
|
*
|
|
* Do Test 13 -- Sturm Sequence Test of Eigenvalues
|
|
* Go up by factors of two until it succeeds
|
|
*
|
|
NTEST = 13
|
|
TEMP1 = THRESH*( HALF-ULP )
|
|
*
|
|
DO 160 J = 0, LOG2UI
|
|
CALL DSTECH( N, SD, SE, D1, TEMP1, WORK, IINFO )
|
|
IF( IINFO.EQ.0 )
|
|
$ GO TO 170
|
|
TEMP1 = TEMP1*TWO
|
|
160 CONTINUE
|
|
*
|
|
170 CONTINUE
|
|
RESULT( 13 ) = TEMP1
|
|
*
|
|
* For positive definite matrices ( JTYPE.GT.15 ) call DPTEQR
|
|
* and do tests 14, 15, and 16 .
|
|
*
|
|
IF( JTYPE.GT.15 ) THEN
|
|
*
|
|
* Compute D4 and Z4
|
|
*
|
|
CALL DCOPY( N, SD, 1, D4, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
|
|
*
|
|
NTEST = 14
|
|
CALL DPTEQR( 'V', N, D4, WORK, Z, LDU, WORK( N+1 ),
|
|
$ IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DPTEQR(V)', IINFO, N,
|
|
$ JTYPE, IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 14 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do Tests 14 and 15
|
|
*
|
|
CALL DSTT21( N, 0, SD, SE, D4, DUMMA, Z, LDU, WORK,
|
|
$ RESULT( 14 ) )
|
|
*
|
|
* Compute D5
|
|
*
|
|
CALL DCOPY( N, SD, 1, D5, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
*
|
|
NTEST = 16
|
|
CALL DPTEQR( 'N', N, D5, WORK, Z, LDU, WORK( N+1 ),
|
|
$ IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DPTEQR(N)', IINFO, N,
|
|
$ JTYPE, IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 16 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do Test 16
|
|
*
|
|
TEMP1 = ZERO
|
|
TEMP2 = ZERO
|
|
DO 180 J = 1, N
|
|
TEMP1 = MAX( TEMP1, ABS( D4( J ) ), ABS( D5( J ) ) )
|
|
TEMP2 = MAX( TEMP2, ABS( D4( J )-D5( J ) ) )
|
|
180 CONTINUE
|
|
*
|
|
RESULT( 16 ) = TEMP2 / MAX( UNFL,
|
|
$ HUN*ULP*MAX( TEMP1, TEMP2 ) )
|
|
ELSE
|
|
RESULT( 14 ) = ZERO
|
|
RESULT( 15 ) = ZERO
|
|
RESULT( 16 ) = ZERO
|
|
END IF
|
|
*
|
|
* Call DSTEBZ with different options and do tests 17-18.
|
|
*
|
|
* If S is positive definite and diagonally dominant,
|
|
* ask for all eigenvalues with high relative accuracy.
|
|
*
|
|
VL = ZERO
|
|
VU = ZERO
|
|
IL = 0
|
|
IU = 0
|
|
IF( JTYPE.EQ.21 ) THEN
|
|
NTEST = 17
|
|
ABSTOL = UNFL + UNFL
|
|
CALL DSTEBZ( 'A', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE,
|
|
$ M, NSPLIT, WR, IWORK( 1 ), IWORK( N+1 ),
|
|
$ WORK, IWORK( 2*N+1 ), IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(A,rel)', IINFO, N,
|
|
$ JTYPE, IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 17 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do test 17
|
|
*
|
|
TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) /
|
|
$ ( ONE-HALF )**4
|
|
*
|
|
TEMP1 = ZERO
|
|
DO 190 J = 1, N
|
|
TEMP1 = MAX( TEMP1, ABS( D4( J )-WR( N-J+1 ) ) /
|
|
$ ( ABSTOL+ABS( D4( J ) ) ) )
|
|
190 CONTINUE
|
|
*
|
|
RESULT( 17 ) = TEMP1 / TEMP2
|
|
ELSE
|
|
RESULT( 17 ) = ZERO
|
|
END IF
|
|
*
|
|
* Now ask for all eigenvalues with high absolute accuracy.
|
|
*
|
|
NTEST = 18
|
|
ABSTOL = UNFL + UNFL
|
|
CALL DSTEBZ( 'A', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE, M,
|
|
$ NSPLIT, WA1, IWORK( 1 ), IWORK( N+1 ), WORK,
|
|
$ IWORK( 2*N+1 ), IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(A)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 18 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do test 18
|
|
*
|
|
TEMP1 = ZERO
|
|
TEMP2 = ZERO
|
|
DO 200 J = 1, N
|
|
TEMP1 = MAX( TEMP1, ABS( D3( J ) ), ABS( WA1( J ) ) )
|
|
TEMP2 = MAX( TEMP2, ABS( D3( J )-WA1( J ) ) )
|
|
200 CONTINUE
|
|
*
|
|
RESULT( 18 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
|
|
*
|
|
* Choose random values for IL and IU, and ask for the
|
|
* IL-th through IU-th eigenvalues.
|
|
*
|
|
NTEST = 19
|
|
IF( N.LE.1 ) THEN
|
|
IL = 1
|
|
IU = N
|
|
ELSE
|
|
IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
|
|
IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
|
|
IF( IU.LT.IL ) THEN
|
|
ITEMP = IU
|
|
IU = IL
|
|
IL = ITEMP
|
|
END IF
|
|
END IF
|
|
*
|
|
CALL DSTEBZ( 'I', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE,
|
|
$ M2, NSPLIT, WA2, IWORK( 1 ), IWORK( N+1 ),
|
|
$ WORK, IWORK( 2*N+1 ), IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(I)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 19 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Determine the values VL and VU of the IL-th and IU-th
|
|
* eigenvalues and ask for all eigenvalues in this range.
|
|
*
|
|
IF( N.GT.0 ) THEN
|
|
IF( IL.NE.1 ) THEN
|
|
VL = WA1( IL ) - MAX( HALF*( WA1( IL )-WA1( IL-1 ) ),
|
|
$ ULP*ANORM, TWO*RTUNFL )
|
|
ELSE
|
|
VL = WA1( 1 ) - MAX( HALF*( WA1( N )-WA1( 1 ) ),
|
|
$ ULP*ANORM, TWO*RTUNFL )
|
|
END IF
|
|
IF( IU.NE.N ) THEN
|
|
VU = WA1( IU ) + MAX( HALF*( WA1( IU+1 )-WA1( IU ) ),
|
|
$ ULP*ANORM, TWO*RTUNFL )
|
|
ELSE
|
|
VU = WA1( N ) + MAX( HALF*( WA1( N )-WA1( 1 ) ),
|
|
$ ULP*ANORM, TWO*RTUNFL )
|
|
END IF
|
|
ELSE
|
|
VL = ZERO
|
|
VU = ONE
|
|
END IF
|
|
*
|
|
CALL DSTEBZ( 'V', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE,
|
|
$ M3, NSPLIT, WA3, IWORK( 1 ), IWORK( N+1 ),
|
|
$ WORK, IWORK( 2*N+1 ), IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(V)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 19 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
IF( M3.EQ.0 .AND. N.NE.0 ) THEN
|
|
RESULT( 19 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
*
|
|
* Do test 19
|
|
*
|
|
TEMP1 = DSXT1( 1, WA2, M2, WA3, M3, ABSTOL, ULP, UNFL )
|
|
TEMP2 = DSXT1( 1, WA3, M3, WA2, M2, ABSTOL, ULP, UNFL )
|
|
IF( N.GT.0 ) THEN
|
|
TEMP3 = MAX( ABS( WA1( N ) ), ABS( WA1( 1 ) ) )
|
|
ELSE
|
|
TEMP3 = ZERO
|
|
END IF
|
|
*
|
|
RESULT( 19 ) = ( TEMP1+TEMP2 ) / MAX( UNFL, TEMP3*ULP )
|
|
*
|
|
* Call DSTEIN to compute eigenvectors corresponding to
|
|
* eigenvalues in WA1. (First call DSTEBZ again, to make sure
|
|
* it returns these eigenvalues in the correct order.)
|
|
*
|
|
NTEST = 21
|
|
CALL DSTEBZ( 'A', 'B', N, VL, VU, IL, IU, ABSTOL, SD, SE, M,
|
|
$ NSPLIT, WA1, IWORK( 1 ), IWORK( N+1 ), WORK,
|
|
$ IWORK( 2*N+1 ), IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(A,B)', IINFO, N,
|
|
$ JTYPE, IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 20 ) = ULPINV
|
|
RESULT( 21 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
CALL DSTEIN( N, SD, SE, M, WA1, IWORK( 1 ), IWORK( N+1 ), Z,
|
|
$ LDU, WORK, IWORK( 2*N+1 ), IWORK( 3*N+1 ),
|
|
$ IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEIN', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 20 ) = ULPINV
|
|
RESULT( 21 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do tests 20 and 21
|
|
*
|
|
CALL DSTT21( N, 0, SD, SE, WA1, DUMMA, Z, LDU, WORK,
|
|
$ RESULT( 20 ) )
|
|
*
|
|
* Call DSTEDC(I) to compute D1 and Z, do tests.
|
|
*
|
|
* Compute D1 and Z
|
|
*
|
|
CALL DCOPY( N, SD, 1, D1, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
|
|
*
|
|
NTEST = 22
|
|
CALL DSTEDC( 'I', N, D1, WORK, Z, LDU, WORK( N+1 ), LWEDC-N,
|
|
$ IWORK, LIWEDC, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEDC(I)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 22 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do Tests 22 and 23
|
|
*
|
|
CALL DSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
|
|
$ RESULT( 22 ) )
|
|
*
|
|
* Call DSTEDC(V) to compute D1 and Z, do tests.
|
|
*
|
|
* Compute D1 and Z
|
|
*
|
|
CALL DCOPY( N, SD, 1, D1, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
|
|
*
|
|
NTEST = 24
|
|
CALL DSTEDC( 'V', N, D1, WORK, Z, LDU, WORK( N+1 ), LWEDC-N,
|
|
$ IWORK, LIWEDC, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEDC(V)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 24 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do Tests 24 and 25
|
|
*
|
|
CALL DSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
|
|
$ RESULT( 24 ) )
|
|
*
|
|
* Call DSTEDC(N) to compute D2, do tests.
|
|
*
|
|
* Compute D2
|
|
*
|
|
CALL DCOPY( N, SD, 1, D2, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
|
|
*
|
|
NTEST = 26
|
|
CALL DSTEDC( 'N', N, D2, WORK, Z, LDU, WORK( N+1 ), LWEDC-N,
|
|
$ IWORK, LIWEDC, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEDC(N)', IINFO, N, JTYPE,
|
|
$ IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 26 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do Test 26
|
|
*
|
|
TEMP1 = ZERO
|
|
TEMP2 = ZERO
|
|
*
|
|
DO 210 J = 1, N
|
|
TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
|
|
TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
|
|
210 CONTINUE
|
|
*
|
|
RESULT( 26 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
|
|
*
|
|
* Only test DSTEMR if IEEE compliant
|
|
*
|
|
IF( ILAENV( 10, 'DSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND.
|
|
$ ILAENV( 11, 'DSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN
|
|
*
|
|
* Call DSTEMR, do test 27 (relative eigenvalue accuracy)
|
|
*
|
|
* If S is positive definite and diagonally dominant,
|
|
* ask for all eigenvalues with high relative accuracy.
|
|
*
|
|
VL = ZERO
|
|
VU = ZERO
|
|
IL = 0
|
|
IU = 0
|
|
IF( JTYPE.EQ.21 .AND. SREL ) THEN
|
|
NTEST = 27
|
|
ABSTOL = UNFL + UNFL
|
|
CALL DSTEMR( 'V', 'A', N, SD, SE, VL, VU, IL, IU,
|
|
$ M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC,
|
|
$ WORK, LWORK, IWORK( 2*N+1 ), LWORK-2*N,
|
|
$ IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,A,rel)',
|
|
$ IINFO, N, JTYPE, IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 27 ) = ULPINV
|
|
GO TO 270
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do test 27
|
|
*
|
|
TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) /
|
|
$ ( ONE-HALF )**4
|
|
*
|
|
TEMP1 = ZERO
|
|
DO 220 J = 1, N
|
|
TEMP1 = MAX( TEMP1, ABS( D4( J )-WR( N-J+1 ) ) /
|
|
$ ( ABSTOL+ABS( D4( J ) ) ) )
|
|
220 CONTINUE
|
|
*
|
|
RESULT( 27 ) = TEMP1 / TEMP2
|
|
*
|
|
IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
|
|
IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
|
|
IF( IU.LT.IL ) THEN
|
|
ITEMP = IU
|
|
IU = IL
|
|
IL = ITEMP
|
|
END IF
|
|
*
|
|
IF( SRANGE ) THEN
|
|
NTEST = 28
|
|
ABSTOL = UNFL + UNFL
|
|
CALL DSTEMR( 'V', 'I', N, SD, SE, VL, VU, IL, IU,
|
|
$ M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC,
|
|
$ WORK, LWORK, IWORK( 2*N+1 ),
|
|
$ LWORK-2*N, IINFO )
|
|
*
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,I,rel)',
|
|
$ IINFO, N, JTYPE, IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 28 ) = ULPINV
|
|
GO TO 270
|
|
END IF
|
|
END IF
|
|
*
|
|
*
|
|
* Do test 28
|
|
*
|
|
TEMP2 = TWO*( TWO*N-ONE )*ULP*
|
|
$ ( ONE+EIGHT*HALF**2 ) / ( ONE-HALF )**4
|
|
*
|
|
TEMP1 = ZERO
|
|
DO 230 J = IL, IU
|
|
TEMP1 = MAX( TEMP1, ABS( WR( J-IL+1 )-D4( N-J+
|
|
$ 1 ) ) / ( ABSTOL+ABS( WR( J-IL+1 ) ) ) )
|
|
230 CONTINUE
|
|
*
|
|
RESULT( 28 ) = TEMP1 / TEMP2
|
|
ELSE
|
|
RESULT( 28 ) = ZERO
|
|
END IF
|
|
ELSE
|
|
RESULT( 27 ) = ZERO
|
|
RESULT( 28 ) = ZERO
|
|
END IF
|
|
*
|
|
* Call DSTEMR(V,I) to compute D1 and Z, do tests.
|
|
*
|
|
* Compute D1 and Z
|
|
*
|
|
CALL DCOPY( N, SD, 1, D5, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
|
|
*
|
|
IF( SRANGE ) THEN
|
|
NTEST = 29
|
|
IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
|
|
IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
|
|
IF( IU.LT.IL ) THEN
|
|
ITEMP = IU
|
|
IU = IL
|
|
IL = ITEMP
|
|
END IF
|
|
CALL DSTEMR( 'V', 'I', N, D5, WORK, VL, VU, IL, IU,
|
|
$ M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC,
|
|
$ WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
|
|
$ LIWORK-2*N, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,I)', IINFO,
|
|
$ N, JTYPE, IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 29 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do Tests 29 and 30
|
|
*
|
|
CALL DSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
|
|
$ M, RESULT( 29 ) )
|
|
*
|
|
* Call DSTEMR to compute D2, do tests.
|
|
*
|
|
* Compute D2
|
|
*
|
|
CALL DCOPY( N, SD, 1, D5, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
*
|
|
NTEST = 31
|
|
CALL DSTEMR( 'N', 'I', N, D5, WORK, VL, VU, IL, IU,
|
|
$ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC,
|
|
$ WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
|
|
$ LIWORK-2*N, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEMR(N,I)', IINFO,
|
|
$ N, JTYPE, IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 31 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do Test 31
|
|
*
|
|
TEMP1 = ZERO
|
|
TEMP2 = ZERO
|
|
*
|
|
DO 240 J = 1, IU - IL + 1
|
|
TEMP1 = MAX( TEMP1, ABS( D1( J ) ),
|
|
$ ABS( D2( J ) ) )
|
|
TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
|
|
240 CONTINUE
|
|
*
|
|
RESULT( 31 ) = TEMP2 / MAX( UNFL,
|
|
$ ULP*MAX( TEMP1, TEMP2 ) )
|
|
*
|
|
*
|
|
* Call DSTEMR(V,V) to compute D1 and Z, do tests.
|
|
*
|
|
* Compute D1 and Z
|
|
*
|
|
CALL DCOPY( N, SD, 1, D5, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
|
|
*
|
|
NTEST = 32
|
|
*
|
|
IF( N.GT.0 ) THEN
|
|
IF( IL.NE.1 ) THEN
|
|
VL = D2( IL ) - MAX( HALF*
|
|
$ ( D2( IL )-D2( IL-1 ) ), ULP*ANORM,
|
|
$ TWO*RTUNFL )
|
|
ELSE
|
|
VL = D2( 1 ) - MAX( HALF*( D2( N )-D2( 1 ) ),
|
|
$ ULP*ANORM, TWO*RTUNFL )
|
|
END IF
|
|
IF( IU.NE.N ) THEN
|
|
VU = D2( IU ) + MAX( HALF*
|
|
$ ( D2( IU+1 )-D2( IU ) ), ULP*ANORM,
|
|
$ TWO*RTUNFL )
|
|
ELSE
|
|
VU = D2( N ) + MAX( HALF*( D2( N )-D2( 1 ) ),
|
|
$ ULP*ANORM, TWO*RTUNFL )
|
|
END IF
|
|
ELSE
|
|
VL = ZERO
|
|
VU = ONE
|
|
END IF
|
|
*
|
|
CALL DSTEMR( 'V', 'V', N, D5, WORK, VL, VU, IL, IU,
|
|
$ M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC,
|
|
$ WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
|
|
$ LIWORK-2*N, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,V)', IINFO,
|
|
$ N, JTYPE, IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 32 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do Tests 32 and 33
|
|
*
|
|
CALL DSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
|
|
$ M, RESULT( 32 ) )
|
|
*
|
|
* Call DSTEMR to compute D2, do tests.
|
|
*
|
|
* Compute D2
|
|
*
|
|
CALL DCOPY( N, SD, 1, D5, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
*
|
|
NTEST = 34
|
|
CALL DSTEMR( 'N', 'V', N, D5, WORK, VL, VU, IL, IU,
|
|
$ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC,
|
|
$ WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
|
|
$ LIWORK-2*N, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEMR(N,V)', IINFO,
|
|
$ N, JTYPE, IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 34 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do Test 34
|
|
*
|
|
TEMP1 = ZERO
|
|
TEMP2 = ZERO
|
|
*
|
|
DO 250 J = 1, IU - IL + 1
|
|
TEMP1 = MAX( TEMP1, ABS( D1( J ) ),
|
|
$ ABS( D2( J ) ) )
|
|
TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
|
|
250 CONTINUE
|
|
*
|
|
RESULT( 34 ) = TEMP2 / MAX( UNFL,
|
|
$ ULP*MAX( TEMP1, TEMP2 ) )
|
|
ELSE
|
|
RESULT( 29 ) = ZERO
|
|
RESULT( 30 ) = ZERO
|
|
RESULT( 31 ) = ZERO
|
|
RESULT( 32 ) = ZERO
|
|
RESULT( 33 ) = ZERO
|
|
RESULT( 34 ) = ZERO
|
|
END IF
|
|
*
|
|
*
|
|
* Call DSTEMR(V,A) to compute D1 and Z, do tests.
|
|
*
|
|
* Compute D1 and Z
|
|
*
|
|
CALL DCOPY( N, SD, 1, D5, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
*
|
|
NTEST = 35
|
|
*
|
|
CALL DSTEMR( 'V', 'A', N, D5, WORK, VL, VU, IL, IU,
|
|
$ M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC,
|
|
$ WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
|
|
$ LIWORK-2*N, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,A)', IINFO, N,
|
|
$ JTYPE, IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 35 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do Tests 35 and 36
|
|
*
|
|
CALL DSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, M,
|
|
$ RESULT( 35 ) )
|
|
*
|
|
* Call DSTEMR to compute D2, do tests.
|
|
*
|
|
* Compute D2
|
|
*
|
|
CALL DCOPY( N, SD, 1, D5, 1 )
|
|
IF( N.GT.0 )
|
|
$ CALL DCOPY( N-1, SE, 1, WORK, 1 )
|
|
*
|
|
NTEST = 37
|
|
CALL DSTEMR( 'N', 'A', N, D5, WORK, VL, VU, IL, IU,
|
|
$ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC,
|
|
$ WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
|
|
$ LIWORK-2*N, IINFO )
|
|
IF( IINFO.NE.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9999 )'DSTEMR(N,A)', IINFO, N,
|
|
$ JTYPE, IOLDSD
|
|
INFO = ABS( IINFO )
|
|
IF( IINFO.LT.0 ) THEN
|
|
RETURN
|
|
ELSE
|
|
RESULT( 37 ) = ULPINV
|
|
GO TO 280
|
|
END IF
|
|
END IF
|
|
*
|
|
* Do Test 34
|
|
*
|
|
TEMP1 = ZERO
|
|
TEMP2 = ZERO
|
|
*
|
|
DO 260 J = 1, N
|
|
TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
|
|
TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
|
|
260 CONTINUE
|
|
*
|
|
RESULT( 37 ) = TEMP2 / MAX( UNFL,
|
|
$ ULP*MAX( TEMP1, TEMP2 ) )
|
|
END IF
|
|
270 CONTINUE
|
|
280 CONTINUE
|
|
NTESTT = NTESTT + NTEST
|
|
*
|
|
* End of Loop -- Check for RESULT(j) > THRESH
|
|
*
|
|
*
|
|
* Print out tests which fail.
|
|
*
|
|
DO 290 JR = 1, NTEST
|
|
IF( RESULT( JR ).GE.THRESH ) THEN
|
|
*
|
|
* If this is the first test to fail,
|
|
* print a header to the data file.
|
|
*
|
|
IF( NERRS.EQ.0 ) THEN
|
|
WRITE( NOUNIT, FMT = 9998 )'DST'
|
|
WRITE( NOUNIT, FMT = 9997 )
|
|
WRITE( NOUNIT, FMT = 9996 )
|
|
WRITE( NOUNIT, FMT = 9995 )'Symmetric'
|
|
WRITE( NOUNIT, FMT = 9994 )
|
|
*
|
|
* Tests performed
|
|
*
|
|
WRITE( NOUNIT, FMT = 9988 )
|
|
END IF
|
|
NERRS = NERRS + 1
|
|
WRITE( NOUNIT, FMT = 9990 )N, IOLDSD, JTYPE, JR,
|
|
$ RESULT( JR )
|
|
END IF
|
|
290 CONTINUE
|
|
300 CONTINUE
|
|
310 CONTINUE
|
|
*
|
|
* Summary
|
|
*
|
|
CALL DLASUM( 'DST', NOUNIT, NERRS, NTESTT )
|
|
RETURN
|
|
*
|
|
9999 FORMAT( ' DCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
|
|
$ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
|
|
*
|
|
9998 FORMAT( / 1X, A3, ' -- Real Symmetric eigenvalue problem' )
|
|
9997 FORMAT( ' Matrix types (see DCHKST for details): ' )
|
|
*
|
|
9996 FORMAT( / ' Special Matrices:',
|
|
$ / ' 1=Zero matrix. ',
|
|
$ ' 5=Diagonal: clustered entries.',
|
|
$ / ' 2=Identity matrix. ',
|
|
$ ' 6=Diagonal: large, evenly spaced.',
|
|
$ / ' 3=Diagonal: evenly spaced entries. ',
|
|
$ ' 7=Diagonal: small, evenly spaced.',
|
|
$ / ' 4=Diagonal: geometr. spaced entries.' )
|
|
9995 FORMAT( ' Dense ', A, ' Matrices:',
|
|
$ / ' 8=Evenly spaced eigenvals. ',
|
|
$ ' 12=Small, evenly spaced eigenvals.',
|
|
$ / ' 9=Geometrically spaced eigenvals. ',
|
|
$ ' 13=Matrix with random O(1) entries.',
|
|
$ / ' 10=Clustered eigenvalues. ',
|
|
$ ' 14=Matrix with large random entries.',
|
|
$ / ' 11=Large, evenly spaced eigenvals. ',
|
|
$ ' 15=Matrix with small random entries.' )
|
|
9994 FORMAT( ' 16=Positive definite, evenly spaced eigenvalues',
|
|
$ / ' 17=Positive definite, geometrically spaced eigenvlaues',
|
|
$ / ' 18=Positive definite, clustered eigenvalues',
|
|
$ / ' 19=Positive definite, small evenly spaced eigenvalues',
|
|
$ / ' 20=Positive definite, large evenly spaced eigenvalues',
|
|
$ / ' 21=Diagonally dominant tridiagonal, geometrically',
|
|
$ ' spaced eigenvalues' )
|
|
*
|
|
9990 FORMAT( ' N=', I5, ', seed=', 4( I4, ',' ), ' type ', I2,
|
|
$ ', test(', I2, ')=', G10.3 )
|
|
*
|
|
9988 FORMAT( / 'Test performed: see DCHKST for details.', / )
|
|
* End of DCHKST
|
|
*
|
|
END
|