floraachy
/
OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
OpenBLAS/lapack-netlib/TESTING/EIG/cchkst2stg.f

2085 lines
71 KiB
Fortran
Raw Blame History

*> \brief \b CCHKST2STG
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE CCHKST2STG( 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
* REAL THRESH
* ..
* .. Array Arguments ..
* LOGICAL DOTYPE( * )
* INTEGER ISEED( 4 ), IWORK( * ), NN( * )
* REAL D1( * ), D2( * ), D3( * ), D4( * ), D5( * ),
* $ RESULT( * ), RWORK( * ), SD( * ), SE( * ),
* $ WA1( * ), WA2( * ), WA3( * ), WR( * )
* COMPLEX A( LDA, * ), AP( * ), TAU( * ), U( LDU, * ),
* $ V( LDU, * ), VP( * ), WORK( * ), Z( LDU, * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> CCHKST2STG 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
*> CHETRD. For that, we call the standard CHETRD and compute D1 using
*> DSTEQR, then we call the 2-stage CHETRD_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 CCHKST in the next
*> release when vectors and generation of Q will be implemented.
*>
*> CHETRD factors A as U S U* , where * means conjugate transpose,
*> S is real symmetric tridiagonal, and U is unitary.
*> CHETRD can use either just the lower or just the upper triangle
*> of A; CCHKST2STG 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.
*>
*> CHPTRD does the same as CHETRD, except that A and V are stored
*> in "packed" format.
*>
*> CUNGTR constructs the matrix U from the contents of V and TAU.
*>
*> CUPGTR constructs the matrix U from the contents of VP and TAU.
*>
*> CSTEQR 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.
*>
*> SSTERF computes D3, the matrix of eigenvalues, by the
*> PWK method, which does not yield eigenvectors.
*>
*> CPTEQR 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.
*>
*> SSTEBZ computes selected eigenvalues. WA1, WA2, and
*> WA3 will denote eigenvalues computed to high
*> absolute accuracy, with different range options.
*> WR will denote eigenvalues computed to high relative
*> accuracy.
*>
*> CSTEIN computes Y, the eigenvectors of S, given the
*> eigenvalues.
*>
*> CSTEDC 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 CHETRD/CUNGTR or CHPTRD/CUPGTR ('V' option). It may
*> also just compute eigenvalues ('N' option).
*>
*> CSTEMR 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). CSTEMR
*> uses the Relatively Robust Representation whenever possible.
*>
*> When CCHKST2STG 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 ) CHETRD( UPLO='U', ... )
*>
*> (2) | I - UV* | / ( n ulp ) CUNGTR( UPLO='U', ... )
*>
*> (3) | A - V S V* | / ( |A| n ulp ) CHETRD( 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
*> CHETRD_2STAGE("N", "U",....). D1 and D2 are computed
*> via DSTEQR('N',...)
*>
*> (4) | I - UV* | / ( n ulp ) CUNGTR( 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
*> CHETRD_2STAGE("N", "L",....). D1 and D3 are computed
*> via DSTEQR('N',...)
*>
*> (5-8) Same as 1-4, but for CHPTRD and CUPGTR.
*>
*> (9) | S - Z D Z* | / ( |S| n ulp ) CSTEQR('V',...)
*>
*> (10) | I - ZZ* | / ( n ulp ) CSTEQR('V',...)
*>
*> (11) | D1 - D2 | / ( |D1| ulp ) CSTEQR('N',...)
*>
*> (12) | D1 - D3 | / ( |D1| ulp ) SSTERF
*>
*> (13) 0 if the true eigenvalues (computed by sturm count)
*> of S are within THRESH of
*> those in D1. 2*THRESH if they are not. (Tested using
*> SSTECH)
*>
*> For S positive definite,
*>
*> (14) | S - Z4 D4 Z4* | / ( |S| n ulp ) CPTEQR('V',...)
*>
*> (15) | I - Z4 Z4* | / ( n ulp ) CPTEQR('V',...)
*>
*> (16) | D4 - D5 | / ( 100 |D4| ulp ) CPTEQR('N',...)
*>
*> When S is also diagonally dominant by the factor gamma < 1,
*>
*> (17) max | D4(i) - WR(i) | / ( |D4(i)| omega ) ,
*> i
*> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
*> SSTEBZ( 'A', 'E', ...)
*>
*> (18) | WA1 - D3 | / ( |D3| ulp ) SSTEBZ( 'A', 'E', ...)
*>
*> (19) ( max { min | WA2(i)-WA3(j) | } +
*> i j
*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
*> i j
*> SSTEBZ( 'I', 'E', ...)
*>
*> (20) | S - Y WA1 Y* | / ( |S| n ulp ) SSTEBZ, CSTEIN
*>
*> (21) | I - Y Y* | / ( n ulp ) SSTEBZ, CSTEIN
*>
*> (22) | S - Z D Z* | / ( |S| n ulp ) CSTEDC('I')
*>
*> (23) | I - ZZ* | / ( n ulp ) CSTEDC('I')
*>
*> (24) | S - Z D Z* | / ( |S| n ulp ) CSTEDC('V')
*>
*> (25) | I - ZZ* | / ( n ulp ) CSTEDC('V')
*>
*> (26) | D1 - D2 | / ( |D1| ulp ) CSTEDC('V') and
*> CSTEDC('N')
*>
*> Test 27 is disabled at the moment because CSTEMR 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
*> CSTEMR('V', 'A')
*>
*> (28) max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
*> i
*> omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
*> CSTEMR('V', 'I')
*>
*> Tests 29 through 34 are disable at present because CSTEMR
*> does not handle partial spectrum requests.
*>
*> (29) | S - Z D Z* | / ( |S| n ulp ) CSTEMR('V', 'I')
*>
*> (30) | I - ZZ* | / ( n ulp ) CSTEMR('V', 'I')
*>
*> (31) ( max { min | WA2(i)-WA3(j) | } +
*> i j
*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
*> i j
*> CSTEMR('N', 'I') vs. CSTEMR('V', 'I')
*>
*> (32) | S - Z D Z* | / ( |S| n ulp ) CSTEMR('V', 'V')
*>
*> (33) | I - ZZ* | / ( n ulp ) CSTEMR('V', 'V')
*>
*> (34) ( max { min | WA2(i)-WA3(j) | } +
*> i j
*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
*> i j
*> CSTEMR('N', 'V') vs. CSTEMR('V', 'V')
*>
*> (35) | S - Z D Z* | / ( |S| n ulp ) CSTEMR('V', 'A')
*>
*> (36) | I - ZZ* | / ( n ulp ) CSTEMR('V', 'A')
*>
*> (37) ( max { min | WA2(i)-WA3(j) | } +
*> i j
*> max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
*> i j
*> CSTEMR('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,
*> CCHKST2STG 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, CCHKST2STG
*> 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 CCHKST2STG to continue the same random number
*> sequence.
*> \endverbatim
*>
*> \param[in] THRESH
*> \verbatim
*> THRESH is REAL
*> A test will count as "failed" if the "error", computed as
*> described above, exceeds THRESH. Note that the error
*> is scaled to be O(1), so THRESH should be a reasonably
*> small multiple of 1, e.g., 10 or 100. In particular,
*> it should not depend on the precision (single vs. double)
*> or the size of the matrix. It must be at least zero.
*> \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 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 array of
*> dimension( max(NN)*max(NN+1)/2 )
*> The matrix A stored in packed format.
*> \endverbatim
*>
*> \param[out] SD
*> \verbatim
*> SD is REAL array of
*> dimension( max(NN) )
*> The diagonal of the tridiagonal matrix computed by CHETRD.
*> On exit, SD and SE contain the tridiagonal form of the
*> matrix in A.
*> \endverbatim
*>
*> \param[out] SE
*> \verbatim
*> SE is REAL array of
*> dimension( max(NN) )
*> The off-diagonal of the tridiagonal matrix computed by
*> CHETRD. On exit, SD and SE contain the tridiagonal form of
*> the matrix in A.
*> \endverbatim
*>
*> \param[out] D1
*> \verbatim
*> D1 is REAL array of
*> dimension( max(NN) )
*> The eigenvalues of A, as computed by CSTEQR simultaneously
*> with Z. On exit, the eigenvalues in D1 correspond with the
*> matrix in A.
*> \endverbatim
*>
*> \param[out] D2
*> \verbatim
*> D2 is REAL array of
*> dimension( max(NN) )
*> The eigenvalues of A, as computed by CSTEQR if Z is not
*> computed. On exit, the eigenvalues in D2 correspond with
*> the matrix in A.
*> \endverbatim
*>
*> \param[out] D3
*> \verbatim
*> D3 is REAL array of
*> dimension( max(NN) )
*> The eigenvalues of A, as computed by SSTERF. On exit, the
*> eigenvalues in D3 correspond with the matrix in A.
*> \endverbatim
*>
*> \param[out] D4
*> \verbatim
*> D4 is REAL array of
*> dimension( max(NN) )
*> The eigenvalues of A, as computed by CPTEQR(V).
*> CPTEQR 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 REAL array of
*> dimension( max(NN) )
*> The eigenvalues of A, as computed by CPTEQR(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 REAL array of
*> dimension( max(NN) )
*> All eigenvalues of A, computed to high
*> absolute accuracy, with different range options.
*> as computed by SSTEBZ.
*> \endverbatim
*>
*> \param[out] WA2
*> \verbatim
*> WA2 is REAL array of
*> dimension( max(NN) )
*> Selected eigenvalues of A, computed to high
*> absolute accuracy, with different range options.
*> as computed by SSTEBZ.
*> Choose random values for IL and IU, and ask for the
*> IL-th through IU-th eigenvalues.
*> \endverbatim
*>
*> \param[out] WA3
*> \verbatim
*> WA3 is REAL array of
*> dimension( max(NN) )
*> Selected eigenvalues of A, computed to high
*> absolute accuracy, with different range options.
*> as computed by SSTEBZ.
*> 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 REAL array of
*> dimension( max(NN) )
*> All eigenvalues of A, computed to high
*> absolute accuracy, with different options.
*> as computed by SSTEBZ.
*> \endverbatim
*>
*> \param[out] U
*> \verbatim
*> U is COMPLEX array of
*> dimension( LDU, max(NN) ).
*> The unitary matrix computed by CHETRD + CUNGTR.
*> \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 array of
*> dimension( LDU, max(NN) ).
*> The Housholder vectors computed by CHETRD 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 CHETRD, 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 CUNGTR, set those entries to
*> 1 before using them, and then restore them later.)
*> \endverbatim
*>
*> \param[out] VP
*> \verbatim
*> VP is COMPLEX array of
*> dimension( max(NN)*max(NN+1)/2 )
*> The matrix V stored in packed format.
*> \endverbatim
*>
*> \param[out] TAU
*> \verbatim
*> TAU is COMPLEX array of
*> dimension( max(NN) )
*> The Householder factors computed by CHETRD in reducing A
*> to tridiagonal form.
*> \endverbatim
*>
*> \param[out] Z
*> \verbatim
*> Z is COMPLEX array of
*> dimension( LDU, max(NN) ).
*> The unitary matrix of eigenvectors computed by CSTEQR,
*> CPTEQR, and CSTEIN.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX 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 REAL array
*> \endverbatim
*>
*> \param[in] LRWORK
*> \verbatim
*> LRWORK is INTEGER
*> The number of entries in LRWORK (dimension( ??? )
*> \endverbatim
*>
*> \param[out] RESULT
*> \verbatim
*> RESULT is REAL 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 CLATMR, CLATMS, CHETRD, CUNGTR, CSTEQR, SSTERF,
*> or CUNMC2 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.
*
*> \ingroup complex_eig
*
* =====================================================================
SUBROUTINE CCHKST2STG( 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 --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LDU, LIWORK, LRWORK, LWORK, NOUNIT,
$ NSIZES, NTYPES
REAL THRESH
* ..
* .. Array Arguments ..
LOGICAL DOTYPE( * )
INTEGER ISEED( 4 ), IWORK( * ), NN( * )
REAL D1( * ), D2( * ), D3( * ), D4( * ), D5( * ),
$ RESULT( * ), RWORK( * ), SD( * ), SE( * ),
$ WA1( * ), WA2( * ), WA3( * ), WR( * )
COMPLEX A( LDA, * ), AP( * ), TAU( * ), U( LDU, * ),
$ V( LDU, * ), VP( * ), WORK( * ), Z( LDU, * )
* ..
*
* =====================================================================
*
* .. Parameters ..
REAL ZERO, ONE, TWO, EIGHT, TEN, HUN
PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
$ EIGHT = 8.0E0, TEN = 10.0E0, HUN = 100.0E0 )
COMPLEX CZERO, CONE
PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
$ CONE = ( 1.0E+0, 0.0E+0 ) )
REAL 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
REAL 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 )
REAL DUMMA( 1 )
* ..
* .. External Functions ..
INTEGER ILAENV
REAL SLAMCH, SLARND, SSXT1
EXTERNAL ILAENV, SLAMCH, SLARND, SSXT1
* ..
* .. External Subroutines ..
EXTERNAL SCOPY, SLASUM, SSTEBZ, SSTECH, SSTERF, XERBLA,
$ CCOPY, CHET21, CHETRD, CHPT21, CHPTRD, CLACPY,
$ CLASET, CLATMR, CLATMS, CPTEQR, CSTEDC, CSTEMR,
$ CSTEIN, CSTEQR, CSTT21, CSTT22, CUNGTR,
$ CUPGTR, CHETRD_2STAGE, SLASET
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, REAL, CONJG, 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, 'CHETRD', '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( 'CCHKST2STG', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
$ RETURN
*
* More Important constants
*
UNFL = SLAMCH( 'Safe minimum' )
OVFL = ONE / UNFL
ULP = SLAMCH( 'Epsilon' )*SLAMCH( '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( REAL( 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 / REAL( 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 CLASET( '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 CLATMS( 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 CLATMS( 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 CLATMR( 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 CLATMR( 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 CLATMS( 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 CLATMS( 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 ) = CONJG( 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 CHETRD and CUNGTR to compute S and U from
* upper triangle.
*
CALL CLACPY( 'U', N, N, A, LDA, V, LDU )
*
NTEST = 1
CALL CHETRD( 'U', N, V, LDU, SD, SE, TAU, WORK, LWORK,
$ IINFO )
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CHETRD(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 CLACPY( 'U', N, N, V, LDU, U, LDU )
*
NTEST = 2
CALL CUNGTR( 'U', N, U, LDU, TAU, WORK, LWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CUNGTR(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 CHET21( 2, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V,
$ LDU, TAU, WORK, RWORK, RESULT( 1 ) )
CALL CHET21( 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 SCOPY( N, SD, 1, D1, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
*
CALL CSTEQR( 'N', N, D1, RWORK, WORK, LDU, RWORK( N+1 ),
$ IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CSTEQR(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 SLASET( 'Full', N, 1, ZERO, ZERO, SD, N )
CALL SLASET( 'Full', N, 1, ZERO, ZERO, SE, N )
CALL CLACPY( 'U', N, N, A, LDA, V, LDU )
LH = MAX(1, 4*N)
LW = LWORK - LH
CALL CHETRD_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 SCOPY( N, SD, 1, D2, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
*
NTEST = 3
CALL CSTEQR( 'N', N, D2, RWORK, WORK, LDU, RWORK( N+1 ),
$ IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CSTEQR(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 SLASET( 'Full', N, 1, ZERO, ZERO, SD, N )
CALL SLASET( 'Full', N, 1, ZERO, ZERO, SE, N )
CALL CLACPY( 'L', N, N, A, LDA, V, LDU )
CALL CHETRD_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 SCOPY( N, SD, 1, D3, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
*
NTEST = 4
CALL CSTEQR( 'N', N, D3, RWORK, WORK, LDU, RWORK( N+1 ),
$ IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CSTEQR(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 ) )
*
* 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 CHPTRD and CUPGTR to compute S and U from AP
*
CALL CCOPY( NAP, AP, 1, VP, 1 )
*
NTEST = 5
CALL CHPTRD( 'U', N, VP, SD, SE, TAU, IINFO )
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CHPTRD(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 CUPGTR( 'U', N, VP, TAU, U, LDU, WORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CUPGTR(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 CHPT21( 2, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU,
$ WORK, RWORK, RESULT( 5 ) )
CALL CHPT21( 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 CHPTRD and CUPGTR to compute S and U from AP
*
CALL CCOPY( NAP, AP, 1, VP, 1 )
*
NTEST = 7
CALL CHPTRD( 'L', N, VP, SD, SE, TAU, IINFO )
*
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CHPTRD(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 CUPGTR( 'L', N, VP, TAU, U, LDU, WORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CUPGTR(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 CHPT21( 2, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU,
$ WORK, RWORK, RESULT( 7 ) )
CALL CHPT21( 3, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU,
$ WORK, RWORK, RESULT( 8 ) )
*
* Call CSTEQR to compute D1, D2, and Z, do tests.
*
* Compute D1 and Z
*
CALL SCOPY( N, SD, 1, D1, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDU )
*
NTEST = 9
CALL CSTEQR( 'V', N, D1, RWORK, Z, LDU, RWORK( N+1 ),
$ IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CSTEQR(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 SCOPY( N, SD, 1, D2, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
*
NTEST = 11
CALL CSTEQR( 'N', N, D2, RWORK, WORK, LDU, RWORK( N+1 ),
$ IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CSTEQR(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 SCOPY( N, SD, 1, D3, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
*
NTEST = 12
CALL SSTERF( N, D3, RWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'SSTERF', 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 CSTT21( 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 SSTECH( 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 CPTEQR
* and do tests 14, 15, and 16 .
*
IF( JTYPE.GT.15 ) THEN
*
* Compute D4 and Z4
*
CALL SCOPY( N, SD, 1, D4, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDU )
*
NTEST = 14
CALL CPTEQR( 'V', N, D4, RWORK, Z, LDU, RWORK( N+1 ),
$ IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CPTEQR(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 CSTT21( N, 0, SD, SE, D4, DUMMA, Z, LDU, WORK,
$ RWORK, RESULT( 14 ) )
*
* Compute D5
*
CALL SCOPY( N, SD, 1, D5, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
*
NTEST = 16
CALL CPTEQR( 'N', N, D5, RWORK, Z, LDU, RWORK( N+1 ),
$ IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CPTEQR(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 SSTEBZ 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 SSTEBZ( '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 )'SSTEBZ(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 SSTEBZ( '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 )'SSTEBZ(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( SLARND( 1, ISEED2 ) )
IU = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) )
IF( IU.LT.IL ) THEN
ITEMP = IU
IU = IL
IL = ITEMP
END IF
END IF
*
CALL SSTEBZ( '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 )'SSTEBZ(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 SSTEBZ( '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 )'SSTEBZ(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 = SSXT1( 1, WA2, M2, WA3, M3, ABSTOL, ULP, UNFL )
TEMP2 = SSXT1( 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 CSTEIN to compute eigenvectors corresponding to
* eigenvalues in WA1. (First call SSTEBZ again, to make sure
* it returns these eigenvalues in the correct order.)
*
NTEST = 21
CALL SSTEBZ( '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 )'SSTEBZ(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 CSTEIN( 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 )'CSTEIN', 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 CSTT21( N, 0, SD, SE, WA1, DUMMA, Z, LDU, WORK, RWORK,
$ RESULT( 20 ) )
*
* Call CSTEDC(I) to compute D1 and Z, do tests.
*
* Compute D1 and Z
*
INDE = 1
INDRWK = INDE + N
CALL SCOPY( N, SD, 1, D1, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK( INDE ), 1 )
CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDU )
*
NTEST = 22
CALL CSTEDC( 'I', N, D1, RWORK( INDE ), Z, LDU, WORK, LWEDC,
$ RWORK( INDRWK ), LRWEDC, IWORK, LIWEDC, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CSTEDC(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 CSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, RWORK,
$ RESULT( 22 ) )
*
* Call CSTEDC(V) to compute D1 and Z, do tests.
*
* Compute D1 and Z
*
CALL SCOPY( N, SD, 1, D1, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK( INDE ), 1 )
CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDU )
*
NTEST = 24
CALL CSTEDC( 'V', N, D1, RWORK( INDE ), Z, LDU, WORK, LWEDC,
$ RWORK( INDRWK ), LRWEDC, IWORK, LIWEDC, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CSTEDC(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 CSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, RWORK,
$ RESULT( 24 ) )
*
* Call CSTEDC(N) to compute D2, do tests.
*
* Compute D2
*
CALL SCOPY( N, SD, 1, D2, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK( INDE ), 1 )
CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDU )
*
NTEST = 26
CALL CSTEDC( 'N', N, D2, RWORK( INDE ), Z, LDU, WORK, LWEDC,
$ RWORK( INDRWK ), LRWEDC, IWORK, LIWEDC, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CSTEDC(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 CSTEMR if IEEE compliant
*
IF( ILAENV( 10, 'CSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND.
$ ILAENV( 11, 'CSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN
*
* Call CSTEMR, 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 CSTEMR( '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 )'CSTEMR(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( SLARND( 1, ISEED2 ) )
IU = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) )
IF( IU.LT.IL ) THEN
ITEMP = IU
IU = IL
IL = ITEMP
END IF
*
IF( CRANGE ) THEN
NTEST = 28
ABSTOL = UNFL + UNFL
CALL CSTEMR( '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 )'CSTEMR(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 CSTEMR(V,I) to compute D1 and Z, do tests.
*
* Compute D1 and Z
*
CALL SCOPY( N, SD, 1, D5, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDU )
*
IF( CRANGE ) THEN
NTEST = 29
IL = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) )
IU = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) )
IF( IU.LT.IL ) THEN
ITEMP = IU
IU = IL
IL = ITEMP
END IF
CALL CSTEMR( '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 )'CSTEMR(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 CSTEMR to compute D2, do tests.
*
* Compute D2
*
CALL SCOPY( N, SD, 1, D5, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
*
NTEST = 31
CALL CSTEMR( '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 )'CSTEMR(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 CSTEMR(V,V) to compute D1 and Z, do tests.
*
* Compute D1 and Z
*
CALL SCOPY( N, SD, 1, D5, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
CALL CLASET( '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 CSTEMR( '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 )'CSTEMR(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 CSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
$ M, RWORK, RESULT( 32 ) )
*
* Call CSTEMR to compute D2, do tests.
*
* Compute D2
*
CALL SCOPY( N, SD, 1, D5, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
*
NTEST = 34
CALL CSTEMR( '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 )'CSTEMR(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 CSTEMR(V,A) to compute D1 and Z, do tests.
*
* Compute D1 and Z
*
CALL SCOPY( N, SD, 1, D5, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
*
NTEST = 35
*
CALL CSTEMR( '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 )'CSTEMR(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 CSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, M,
$ RWORK, RESULT( 35 ) )
*
* Call CSTEMR to compute D2, do tests.
*
* Compute D2
*
CALL SCOPY( N, SD, 1, D5, 1 )
IF( N.GT.0 )
$ CALL SCOPY( N-1, SE, 1, RWORK, 1 )
*
NTEST = 37
CALL CSTEMR( '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 )'CSTEMR(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 37
*
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 )'CST'
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.0E0 ) 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 SLASUM( 'CST', NOUNIT, NERRS, NTESTT )
RETURN
*
9999 FORMAT( ' CCHKST2STG: ', 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 CCHKST2STG 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, E10.3 )
*
9987 FORMAT( / 'Test performed: see CCHKST2STG for details.', / )
*
* End of CCHKST2STG
*
END
Reference in New Issue View Git Blame Copy Permalink