diff options
Diffstat (limited to 'TESTING/EIG/zchkst2stg.f')
-rw-r--r-- | TESTING/EIG/zchkst2stg.f | 2145 |
1 files changed, 2145 insertions, 0 deletions
diff --git a/TESTING/EIG/zchkst2stg.f b/TESTING/EIG/zchkst2stg.f new file mode 100644 index 00000000..a1aaffbc --- /dev/null +++ b/TESTING/EIG/zchkst2stg.f @@ -0,0 +1,2145 @@ +*> \brief \b ZCHKST2STG +* +* @precisions fortran z -> c +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE ZCHKST2STG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, +* NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5, +* WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK, +* LWORK, RWORK, LRWORK, IWORK, LIWORK, RESULT, +* INFO ) +* +* .. Scalar Arguments .. +* INTEGER INFO, LDA, LDU, LIWORK, LRWORK, LWORK, NOUNIT, +* $ NSIZES, NTYPES +* DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER ISEED( 4 ), IWORK( * ), NN( * ) +* DOUBLE PRECISION D1( * ), D2( * ), D3( * ), D4( * ), D5( * ), +* $ RESULT( * ), RWORK( * ), SD( * ), SE( * ), +* $ WA1( * ), WA2( * ), WA3( * ), WR( * ) +* COMPLEX*16 A( LDA, * ), AP( * ), TAU( * ), U( LDU, * ), +* $ V( LDU, * ), VP( * ), WORK( * ), Z( LDU, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZCHKST2STG checks the Hermitian eigenvalue problem routines +*> using the 2-stage reduction techniques. Since the generation +*> of Q or the vectors is not available in this release, we only +*> compare the eigenvalue resulting when using the 2-stage to the +*> one considered as reference using the standard 1-stage reduction +*> ZHETRD. For that, we call the standard ZHETRD and compute D1 using +*> DSTEQR, then we call the 2-stage ZHETRD_2STAGE with Upper and Lower +*> and we compute D2 and D3 using DSTEQR and then we replaced tests +*> 3 and 4 by tests 11 and 12. test 1 and 2 remain to verify that +*> the 1-stage results are OK and can be trusted. +*> This testing routine will converge to the ZCHKST in the next +*> release when vectors and generation of Q will be implemented. +*> +*> ZHETRD factors A as U S U* , where * means conjugate transpose, +*> S is real symmetric tridiagonal, and U is unitary. +*> ZHETRD can use either just the lower or just the upper triangle +*> of A; ZCHKST2STG checks both cases. +*> U is represented as a product of Householder +*> transformations, whose vectors are stored in the first +*> n-1 columns of V, and whose scale factors are in TAU. +*> +*> ZHPTRD does the same as ZHETRD, except that A and V are stored +*> in "packed" format. +*> +*> ZUNGTR constructs the matrix U from the contents of V and TAU. +*> +*> ZUPGTR constructs the matrix U from the contents of VP and TAU. +*> +*> ZSTEQR factors S as Z D1 Z* , where Z is the unitary +*> matrix of eigenvectors and D1 is a diagonal matrix with +*> the eigenvalues on the diagonal. D2 is the matrix of +*> eigenvalues computed when Z is not computed. +*> +*> DSTERF computes D3, the matrix of eigenvalues, by the +*> PWK method, which does not yield eigenvectors. +*> +*> ZPTEQR factors S as Z4 D4 Z4* , for a +*> Hermitian positive definite tridiagonal matrix. +*> D5 is the matrix of eigenvalues computed when Z is not +*> computed. +*> +*> DSTEBZ computes selected eigenvalues. WA1, WA2, and +*> WA3 will denote eigenvalues computed to high +*> absolute accuracy, with different range options. +*> WR will denote eigenvalues computed to high relative +*> accuracy. +*> +*> ZSTEIN computes Y, the eigenvectors of S, given the +*> eigenvalues. +*> +*> ZSTEDC factors S as Z D1 Z* , where Z is the unitary +*> matrix of eigenvectors and D1 is a diagonal matrix with +*> the eigenvalues on the diagonal ('I' option). It may also +*> update an input unitary matrix, usually the output +*> from ZHETRD/ZUNGTR or ZHPTRD/ZUPGTR ('V' option). It may +*> also just compute eigenvalues ('N' option). +*> +*> ZSTEMR factors S as Z D1 Z* , where Z is the unitary +*> matrix of eigenvectors and D1 is a diagonal matrix with +*> the eigenvalues on the diagonal ('I' option). ZSTEMR +*> uses the Relatively Robust Representation whenever possible. +*> +*> When ZCHKST2STG is called, a number of matrix "sizes" ("n's") and a +*> number of matrix "types" are specified. For each size ("n") +*> and each type of matrix, one matrix will be generated and used +*> to test the Hermitian eigenroutines. For each matrix, a number +*> of tests will be performed: +*> +*> (1) | A - V S V* | / ( |A| n ulp ) ZHETRD( UPLO='U', ... ) +*> +*> (2) | I - UV* | / ( n ulp ) ZUNGTR( UPLO='U', ... ) +*> +*> (3) | A - V S V* | / ( |A| n ulp ) ZHETRD( UPLO='L', ... ) +*> replaced by | D1 - D2 | / ( |D1| ulp ) where D1 is the +*> eigenvalue matrix computed using S and D2 is the +*> eigenvalue matrix computed using S_2stage the output of +*> ZHETRD_2STAGE("N", "U",....). D1 and D2 are computed +*> via DSTEQR('N',...) +*> +*> (4) | I - UV* | / ( n ulp ) ZUNGTR( UPLO='L', ... ) +*> replaced by | D1 - D3 | / ( |D1| ulp ) where D1 is the +*> eigenvalue matrix computed using S and D3 is the +*> eigenvalue matrix computed using S_2stage the output of +*> ZHETRD_2STAGE("N", "L",....). D1 and D3 are computed +*> via DSTEQR('N',...) +*> +*> (5-8) Same as 1-4, but for ZHPTRD and ZUPGTR. +*> +*> (9) | S - Z D Z* | / ( |S| n ulp ) ZSTEQR('V',...) +*> +*> (10) | I - ZZ* | / ( n ulp ) ZSTEQR('V',...) +*> +*> (11) | D1 - D2 | / ( |D1| ulp ) ZSTEQR('N',...) +*> +*> (12) | D1 - D3 | / ( |D1| ulp ) DSTERF +*> +*> (13) 0 if the true eigenvalues (computed by sturm count) +*> of S are within THRESH of +*> those in D1. 2*THRESH if they are not. (Tested using +*> DSTECH) +*> +*> For S positive definite, +*> +*> (14) | S - Z4 D4 Z4* | / ( |S| n ulp ) ZPTEQR('V',...) +*> +*> (15) | I - Z4 Z4* | / ( n ulp ) ZPTEQR('V',...) +*> +*> (16) | D4 - D5 | / ( 100 |D4| ulp ) ZPTEQR('N',...) +*> +*> When S is also diagonally dominant by the factor gamma < 1, +*> +*> (17) max | D4(i) - WR(i) | / ( |D4(i)| omega ) , +*> i +*> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4 +*> DSTEBZ( 'A', 'E', ...) +*> +*> (18) | WA1 - D3 | / ( |D3| ulp ) DSTEBZ( 'A', 'E', ...) +*> +*> (19) ( max { min | WA2(i)-WA3(j) | } + +*> i j +*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) +*> i j +*> DSTEBZ( 'I', 'E', ...) +*> +*> (20) | S - Y WA1 Y* | / ( |S| n ulp ) DSTEBZ, ZSTEIN +*> +*> (21) | I - Y Y* | / ( n ulp ) DSTEBZ, ZSTEIN +*> +*> (22) | S - Z D Z* | / ( |S| n ulp ) ZSTEDC('I') +*> +*> (23) | I - ZZ* | / ( n ulp ) ZSTEDC('I') +*> +*> (24) | S - Z D Z* | / ( |S| n ulp ) ZSTEDC('V') +*> +*> (25) | I - ZZ* | / ( n ulp ) ZSTEDC('V') +*> +*> (26) | D1 - D2 | / ( |D1| ulp ) ZSTEDC('V') and +*> ZSTEDC('N') +*> +*> Test 27 is disabled at the moment because ZSTEMR does not +*> guarantee high relatvie accuracy. +*> +*> (27) max | D6(i) - WR(i) | / ( |D6(i)| omega ) , +*> i +*> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4 +*> ZSTEMR('V', 'A') +*> +*> (28) max | D6(i) - WR(i) | / ( |D6(i)| omega ) , +*> i +*> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4 +*> ZSTEMR('V', 'I') +*> +*> Tests 29 through 34 are disable at present because ZSTEMR +*> does not handle partial specturm requests. +*> +*> (29) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'I') +*> +*> (30) | I - ZZ* | / ( n ulp ) ZSTEMR('V', 'I') +*> +*> (31) ( max { min | WA2(i)-WA3(j) | } + +*> i j +*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) +*> i j +*> ZSTEMR('N', 'I') vs. CSTEMR('V', 'I') +*> +*> (32) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'V') +*> +*> (33) | I - ZZ* | / ( n ulp ) ZSTEMR('V', 'V') +*> +*> (34) ( max { min | WA2(i)-WA3(j) | } + +*> i j +*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) +*> i j +*> ZSTEMR('N', 'V') vs. CSTEMR('V', 'V') +*> +*> (35) | S - Z D Z* | / ( |S| n ulp ) ZSTEMR('V', 'A') +*> +*> (36) | I - ZZ* | / ( n ulp ) ZSTEMR('V', 'A') +*> +*> (37) ( max { min | WA2(i)-WA3(j) | } + +*> i j +*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) +*> i j +*> ZSTEMR('N', 'A') vs. CSTEMR('V', 'A') +*> +*> The "sizes" are specified by an array NN(1:NSIZES); the value of +*> each element NN(j) specifies one size. +*> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); +*> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. +*> Currently, the list of possible types is: +*> +*> (1) The zero matrix. +*> (2) The identity matrix. +*> +*> (3) A diagonal matrix with evenly spaced entries +*> 1, ..., ULP and random signs. +*> (ULP = (first number larger than 1) - 1 ) +*> (4) A diagonal matrix with geometrically spaced entries +*> 1, ..., ULP and random signs. +*> (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP +*> and random signs. +*> +*> (6) Same as (4), but multiplied by SQRT( overflow threshold ) +*> (7) Same as (4), but multiplied by SQRT( underflow threshold ) +*> +*> (8) A matrix of the form U* D U, where U is unitary and +*> D has evenly spaced entries 1, ..., ULP with random signs +*> on the diagonal. +*> +*> (9) A matrix of the form U* D U, where U is unitary and +*> D has geometrically spaced entries 1, ..., ULP with random +*> signs on the diagonal. +*> +*> (10) A matrix of the form U* D U, where U is unitary and +*> D has "clustered" entries 1, ULP,..., ULP with random +*> signs on the diagonal. +*> +*> (11) Same as (8), but multiplied by SQRT( overflow threshold ) +*> (12) Same as (8), but multiplied by SQRT( underflow threshold ) +*> +*> (13) Hermitian matrix with random entries chosen from (-1,1). +*> (14) Same as (13), but multiplied by SQRT( overflow threshold ) +*> (15) Same as (13), but multiplied by SQRT( underflow threshold ) +*> (16) Same as (8), but diagonal elements are all positive. +*> (17) Same as (9), but diagonal elements are all positive. +*> (18) Same as (10), but diagonal elements are all positive. +*> (19) Same as (16), but multiplied by SQRT( overflow threshold ) +*> (20) Same as (16), but multiplied by SQRT( underflow threshold ) +*> (21) A diagonally dominant tridiagonal matrix with geometrically +*> spaced diagonal entries 1, ..., ULP. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] NSIZES +*> \verbatim +*> NSIZES is INTEGER +*> The number of sizes of matrices to use. If it is zero, +*> ZCHKST2STG 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, ZCHKST2STG +*> 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 ZCHKST2STG 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 COMPLEX*16 array of +*> dimension ( LDA , max(NN) ) +*> Used to hold the matrix whose eigenvalues are to be +*> computed. On exit, A contains the last matrix actually +*> used. +*> \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 COMPLEX*16 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 ZHETRD. +*> 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 +*> ZHETRD. 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 ZSTEQR 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 ZSTEQR 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 ZPTEQR(V). +*> ZPTEQR 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 ZPTEQR(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 COMPLEX*16 array of +*> dimension( LDU, max(NN) ). +*> The unitary matrix computed by ZHETRD + ZUNGTR. +*> \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 COMPLEX*16 array of +*> dimension( LDU, max(NN) ). +*> The Housholder vectors computed by ZHETRD in reducing A to +*> tridiagonal form. The vectors computed with UPLO='U' are +*> in the upper triangle, and the vectors computed with UPLO='L' +*> are in the lower triangle. (As described in ZHETRD, the +*> sub- and superdiagonal are not set to 1, although the +*> true Householder vector has a 1 in that position. The +*> routines that use V, such as ZUNGTR, set those entries to +*> 1 before using them, and then restore them later.) +*> \endverbatim +*> +*> \param[out] VP +*> \verbatim +*> VP is COMPLEX*16 array of +*> dimension( max(NN)*max(NN+1)/2 ) +*> The matrix V stored in packed format. +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is COMPLEX*16 array of +*> dimension( max(NN) ) +*> The Householder factors computed by ZHETRD in reducing A +*> to tridiagonal form. +*> \endverbatim +*> +*> \param[out] Z +*> \verbatim +*> Z is COMPLEX*16 array of +*> dimension( LDU, max(NN) ). +*> The unitary matrix of eigenvectors computed by ZSTEQR, +*> ZPTEQR, and ZSTEIN. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 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] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array +*> \endverbatim +*> +*> \param[in] LRWORK +*> \verbatim +*> LRWORK is INTEGER +*> The number of entries in LRWORK (dimension( ??? ) +*> \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 ZLATMR, CLATMS, ZHETRD, ZUNGTR, ZSTEQR, DSTERF, +*> or ZUNMC2 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 November 2011 +* +*> \ingroup complex16_eig +* +* ===================================================================== + SUBROUTINE ZCHKST2STG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, + $ NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5, + $ WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK, + $ LWORK, RWORK, LRWORK, IWORK, LIWORK, RESULT, + $ INFO ) +* +* -- LAPACK test routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, LDU, LIWORK, LRWORK, LWORK, NOUNIT, + $ NSIZES, NTYPES + DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER ISEED( 4 ), IWORK( * ), NN( * ) + DOUBLE PRECISION D1( * ), D2( * ), D3( * ), D4( * ), D5( * ), + $ RESULT( * ), RWORK( * ), SD( * ), SE( * ), + $ WA1( * ), WA2( * ), WA3( * ), WR( * ) + COMPLEX*16 A( LDA, * ), AP( * ), TAU( * ), U( LDU, * ), + $ V( LDU, * ), VP( * ), WORK( * ), Z( LDU, * ) +* .. +* +* ===================================================================== +* +* .. 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 ) + COMPLEX*16 CZERO, CONE + PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), + $ CONE = ( 1.0D+0, 0.0D+0 ) ) + DOUBLE PRECISION HALF + PARAMETER ( HALF = ONE / TWO ) + INTEGER MAXTYP + PARAMETER ( MAXTYP = 21 ) + LOGICAL CRANGE + PARAMETER ( CRANGE = .FALSE. ) + LOGICAL CREL + PARAMETER ( CREL = .FALSE. ) +* .. +* .. Local Scalars .. + LOGICAL BADNN, TRYRAC + INTEGER I, IINFO, IL, IMODE, INDE, INDRWK, ITEMP, + $ ITYPE, IU, J, JC, JR, JSIZE, JTYPE, LGN, + $ LIWEDC, LOG2UI, LRWEDC, LWEDC, M, M2, M3, + $ MTYPES, N, NAP, NBLOCK, NERRS, NMATS, NMAX, + $ NSPLIT, NTEST, NTESTT, LH, LW + 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, DLASUM, DSTEBZ, DSTECH, DSTERF, + $ XERBLA, ZCOPY, ZHET21, ZHETRD, ZHPT21, ZHPTRD, + $ ZLACPY, ZLASET, ZLATMR, ZLATMS, ZPTEQR, ZSTEDC, + $ ZSTEMR, ZSTEIN, ZSTEQR, ZSTT21, ZSTT22, ZUNGTR, + $ ZUPGTR, ZHETRD_2STAGE +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, DCONJG, 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, 'ZHETRD', '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( 'ZCHKST2STG', -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 + LRWEDC = 1 + 3*N + 2*N*LGN + 4*N**2 + LIWEDC = 6 + 6*N + 5*N*LGN + ELSE + LWEDC = 8 + LRWEDC = 7 + 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 Hermitian, w/ eigenvalues +* =6 random (none) +* =7 random diagonal +* =8 random Hermitian +* =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 ZLASET( 'Full', LDA, N, CZERO, CZERO, 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 ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND, + $ ANORM, 0, 0, 'N', A, LDA, WORK, IINFO ) +* +* + ELSE IF( ITYPE.EQ.5 ) THEN +* +* Hermitian, eigenvalues specified +* + CALL ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND, + $ ANORM, N, N, 'N', A, LDA, WORK, IINFO ) +* + ELSE IF( ITYPE.EQ.7 ) THEN +* +* Diagonal, random eigenvalues +* + CALL ZLATMR( N, N, 'S', ISEED, 'H', WORK, 6, ONE, CONE, + $ 'T', 'N', WORK( N+1 ), 1, ONE, + $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0, + $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) +* + ELSE IF( ITYPE.EQ.8 ) THEN +* +* Hermitian, random eigenvalues +* + CALL ZLATMR( N, N, 'S', ISEED, 'H', WORK, 6, ONE, CONE, + $ 'T', 'N', WORK( N+1 ), 1, ONE, + $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N, + $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO ) +* + ELSE IF( ITYPE.EQ.9 ) THEN +* +* Positive definite, eigenvalues specified. +* + CALL ZLATMS( N, N, 'S', ISEED, 'P', RWORK, IMODE, COND, + $ ANORM, N, N, 'N', A, LDA, WORK, IINFO ) +* + ELSE IF( ITYPE.EQ.10 ) THEN +* +* Positive definite tridiagonal, eigenvalues specified. +* + CALL ZLATMS( N, N, 'S', ISEED, 'P', RWORK, IMODE, COND, + $ ANORM, 1, 1, 'N', A, LDA, WORK, IINFO ) + DO 90 I = 2, N + TEMP1 = ABS( A( I-1, I ) ) + TEMP2 = SQRT( ABS( A( I-1, I-1 )*A( I, I ) ) ) + IF( TEMP1.GT.HALF*TEMP2 ) THEN + A( I-1, I ) = A( I-1, I )* + $ ( HALF*TEMP2 / ( UNFL+TEMP1 ) ) + A( I, I-1 ) = DCONJG( 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 ZHETRD and ZUNGTR to compute S and U from +* upper triangle. +* + CALL ZLACPY( 'U', N, N, A, LDA, V, LDU ) +* + NTEST = 1 + CALL ZHETRD( 'U', N, V, LDU, SD, SE, TAU, WORK, LWORK, + $ IINFO ) +* + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZHETRD(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 ZLACPY( 'U', N, N, V, LDU, U, LDU ) +* + NTEST = 2 + CALL ZUNGTR( 'U', N, U, LDU, TAU, WORK, LWORK, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZUNGTR(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 ZHET21( 2, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V, + $ LDU, TAU, WORK, RWORK, RESULT( 1 ) ) + CALL ZHET21( 3, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V, + $ LDU, TAU, WORK, RWORK, RESULT( 2 ) ) +* +* Compute D1 the eigenvalues resulting from the tridiagonal +* form using the standard 1-stage algorithm and use it as a +* reference to compare with the 2-stage technique +* +* Compute D1 from the 1-stage and used as reference for the +* 2-stage +* + CALL DCOPY( N, SD, 1, D1, 1 ) + IF( N.GT.0 ) + $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) +* + CALL ZSTEQR( 'N', N, D1, RWORK, WORK, LDU, RWORK( N+1 ), + $ IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEQR(N)', IINFO, N, JTYPE, + $ IOLDSD + INFO = ABS( IINFO ) + IF( IINFO.LT.0 ) THEN + RETURN + ELSE + RESULT( 3 ) = ULPINV + GO TO 280 + END IF + END IF +* +* 2-STAGE TRD Upper case is used to compute D2. +* Note to set SD and SE to zero to be sure not reusing +* the one from above. Compare it with D1 computed +* using the 1-stage. +* + CALL DLASET( 'Full', N, 1, ZERO, ZERO, SD, 1 ) + CALL DLASET( 'Full', N, 1, ZERO, ZERO, SE, 1 ) + CALL ZLACPY( 'U', N, N, A, LDA, V, LDU ) + LH = MAX(1, 4*N) + LW = LWORK - LH + CALL ZHETRD_2STAGE( 'N', "U", N, V, LDU, SD, SE, TAU, + $ WORK, LH, WORK( LH+1 ), LW, IINFO ) +* +* Compute D2 from the 2-stage Upper case +* + CALL DCOPY( N, SD, 1, D2, 1 ) + IF( N.GT.0 ) + $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) +* + NTEST = 3 + CALL ZSTEQR( 'N', N, D2, RWORK, WORK, LDU, RWORK( N+1 ), + $ IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEQR(N)', IINFO, N, JTYPE, + $ IOLDSD + INFO = ABS( IINFO ) + IF( IINFO.LT.0 ) THEN + RETURN + ELSE + RESULT( 3 ) = ULPINV + GO TO 280 + END IF + END IF +* +* 2-STAGE TRD Lower case is used to compute D3. +* Note to set SD and SE to zero to be sure not reusing +* the one from above. Compare it with D1 computed +* using the 1-stage. +* + CALL DLASET( 'Full', N, 1, ZERO, ZERO, SD, 1 ) + CALL DLASET( 'Full', N, 1, ZERO, ZERO, SE, 1 ) + CALL ZLACPY( 'L', N, N, A, LDA, V, LDU ) + CALL ZHETRD_2STAGE( 'N', "L", N, V, LDU, SD, SE, TAU, + $ WORK, LH, WORK( LH+1 ), LW, IINFO ) +* +* Compute D3 from the 2-stage Upper case +* + CALL DCOPY( N, SD, 1, D3, 1 ) + IF( N.GT.0 ) + $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) +* + NTEST = 4 + CALL ZSTEQR( 'N', N, D3, RWORK, WORK, LDU, RWORK( N+1 ), + $ IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEQR(N)', IINFO, N, JTYPE, + $ IOLDSD + INFO = ABS( IINFO ) + IF( IINFO.LT.0 ) THEN + RETURN + ELSE + RESULT( 4 ) = ULPINV + GO TO 280 + END IF + END IF +* +* +* Do Tests 3 and 4 which are similar to 11 and 12 but with the +* D1 computed using the standard 1-stage reduction as reference +* + NTEST = 4 + TEMP1 = ZERO + TEMP2 = ZERO + TEMP3 = ZERO + TEMP4 = ZERO +* + DO 151 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 ) ) ) + 151 CONTINUE +* + RESULT( 3 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) + RESULT( 4 ) = TEMP4 / MAX( UNFL, ULP*MAX( TEMP3, TEMP4 ) ) +* +* Skip the DSYTRD for lower that since we replaced its testing +* 3 and 4 by the 2-stage one. + GOTO 101 +* +* Call ZHETRD and ZUNGTR to compute S and U from +* lower triangle, do tests. +* + CALL ZLACPY( 'L', N, N, A, LDA, V, LDU ) +* + NTEST = 3 + CALL ZHETRD( 'L', N, V, LDU, SD, SE, TAU, WORK, LWORK, + $ IINFO ) +* + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZHETRD(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 ZLACPY( 'L', N, N, V, LDU, U, LDU ) +* + NTEST = 4 + CALL ZUNGTR( 'L', N, U, LDU, TAU, WORK, LWORK, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZUNGTR(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 ZHET21( 2, 'Lower', N, 1, A, LDA, SD, SE, U, LDU, V, + $ LDU, TAU, WORK, RWORK, RESULT( 3 ) ) + CALL ZHET21( 3, 'Lower', N, 1, A, LDA, SD, SE, U, LDU, V, + $ LDU, TAU, WORK, RWORK, RESULT( 4 ) ) +* +*after skipping old tests 3 4 back to the normal +* + 101 CONTINUE +* +* 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 ZHPTRD and ZUPGTR to compute S and U from AP +* + CALL ZCOPY( NAP, AP, 1, VP, 1 ) +* + NTEST = 5 + CALL ZHPTRD( 'U', N, VP, SD, SE, TAU, IINFO ) +* + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZHPTRD(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 ZUPGTR( 'U', N, VP, TAU, U, LDU, WORK, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZUPGTR(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 ZHPT21( 2, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU, + $ WORK, RWORK, RESULT( 5 ) ) + CALL ZHPT21( 3, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU, + $ WORK, RWORK, 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 ZHPTRD and ZUPGTR to compute S and U from AP +* + CALL ZCOPY( NAP, AP, 1, VP, 1 ) +* + NTEST = 7 + CALL ZHPTRD( 'L', N, VP, SD, SE, TAU, IINFO ) +* + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZHPTRD(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 ZUPGTR( 'L', N, VP, TAU, U, LDU, WORK, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZUPGTR(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 ZHPT21( 2, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU, + $ WORK, RWORK, RESULT( 7 ) ) + CALL ZHPT21( 3, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU, + $ WORK, RWORK, RESULT( 8 ) ) +* +* Call ZSTEQR 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, RWORK, 1 ) + CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) +* + NTEST = 9 + CALL ZSTEQR( 'V', N, D1, RWORK, Z, LDU, RWORK( N+1 ), + $ IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEQR(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, RWORK, 1 ) +* + NTEST = 11 + CALL ZSTEQR( 'N', N, D2, RWORK, WORK, LDU, RWORK( N+1 ), + $ IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEQR(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, RWORK, 1 ) +* + NTEST = 12 + CALL DSTERF( N, D3, RWORK, 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 ZSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, RWORK, + $ 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, RWORK, 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 ZPTEQR +* 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, RWORK, 1 ) + CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) +* + NTEST = 14 + CALL ZPTEQR( 'V', N, D4, RWORK, Z, LDU, RWORK( N+1 ), + $ IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZPTEQR(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 ZSTT21( N, 0, SD, SE, D4, DUMMA, Z, LDU, WORK, + $ RWORK, RESULT( 14 ) ) +* +* Compute D5 +* + CALL DCOPY( N, SD, 1, D5, 1 ) + IF( N.GT.0 ) + $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) +* + NTEST = 16 + CALL ZPTEQR( 'N', N, D5, RWORK, Z, LDU, RWORK( N+1 ), + $ IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZPTEQR(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 ), + $ RWORK, 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 ), RWORK, + $ 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 ), + $ RWORK, 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 ), + $ RWORK, 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 ZSTEIN 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 ), RWORK, + $ 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 ZSTEIN( N, SD, SE, M, WA1, IWORK( 1 ), IWORK( N+1 ), Z, + $ LDU, RWORK, IWORK( 2*N+1 ), IWORK( 3*N+1 ), + $ IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEIN', 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 ZSTT21( N, 0, SD, SE, WA1, DUMMA, Z, LDU, WORK, RWORK, + $ RESULT( 20 ) ) +* +* Call ZSTEDC(I) to compute D1 and Z, do tests. +* +* Compute D1 and Z +* + INDE = 1 + INDRWK = INDE + N + CALL DCOPY( N, SD, 1, D1, 1 ) + IF( N.GT.0 ) + $ CALL DCOPY( N-1, SE, 1, RWORK( INDE ), 1 ) + CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) +* + NTEST = 22 + CALL ZSTEDC( 'I', N, D1, RWORK( INDE ), Z, LDU, WORK, LWEDC, + $ RWORK( INDRWK ), LRWEDC, IWORK, LIWEDC, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEDC(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 ZSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, RWORK, + $ RESULT( 22 ) ) +* +* Call ZSTEDC(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, RWORK( INDE ), 1 ) + CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) +* + NTEST = 24 + CALL ZSTEDC( 'V', N, D1, RWORK( INDE ), Z, LDU, WORK, LWEDC, + $ RWORK( INDRWK ), LRWEDC, IWORK, LIWEDC, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEDC(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 ZSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, RWORK, + $ RESULT( 24 ) ) +* +* Call ZSTEDC(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, RWORK( INDE ), 1 ) + CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) +* + NTEST = 26 + CALL ZSTEDC( 'N', N, D2, RWORK( INDE ), Z, LDU, WORK, LWEDC, + $ RWORK( INDRWK ), LRWEDC, IWORK, LIWEDC, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEDC(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 ZSTEMR if IEEE compliant +* + IF( ILAENV( 10, 'ZSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'ZSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN +* +* Call ZSTEMR, 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. CREL ) THEN + NTEST = 27 + ABSTOL = UNFL + UNFL + CALL ZSTEMR( 'V', 'A', N, SD, SE, VL, VU, IL, IU, + $ M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC, + $ RWORK, LRWORK, IWORK( 2*N+1 ), LWORK-2*N, + $ IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(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( CRANGE ) THEN + NTEST = 28 + ABSTOL = UNFL + UNFL + CALL ZSTEMR( 'V', 'I', N, SD, SE, VL, VU, IL, IU, + $ M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC, + $ RWORK, LRWORK, IWORK( 2*N+1 ), + $ LWORK-2*N, IINFO ) +* + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(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 ZSTEMR(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, RWORK, 1 ) + CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) +* + IF( CRANGE ) 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 ZSTEMR( 'V', 'I', N, D5, RWORK, VL, VU, IL, IU, + $ M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC, + $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ), + $ LIWORK-2*N, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(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 ZSTEMR to compute D2, do tests. +* +* Compute D2 +* + CALL DCOPY( N, SD, 1, D5, 1 ) + IF( N.GT.0 ) + $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) +* + NTEST = 31 + CALL ZSTEMR( 'N', 'I', N, D5, RWORK, VL, VU, IL, IU, + $ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC, + $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ), + $ LIWORK-2*N, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(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 ZSTEMR(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, RWORK, 1 ) + CALL ZLASET( 'Full', N, N, CZERO, CONE, 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 ZSTEMR( 'V', 'V', N, D5, RWORK, VL, VU, IL, IU, + $ M, D1, Z, LDU, M, IWORK( 1 ), TRYRAC, + $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ), + $ LIWORK-2*N, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(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 ZSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, + $ M, RWORK, RESULT( 32 ) ) +* +* Call ZSTEMR to compute D2, do tests. +* +* Compute D2 +* + CALL DCOPY( N, SD, 1, D5, 1 ) + IF( N.GT.0 ) + $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) +* + NTEST = 34 + CALL ZSTEMR( 'N', 'V', N, D5, RWORK, VL, VU, IL, IU, + $ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC, + $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ), + $ LIWORK-2*N, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(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 ZSTEMR(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, RWORK, 1 ) +* + NTEST = 35 +* + CALL ZSTEMR( 'V', 'A', N, D5, RWORK, VL, VU, IL, IU, + $ M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC, + $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ), + $ LIWORK-2*N, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(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 ZSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, M, + $ RWORK, RESULT( 35 ) ) +* +* Call ZSTEMR to compute D2, do tests. +* +* Compute D2 +* + CALL DCOPY( N, SD, 1, D5, 1 ) + IF( N.GT.0 ) + $ CALL DCOPY( N-1, SE, 1, RWORK, 1 ) +* + NTEST = 37 + CALL ZSTEMR( 'N', 'A', N, D5, RWORK, VL, VU, IL, IU, + $ M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC, + $ RWORK( N+1 ), LRWORK-N, IWORK( 2*N+1 ), + $ LIWORK-2*N, IINFO ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'ZSTEMR(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 )'ZST' + WRITE( NOUNIT, FMT = 9997 ) + WRITE( NOUNIT, FMT = 9996 ) + WRITE( NOUNIT, FMT = 9995 )'Hermitian' + WRITE( NOUNIT, FMT = 9994 ) +* +* Tests performed +* + WRITE( NOUNIT, FMT = 9987 ) + END IF + NERRS = NERRS + 1 + IF( RESULT( JR ).LT.10000.0D0 ) THEN + WRITE( NOUNIT, FMT = 9989 )N, JTYPE, IOLDSD, JR, + $ RESULT( JR ) + ELSE + WRITE( NOUNIT, FMT = 9988 )N, JTYPE, IOLDSD, JR, + $ RESULT( JR ) + END IF + END IF + 290 CONTINUE + 300 CONTINUE + 310 CONTINUE +* +* Summary +* + CALL DLASUM( 'ZST', NOUNIT, NERRS, NTESTT ) + RETURN +* + 9999 FORMAT( ' ZCHKST2STG: ', A, ' returned INFO=', I6, '.', / 9X, + $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) +* + 9998 FORMAT( / 1X, A3, ' -- Complex Hermitian eigenvalue problem' ) + 9997 FORMAT( ' Matrix types (see ZCHKST2STG 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' ) +* + 9989 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', + $ 4( I4, ',' ), ' result ', I3, ' is', 0P, F8.2 ) + 9988 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=', + $ 4( I4, ',' ), ' result ', I3, ' is', 1P, D10.3 ) +* + 9987 FORMAT( / 'Test performed: see ZCHKST2STG for details.', / ) +* End of ZCHKST2STG +* + END |