*> \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