LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dchkst.f
Go to the documentation of this file.
00001 *> \brief \b DCHKST
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
00012 *                          NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5,
00013 *                          WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK,
00014 *                          LWORK, IWORK, LIWORK, RESULT, INFO )
00015 * 
00016 *       .. Scalar Arguments ..
00017 *       INTEGER            INFO, LDA, LDU, LIWORK, LWORK, NOUNIT, NSIZES,
00018 *      $                   NTYPES
00019 *       DOUBLE PRECISION   THRESH
00020 *       ..
00021 *       .. Array Arguments ..
00022 *       LOGICAL            DOTYPE( * )
00023 *       INTEGER            ISEED( 4 ), IWORK( * ), NN( * )
00024 *       DOUBLE PRECISION   A( LDA, * ), AP( * ), D1( * ), D2( * ),
00025 *      $                   D3( * ), D4( * ), D5( * ), RESULT( * ),
00026 *      $                   SD( * ), SE( * ), TAU( * ), U( LDU, * ),
00027 *      $                   V( LDU, * ), VP( * ), WA1( * ), WA2( * ),
00028 *      $                   WA3( * ), WORK( * ), WR( * ), Z( LDU, * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> DCHKST  checks the symmetric eigenvalue problem routines.
00038 *>
00039 *>    DSYTRD factors A as  U S U' , where ' means transpose,
00040 *>    S is symmetric tridiagonal, and U is orthogonal.
00041 *>    DSYTRD can use either just the lower or just the upper triangle
00042 *>    of A; DCHKST checks both cases.
00043 *>    U is represented as a product of Householder
00044 *>    transformations, whose vectors are stored in the first
00045 *>    n-1 columns of V, and whose scale factors are in TAU.
00046 *>
00047 *>    DSPTRD does the same as DSYTRD, except that A and V are stored
00048 *>    in "packed" format.
00049 *>
00050 *>    DORGTR constructs the matrix U from the contents of V and TAU.
00051 *>
00052 *>    DOPGTR constructs the matrix U from the contents of VP and TAU.
00053 *>
00054 *>    DSTEQR factors S as  Z D1 Z' , where Z is the orthogonal
00055 *>    matrix of eigenvectors and D1 is a diagonal matrix with
00056 *>    the eigenvalues on the diagonal.  D2 is the matrix of
00057 *>    eigenvalues computed when Z is not computed.
00058 *>
00059 *>    DSTERF computes D3, the matrix of eigenvalues, by the
00060 *>    PWK method, which does not yield eigenvectors.
00061 *>
00062 *>    DPTEQR factors S as  Z4 D4 Z4' , for a
00063 *>    symmetric positive definite tridiagonal matrix.
00064 *>    D5 is the matrix of eigenvalues computed when Z is not
00065 *>    computed.
00066 *>
00067 *>    DSTEBZ computes selected eigenvalues.  WA1, WA2, and
00068 *>    WA3 will denote eigenvalues computed to high
00069 *>    absolute accuracy, with different range options.
00070 *>    WR will denote eigenvalues computed to high relative
00071 *>    accuracy.
00072 *>
00073 *>    DSTEIN computes Y, the eigenvectors of S, given the
00074 *>    eigenvalues.
00075 *>
00076 *>    DSTEDC factors S as Z D1 Z' , where Z is the orthogonal
00077 *>    matrix of eigenvectors and D1 is a diagonal matrix with
00078 *>    the eigenvalues on the diagonal ('I' option). It may also
00079 *>    update an input orthogonal matrix, usually the output
00080 *>    from DSYTRD/DORGTR or DSPTRD/DOPGTR ('V' option). It may
00081 *>    also just compute eigenvalues ('N' option).
00082 *>
00083 *>    DSTEMR factors S as Z D1 Z' , where Z is the orthogonal
00084 *>    matrix of eigenvectors and D1 is a diagonal matrix with
00085 *>    the eigenvalues on the diagonal ('I' option).  DSTEMR
00086 *>    uses the Relatively Robust Representation whenever possible.
00087 *>
00088 *> When DCHKST is called, a number of matrix "sizes" ("n's") and a
00089 *> number of matrix "types" are specified.  For each size ("n")
00090 *> and each type of matrix, one matrix will be generated and used
00091 *> to test the symmetric eigenroutines.  For each matrix, a number
00092 *> of tests will be performed:
00093 *>
00094 *> (1)     | A - V S V' | / ( |A| n ulp ) DSYTRD( UPLO='U', ... )
00095 *>
00096 *> (2)     | I - UV' | / ( n ulp )        DORGTR( UPLO='U', ... )
00097 *>
00098 *> (3)     | A - V S V' | / ( |A| n ulp ) DSYTRD( UPLO='L', ... )
00099 *>
00100 *> (4)     | I - UV' | / ( n ulp )        DORGTR( UPLO='L', ... )
00101 *>
00102 *> (5-8)   Same as 1-4, but for DSPTRD and DOPGTR.
00103 *>
00104 *> (9)     | S - Z D Z' | / ( |S| n ulp ) DSTEQR('V',...)
00105 *>
00106 *> (10)    | I - ZZ' | / ( n ulp )        DSTEQR('V',...)
00107 *>
00108 *> (11)    | D1 - D2 | / ( |D1| ulp )        DSTEQR('N',...)
00109 *>
00110 *> (12)    | D1 - D3 | / ( |D1| ulp )        DSTERF
00111 *>
00112 *> (13)    0 if the true eigenvalues (computed by sturm count)
00113 *>         of S are within THRESH of
00114 *>         those in D1.  2*THRESH if they are not.  (Tested using
00115 *>         DSTECH)
00116 *>
00117 *> For S positive definite,
00118 *>
00119 *> (14)    | S - Z4 D4 Z4' | / ( |S| n ulp ) DPTEQR('V',...)
00120 *>
00121 *> (15)    | I - Z4 Z4' | / ( n ulp )        DPTEQR('V',...)
00122 *>
00123 *> (16)    | D4 - D5 | / ( 100 |D4| ulp )       DPTEQR('N',...)
00124 *>
00125 *> When S is also diagonally dominant by the factor gamma < 1,
00126 *>
00127 *> (17)    max | D4(i) - WR(i) | / ( |D4(i)| omega ) ,
00128 *>          i
00129 *>         omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
00130 *>                                              DSTEBZ( 'A', 'E', ...)
00131 *>
00132 *> (18)    | WA1 - D3 | / ( |D3| ulp )          DSTEBZ( 'A', 'E', ...)
00133 *>
00134 *> (19)    ( max { min | WA2(i)-WA3(j) | } +
00135 *>            i     j
00136 *>           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
00137 *>            i     j
00138 *>                                              DSTEBZ( 'I', 'E', ...)
00139 *>
00140 *> (20)    | S - Y WA1 Y' | / ( |S| n ulp )  DSTEBZ, SSTEIN
00141 *>
00142 *> (21)    | I - Y Y' | / ( n ulp )          DSTEBZ, SSTEIN
00143 *>
00144 *> (22)    | S - Z D Z' | / ( |S| n ulp )    DSTEDC('I')
00145 *>
00146 *> (23)    | I - ZZ' | / ( n ulp )           DSTEDC('I')
00147 *>
00148 *> (24)    | S - Z D Z' | / ( |S| n ulp )    DSTEDC('V')
00149 *>
00150 *> (25)    | I - ZZ' | / ( n ulp )           DSTEDC('V')
00151 *>
00152 *> (26)    | D1 - D2 | / ( |D1| ulp )           DSTEDC('V') and
00153 *>                                              DSTEDC('N')
00154 *>
00155 *> Test 27 is disabled at the moment because DSTEMR does not
00156 *> guarantee high relatvie accuracy.
00157 *>
00158 *> (27)    max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
00159 *>          i
00160 *>         omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
00161 *>                                              DSTEMR('V', 'A')
00162 *>
00163 *> (28)    max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
00164 *>          i
00165 *>         omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
00166 *>                                              DSTEMR('V', 'I')
00167 *>
00168 *> Tests 29 through 34 are disable at present because DSTEMR
00169 *> does not handle partial specturm requests.
00170 *>
00171 *> (29)    | S - Z D Z' | / ( |S| n ulp )    DSTEMR('V', 'I')
00172 *>
00173 *> (30)    | I - ZZ' | / ( n ulp )           DSTEMR('V', 'I')
00174 *>
00175 *> (31)    ( max { min | WA2(i)-WA3(j) | } +
00176 *>            i     j
00177 *>           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
00178 *>            i     j
00179 *>         DSTEMR('N', 'I') vs. SSTEMR('V', 'I')
00180 *>
00181 *> (32)    | S - Z D Z' | / ( |S| n ulp )    DSTEMR('V', 'V')
00182 *>
00183 *> (33)    | I - ZZ' | / ( n ulp )           DSTEMR('V', 'V')
00184 *>
00185 *> (34)    ( max { min | WA2(i)-WA3(j) | } +
00186 *>            i     j
00187 *>           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
00188 *>            i     j
00189 *>         DSTEMR('N', 'V') vs. SSTEMR('V', 'V')
00190 *>
00191 *> (35)    | S - Z D Z' | / ( |S| n ulp )    DSTEMR('V', 'A')
00192 *>
00193 *> (36)    | I - ZZ' | / ( n ulp )           DSTEMR('V', 'A')
00194 *>
00195 *> (37)    ( max { min | WA2(i)-WA3(j) | } +
00196 *>            i     j
00197 *>           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
00198 *>            i     j
00199 *>         DSTEMR('N', 'A') vs. SSTEMR('V', 'A')
00200 *>
00201 *> The "sizes" are specified by an array NN(1:NSIZES); the value of
00202 *> each element NN(j) specifies one size.
00203 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
00204 *> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
00205 *> Currently, the list of possible types is:
00206 *>
00207 *> (1)  The zero matrix.
00208 *> (2)  The identity matrix.
00209 *>
00210 *> (3)  A diagonal matrix with evenly spaced entries
00211 *>      1, ..., ULP  and random signs.
00212 *>      (ULP = (first number larger than 1) - 1 )
00213 *> (4)  A diagonal matrix with geometrically spaced entries
00214 *>      1, ..., ULP  and random signs.
00215 *> (5)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
00216 *>      and random signs.
00217 *>
00218 *> (6)  Same as (4), but multiplied by SQRT( overflow threshold )
00219 *> (7)  Same as (4), but multiplied by SQRT( underflow threshold )
00220 *>
00221 *> (8)  A matrix of the form  U' D U, where U is orthogonal and
00222 *>      D has evenly spaced entries 1, ..., ULP with random signs
00223 *>      on the diagonal.
00224 *>
00225 *> (9)  A matrix of the form  U' D U, where U is orthogonal and
00226 *>      D has geometrically spaced entries 1, ..., ULP with random
00227 *>      signs on the diagonal.
00228 *>
00229 *> (10) A matrix of the form  U' D U, where U is orthogonal and
00230 *>      D has "clustered" entries 1, ULP,..., ULP with random
00231 *>      signs on the diagonal.
00232 *>
00233 *> (11) Same as (8), but multiplied by SQRT( overflow threshold )
00234 *> (12) Same as (8), but multiplied by SQRT( underflow threshold )
00235 *>
00236 *> (13) Symmetric matrix with random entries chosen from (-1,1).
00237 *> (14) Same as (13), but multiplied by SQRT( overflow threshold )
00238 *> (15) Same as (13), but multiplied by SQRT( underflow threshold )
00239 *> (16) Same as (8), but diagonal elements are all positive.
00240 *> (17) Same as (9), but diagonal elements are all positive.
00241 *> (18) Same as (10), but diagonal elements are all positive.
00242 *> (19) Same as (16), but multiplied by SQRT( overflow threshold )
00243 *> (20) Same as (16), but multiplied by SQRT( underflow threshold )
00244 *> (21) A diagonally dominant tridiagonal matrix with geometrically
00245 *>      spaced diagonal entries 1, ..., ULP.
00246 *> \endverbatim
00247 *
00248 *  Arguments:
00249 *  ==========
00250 *
00251 *> \param[in] NSIZES
00252 *> \verbatim
00253 *>          NSIZES is INTEGER
00254 *>          The number of sizes of matrices to use.  If it is zero,
00255 *>          DCHKST does nothing.  It must be at least zero.
00256 *> \endverbatim
00257 *>
00258 *> \param[in] NN
00259 *> \verbatim
00260 *>          NN is INTEGER array, dimension (NSIZES)
00261 *>          An array containing the sizes to be used for the matrices.
00262 *>          Zero values will be skipped.  The values must be at least
00263 *>          zero.
00264 *> \endverbatim
00265 *>
00266 *> \param[in] NTYPES
00267 *> \verbatim
00268 *>          NTYPES is INTEGER
00269 *>          The number of elements in DOTYPE.   If it is zero, DCHKST
00270 *>          does nothing.  It must be at least zero.  If it is MAXTYP+1
00271 *>          and NSIZES is 1, then an additional type, MAXTYP+1 is
00272 *>          defined, which is to use whatever matrix is in A.  This
00273 *>          is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
00274 *>          DOTYPE(MAXTYP+1) is .TRUE. .
00275 *> \endverbatim
00276 *>
00277 *> \param[in] DOTYPE
00278 *> \verbatim
00279 *>          DOTYPE is LOGICAL array, dimension (NTYPES)
00280 *>          If DOTYPE(j) is .TRUE., then for each size in NN a
00281 *>          matrix of that size and of type j will be generated.
00282 *>          If NTYPES is smaller than the maximum number of types
00283 *>          defined (PARAMETER MAXTYP), then types NTYPES+1 through
00284 *>          MAXTYP will not be generated.  If NTYPES is larger
00285 *>          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
00286 *>          will be ignored.
00287 *> \endverbatim
00288 *>
00289 *> \param[in,out] ISEED
00290 *> \verbatim
00291 *>          ISEED is INTEGER array, dimension (4)
00292 *>          On entry ISEED specifies the seed of the random number
00293 *>          generator. The array elements should be between 0 and 4095;
00294 *>          if not they will be reduced mod 4096.  Also, ISEED(4) must
00295 *>          be odd.  The random number generator uses a linear
00296 *>          congruential sequence limited to small integers, and so
00297 *>          should produce machine independent random numbers. The
00298 *>          values of ISEED are changed on exit, and can be used in the
00299 *>          next call to DCHKST to continue the same random number
00300 *>          sequence.
00301 *> \endverbatim
00302 *>
00303 *> \param[in] THRESH
00304 *> \verbatim
00305 *>          THRESH is DOUBLE PRECISION
00306 *>          A test will count as "failed" if the "error", computed as
00307 *>          described above, exceeds THRESH.  Note that the error
00308 *>          is scaled to be O(1), so THRESH should be a reasonably
00309 *>          small multiple of 1, e.g., 10 or 100.  In particular,
00310 *>          it should not depend on the precision (single vs. double)
00311 *>          or the size of the matrix.  It must be at least zero.
00312 *> \endverbatim
00313 *>
00314 *> \param[in] NOUNIT
00315 *> \verbatim
00316 *>          NOUNIT is INTEGER
00317 *>          The FORTRAN unit number for printing out error messages
00318 *>          (e.g., if a routine returns IINFO not equal to 0.)
00319 *> \endverbatim
00320 *>
00321 *> \param[in,out] A
00322 *> \verbatim
00323 *>          A is DOUBLE PRECISION array of
00324 *>                                  dimension ( LDA , max(NN) )
00325 *>          Used to hold the matrix whose eigenvalues are to be
00326 *>          computed.  On exit, A contains the last matrix actually
00327 *>          used.
00328 *> \endverbatim
00329 *>
00330 *> \param[in] LDA
00331 *> \verbatim
00332 *>          LDA is INTEGER
00333 *>          The leading dimension of A.  It must be at
00334 *>          least 1 and at least max( NN ).
00335 *> \endverbatim
00336 *>
00337 *> \param[out] AP
00338 *> \verbatim
00339 *>          AP is DOUBLE PRECISION array of
00340 *>                      dimension( max(NN)*max(NN+1)/2 )
00341 *>          The matrix A stored in packed format.
00342 *> \endverbatim
00343 *>
00344 *> \param[out] SD
00345 *> \verbatim
00346 *>          SD is DOUBLE PRECISION array of
00347 *>                             dimension( max(NN) )
00348 *>          The diagonal of the tridiagonal matrix computed by DSYTRD.
00349 *>          On exit, SD and SE contain the tridiagonal form of the
00350 *>          matrix in A.
00351 *> \endverbatim
00352 *>
00353 *> \param[out] SE
00354 *> \verbatim
00355 *>          SE is DOUBLE PRECISION array of
00356 *>                             dimension( max(NN) )
00357 *>          The off-diagonal of the tridiagonal matrix computed by
00358 *>          DSYTRD.  On exit, SD and SE contain the tridiagonal form of
00359 *>          the matrix in A.
00360 *> \endverbatim
00361 *>
00362 *> \param[out] D1
00363 *> \verbatim
00364 *>          D1 is DOUBLE PRECISION array of
00365 *>                             dimension( max(NN) )
00366 *>          The eigenvalues of A, as computed by DSTEQR simlutaneously
00367 *>          with Z.  On exit, the eigenvalues in D1 correspond with the
00368 *>          matrix in A.
00369 *> \endverbatim
00370 *>
00371 *> \param[out] D2
00372 *> \verbatim
00373 *>          D2 is DOUBLE PRECISION array of
00374 *>                             dimension( max(NN) )
00375 *>          The eigenvalues of A, as computed by DSTEQR if Z is not
00376 *>          computed.  On exit, the eigenvalues in D2 correspond with
00377 *>          the matrix in A.
00378 *> \endverbatim
00379 *>
00380 *> \param[out] D3
00381 *> \verbatim
00382 *>          D3 is DOUBLE PRECISION array of
00383 *>                             dimension( max(NN) )
00384 *>          The eigenvalues of A, as computed by DSTERF.  On exit, the
00385 *>          eigenvalues in D3 correspond with the matrix in A.
00386 *> \endverbatim
00387 *>
00388 *> \param[out] D4
00389 *> \verbatim
00390 *>          D4 is DOUBLE PRECISION array of
00391 *>                             dimension( max(NN) )
00392 *>          The eigenvalues of A, as computed by DPTEQR(V).
00393 *>          DPTEQR factors S as  Z4 D4 Z4*
00394 *>          On exit, the eigenvalues in D4 correspond with the matrix in A.
00395 *> \endverbatim
00396 *>
00397 *> \param[out] D5
00398 *> \verbatim
00399 *>          D5 is DOUBLE PRECISION array of
00400 *>                             dimension( max(NN) )
00401 *>          The eigenvalues of A, as computed by DPTEQR(N)
00402 *>          when Z is not computed. On exit, the
00403 *>          eigenvalues in D4 correspond with the matrix in A.
00404 *> \endverbatim
00405 *>
00406 *> \param[out] WA1
00407 *> \verbatim
00408 *>          WA1 is DOUBLE PRECISION array of
00409 *>                             dimension( max(NN) )
00410 *>          All eigenvalues of A, computed to high
00411 *>          absolute accuracy, with different range options.
00412 *>          as computed by DSTEBZ.
00413 *> \endverbatim
00414 *>
00415 *> \param[out] WA2
00416 *> \verbatim
00417 *>          WA2 is DOUBLE PRECISION array of
00418 *>                             dimension( max(NN) )
00419 *>          Selected eigenvalues of A, computed to high
00420 *>          absolute accuracy, with different range options.
00421 *>          as computed by DSTEBZ.
00422 *>          Choose random values for IL and IU, and ask for the
00423 *>          IL-th through IU-th eigenvalues.
00424 *> \endverbatim
00425 *>
00426 *> \param[out] WA3
00427 *> \verbatim
00428 *>          WA3 is DOUBLE PRECISION array of
00429 *>                             dimension( max(NN) )
00430 *>          Selected eigenvalues of A, computed to high
00431 *>          absolute accuracy, with different range options.
00432 *>          as computed by DSTEBZ.
00433 *>          Determine the values VL and VU of the IL-th and IU-th
00434 *>          eigenvalues and ask for all eigenvalues in this range.
00435 *> \endverbatim
00436 *>
00437 *> \param[out] WR
00438 *> \verbatim
00439 *>          WR is DOUBLE PRECISION array of
00440 *>                             dimension( max(NN) )
00441 *>          All eigenvalues of A, computed to high
00442 *>          absolute accuracy, with different options.
00443 *>          as computed by DSTEBZ.
00444 *> \endverbatim
00445 *>
00446 *> \param[out] U
00447 *> \verbatim
00448 *>          U is DOUBLE PRECISION array of
00449 *>                             dimension( LDU, max(NN) ).
00450 *>          The orthogonal matrix computed by DSYTRD + DORGTR.
00451 *> \endverbatim
00452 *>
00453 *> \param[in] LDU
00454 *> \verbatim
00455 *>          LDU is INTEGER
00456 *>          The leading dimension of U, Z, and V.  It must be at least 1
00457 *>          and at least max( NN ).
00458 *> \endverbatim
00459 *>
00460 *> \param[out] V
00461 *> \verbatim
00462 *>          V is DOUBLE PRECISION array of
00463 *>                             dimension( LDU, max(NN) ).
00464 *>          The Housholder vectors computed by DSYTRD in reducing A to
00465 *>          tridiagonal form.  The vectors computed with UPLO='U' are
00466 *>          in the upper triangle, and the vectors computed with UPLO='L'
00467 *>          are in the lower triangle.  (As described in DSYTRD, the
00468 *>          sub- and superdiagonal are not set to 1, although the
00469 *>          true Householder vector has a 1 in that position.  The
00470 *>          routines that use V, such as DORGTR, set those entries to
00471 *>          1 before using them, and then restore them later.)
00472 *> \endverbatim
00473 *>
00474 *> \param[out] VP
00475 *> \verbatim
00476 *>          VP is DOUBLE PRECISION array of
00477 *>                      dimension( max(NN)*max(NN+1)/2 )
00478 *>          The matrix V stored in packed format.
00479 *> \endverbatim
00480 *>
00481 *> \param[out] TAU
00482 *> \verbatim
00483 *>          TAU is DOUBLE PRECISION array of
00484 *>                             dimension( max(NN) )
00485 *>          The Householder factors computed by DSYTRD in reducing A
00486 *>          to tridiagonal form.
00487 *> \endverbatim
00488 *>
00489 *> \param[out] Z
00490 *> \verbatim
00491 *>          Z is DOUBLE PRECISION array of
00492 *>                             dimension( LDU, max(NN) ).
00493 *>          The orthogonal matrix of eigenvectors computed by DSTEQR,
00494 *>          DPTEQR, and DSTEIN.
00495 *> \endverbatim
00496 *>
00497 *> \param[out] WORK
00498 *> \verbatim
00499 *>          WORK is DOUBLE PRECISION array of
00500 *>                      dimension( LWORK )
00501 *> \endverbatim
00502 *>
00503 *> \param[in] LWORK
00504 *> \verbatim
00505 *>          LWORK is INTEGER
00506 *>          The number of entries in WORK.  This must be at least
00507 *>          1 + 4 * Nmax + 2 * Nmax * lg Nmax + 3 * Nmax**2
00508 *>          where Nmax = max( NN(j), 2 ) and lg = log base 2.
00509 *> \endverbatim
00510 *>
00511 *> \param[out] IWORK
00512 *> \verbatim
00513 *>          IWORK is INTEGER array,
00514 *>          Workspace.
00515 *> \endverbatim
00516 *>
00517 *> \param[out] LIWORK
00518 *> \verbatim
00519 *>          LIWORK is INTEGER
00520 *>          The number of entries in IWORK.  This must be at least
00521 *>                  6 + 6*Nmax + 5 * Nmax * lg Nmax 
00522 *>          where Nmax = max( NN(j), 2 ) and lg = log base 2.
00523 *> \endverbatim
00524 *>
00525 *> \param[out] RESULT
00526 *> \verbatim
00527 *>          RESULT is DOUBLE PRECISION array, dimension (26)
00528 *>          The values computed by the tests described above.
00529 *>          The values are currently limited to 1/ulp, to avoid
00530 *>          overflow.
00531 *> \endverbatim
00532 *>
00533 *> \param[out] INFO
00534 *> \verbatim
00535 *>          INFO is INTEGER
00536 *>          If 0, then everything ran OK.
00537 *>           -1: NSIZES < 0
00538 *>           -2: Some NN(j) < 0
00539 *>           -3: NTYPES < 0
00540 *>           -5: THRESH < 0
00541 *>           -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
00542 *>          -23: LDU < 1 or LDU < NMAX.
00543 *>          -29: LWORK too small.
00544 *>          If  DLATMR, SLATMS, DSYTRD, DORGTR, DSTEQR, SSTERF,
00545 *>              or DORMC2 returns an error code, the
00546 *>              absolute value of it is returned.
00547 *>
00548 *>-----------------------------------------------------------------------
00549 *>
00550 *>       Some Local Variables and Parameters:
00551 *>       ---- ----- --------- --- ----------
00552 *>       ZERO, ONE       Real 0 and 1.
00553 *>       MAXTYP          The number of types defined.
00554 *>       NTEST           The number of tests performed, or which can
00555 *>                       be performed so far, for the current matrix.
00556 *>       NTESTT          The total number of tests performed so far.
00557 *>       NBLOCK          Blocksize as returned by ENVIR.
00558 *>       NMAX            Largest value in NN.
00559 *>       NMATS           The number of matrices generated so far.
00560 *>       NERRS           The number of tests which have exceeded THRESH
00561 *>                       so far.
00562 *>       COND, IMODE     Values to be passed to the matrix generators.
00563 *>       ANORM           Norm of A; passed to matrix generators.
00564 *>
00565 *>       OVFL, UNFL      Overflow and underflow thresholds.
00566 *>       ULP, ULPINV     Finest relative precision and its inverse.
00567 *>       RTOVFL, RTUNFL  Square roots of the previous 2 values.
00568 *>               The following four arrays decode JTYPE:
00569 *>       KTYPE(j)        The general type (1-10) for type "j".
00570 *>       KMODE(j)        The MODE value to be passed to the matrix
00571 *>                       generator for type "j".
00572 *>       KMAGN(j)        The order of magnitude ( O(1),
00573 *>                       O(overflow^(1/2) ), O(underflow^(1/2) )
00574 *> \endverbatim
00575 *
00576 *  Authors:
00577 *  ========
00578 *
00579 *> \author Univ. of Tennessee 
00580 *> \author Univ. of California Berkeley 
00581 *> \author Univ. of Colorado Denver 
00582 *> \author NAG Ltd. 
00583 *
00584 *> \date November 2011
00585 *
00586 *> \ingroup double_eig
00587 *
00588 *  =====================================================================
00589       SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
00590      $                   NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5,
00591      $                   WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK,
00592      $                   LWORK, IWORK, LIWORK, RESULT, INFO )
00593 *
00594 *  -- LAPACK test routine (version 3.4.0) --
00595 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00596 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00597 *     November 2011
00598 *
00599 *     .. Scalar Arguments ..
00600       INTEGER            INFO, LDA, LDU, LIWORK, LWORK, NOUNIT, NSIZES,
00601      $                   NTYPES
00602       DOUBLE PRECISION   THRESH
00603 *     ..
00604 *     .. Array Arguments ..
00605       LOGICAL            DOTYPE( * )
00606       INTEGER            ISEED( 4 ), IWORK( * ), NN( * )
00607       DOUBLE PRECISION   A( LDA, * ), AP( * ), D1( * ), D2( * ),
00608      $                   D3( * ), D4( * ), D5( * ), RESULT( * ),
00609      $                   SD( * ), SE( * ), TAU( * ), U( LDU, * ),
00610      $                   V( LDU, * ), VP( * ), WA1( * ), WA2( * ),
00611      $                   WA3( * ), WORK( * ), WR( * ), Z( LDU, * )
00612 *     ..
00613 *
00614 *  =====================================================================
00615 *
00616 *     .. Parameters ..
00617       DOUBLE PRECISION   ZERO, ONE, TWO, EIGHT, TEN, HUN
00618       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
00619      $                   EIGHT = 8.0D0, TEN = 10.0D0, HUN = 100.0D0 )
00620       DOUBLE PRECISION   HALF
00621       PARAMETER          ( HALF = ONE / TWO )
00622       INTEGER            MAXTYP
00623       PARAMETER          ( MAXTYP = 21 )
00624       LOGICAL            SRANGE
00625       PARAMETER          ( SRANGE = .FALSE. )
00626       LOGICAL            SREL
00627       PARAMETER          ( SREL = .FALSE. )
00628 *     ..
00629 *     .. Local Scalars ..
00630       LOGICAL            BADNN, TRYRAC
00631       INTEGER            I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, J, JC,
00632      $                   JR, JSIZE, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC,
00633      $                   M, M2, M3, MTYPES, N, NAP, NBLOCK, NERRS,
00634      $                   NMATS, NMAX, NSPLIT, NTEST, NTESTT
00635       DOUBLE PRECISION   ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL,
00636      $                   RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP,
00637      $                   ULPINV, UNFL, VL, VU
00638 *     ..
00639 *     .. Local Arrays ..
00640       INTEGER            IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
00641      $                   KMAGN( MAXTYP ), KMODE( MAXTYP ),
00642      $                   KTYPE( MAXTYP )
00643       DOUBLE PRECISION   DUMMA( 1 )
00644 *     ..
00645 *     .. External Functions ..
00646       INTEGER            ILAENV
00647       DOUBLE PRECISION   DLAMCH, DLARND, DSXT1
00648       EXTERNAL           ILAENV, DLAMCH, DLARND, DSXT1
00649 *     ..
00650 *     .. External Subroutines ..
00651       EXTERNAL           DCOPY, DLABAD, DLACPY, DLASET, DLASUM, DLATMR,
00652      $                   DLATMS, DOPGTR, DORGTR, DPTEQR, DSPT21, DSPTRD,
00653      $                   DSTEBZ, DSTECH, DSTEDC, DSTEMR, DSTEIN, DSTEQR,
00654      $                   DSTERF, DSTT21, DSTT22, DSYT21, DSYTRD, XERBLA
00655 *     ..
00656 *     .. Intrinsic Functions ..
00657       INTRINSIC          ABS, DBLE, INT, LOG, MAX, MIN, SQRT
00658 *     ..
00659 *     .. Data statements ..
00660       DATA               KTYPE / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8,
00661      $                   8, 8, 9, 9, 9, 9, 9, 10 /
00662       DATA               KMAGN / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
00663      $                   2, 3, 1, 1, 1, 2, 3, 1 /
00664       DATA               KMODE / 0, 0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
00665      $                   0, 0, 4, 3, 1, 4, 4, 3 /
00666 *     ..
00667 *     .. Executable Statements ..
00668 *
00669 *     Keep ftnchek happy
00670       IDUMMA( 1 ) = 1
00671 *
00672 *     Check for errors
00673 *
00674       NTESTT = 0
00675       INFO = 0
00676 *
00677 *     Important constants
00678 *
00679       BADNN = .FALSE.
00680       TRYRAC = .TRUE.
00681       NMAX = 1
00682       DO 10 J = 1, NSIZES
00683          NMAX = MAX( NMAX, NN( J ) )
00684          IF( NN( J ).LT.0 )
00685      $      BADNN = .TRUE.
00686    10 CONTINUE
00687 *
00688       NBLOCK = ILAENV( 1, 'DSYTRD', 'L', NMAX, -1, -1, -1 )
00689       NBLOCK = MIN( NMAX, MAX( 1, NBLOCK ) )
00690 *
00691 *     Check for errors
00692 *
00693       IF( NSIZES.LT.0 ) THEN
00694          INFO = -1
00695       ELSE IF( BADNN ) THEN
00696          INFO = -2
00697       ELSE IF( NTYPES.LT.0 ) THEN
00698          INFO = -3
00699       ELSE IF( LDA.LT.NMAX ) THEN
00700          INFO = -9
00701       ELSE IF( LDU.LT.NMAX ) THEN
00702          INFO = -23
00703       ELSE IF( 2*MAX( 2, NMAX )**2.GT.LWORK ) THEN
00704          INFO = -29
00705       END IF
00706 *
00707       IF( INFO.NE.0 ) THEN
00708          CALL XERBLA( 'DCHKST', -INFO )
00709          RETURN
00710       END IF
00711 *
00712 *     Quick return if possible
00713 *
00714       IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
00715      $   RETURN
00716 *
00717 *     More Important constants
00718 *
00719       UNFL = DLAMCH( 'Safe minimum' )
00720       OVFL = ONE / UNFL
00721       CALL DLABAD( UNFL, OVFL )
00722       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00723       ULPINV = ONE / ULP
00724       LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) )
00725       RTUNFL = SQRT( UNFL )
00726       RTOVFL = SQRT( OVFL )
00727 *
00728 *     Loop over sizes, types
00729 *
00730       DO 20 I = 1, 4
00731          ISEED2( I ) = ISEED( I )
00732    20 CONTINUE
00733       NERRS = 0
00734       NMATS = 0
00735 *
00736       DO 310 JSIZE = 1, NSIZES
00737          N = NN( JSIZE )
00738          IF( N.GT.0 ) THEN
00739             LGN = INT( LOG( DBLE( N ) ) / LOG( TWO ) )
00740             IF( 2**LGN.LT.N )
00741      $         LGN = LGN + 1
00742             IF( 2**LGN.LT.N )
00743      $         LGN = LGN + 1
00744             LWEDC = 1 + 4*N + 2*N*LGN + 4*N**2
00745             LIWEDC = 6 + 6*N + 5*N*LGN
00746          ELSE
00747             LWEDC = 8
00748             LIWEDC = 12
00749          END IF
00750          NAP = ( N*( N+1 ) ) / 2
00751          ANINV = ONE / DBLE( MAX( 1, N ) )
00752 *
00753          IF( NSIZES.NE.1 ) THEN
00754             MTYPES = MIN( MAXTYP, NTYPES )
00755          ELSE
00756             MTYPES = MIN( MAXTYP+1, NTYPES )
00757          END IF
00758 *
00759          DO 300 JTYPE = 1, MTYPES
00760             IF( .NOT.DOTYPE( JTYPE ) )
00761      $         GO TO 300
00762             NMATS = NMATS + 1
00763             NTEST = 0
00764 *
00765             DO 30 J = 1, 4
00766                IOLDSD( J ) = ISEED( J )
00767    30       CONTINUE
00768 *
00769 *           Compute "A"
00770 *
00771 *           Control parameters:
00772 *
00773 *               KMAGN  KMODE        KTYPE
00774 *           =1  O(1)   clustered 1  zero
00775 *           =2  large  clustered 2  identity
00776 *           =3  small  exponential  (none)
00777 *           =4         arithmetic   diagonal, (w/ eigenvalues)
00778 *           =5         random log   symmetric, w/ eigenvalues
00779 *           =6         random       (none)
00780 *           =7                      random diagonal
00781 *           =8                      random symmetric
00782 *           =9                      positive definite
00783 *           =10                     diagonally dominant tridiagonal
00784 *
00785             IF( MTYPES.GT.MAXTYP )
00786      $         GO TO 100
00787 *
00788             ITYPE = KTYPE( JTYPE )
00789             IMODE = KMODE( JTYPE )
00790 *
00791 *           Compute norm
00792 *
00793             GO TO ( 40, 50, 60 )KMAGN( JTYPE )
00794 *
00795    40       CONTINUE
00796             ANORM = ONE
00797             GO TO 70
00798 *
00799    50       CONTINUE
00800             ANORM = ( RTOVFL*ULP )*ANINV
00801             GO TO 70
00802 *
00803    60       CONTINUE
00804             ANORM = RTUNFL*N*ULPINV
00805             GO TO 70
00806 *
00807    70       CONTINUE
00808 *
00809             CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
00810             IINFO = 0
00811             IF( JTYPE.LE.15 ) THEN
00812                COND = ULPINV
00813             ELSE
00814                COND = ULPINV*ANINV / TEN
00815             END IF
00816 *
00817 *           Special Matrices -- Identity & Jordan block
00818 *
00819 *              Zero
00820 *
00821             IF( ITYPE.EQ.1 ) THEN
00822                IINFO = 0
00823 *
00824             ELSE IF( ITYPE.EQ.2 ) THEN
00825 *
00826 *              Identity
00827 *
00828                DO 80 JC = 1, N
00829                   A( JC, JC ) = ANORM
00830    80          CONTINUE
00831 *
00832             ELSE IF( ITYPE.EQ.4 ) THEN
00833 *
00834 *              Diagonal Matrix, [Eigen]values Specified
00835 *
00836                CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
00837      $                      ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
00838      $                      IINFO )
00839 *
00840 *
00841             ELSE IF( ITYPE.EQ.5 ) THEN
00842 *
00843 *              Symmetric, eigenvalues specified
00844 *
00845                CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
00846      $                      ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
00847      $                      IINFO )
00848 *
00849             ELSE IF( ITYPE.EQ.7 ) THEN
00850 *
00851 *              Diagonal, random eigenvalues
00852 *
00853                CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
00854      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00855      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
00856      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00857 *
00858             ELSE IF( ITYPE.EQ.8 ) THEN
00859 *
00860 *              Symmetric, random eigenvalues
00861 *
00862                CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
00863      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00864      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
00865      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00866 *
00867             ELSE IF( ITYPE.EQ.9 ) THEN
00868 *
00869 *              Positive definite, eigenvalues specified.
00870 *
00871                CALL DLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND,
00872      $                      ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
00873      $                      IINFO )
00874 *
00875             ELSE IF( ITYPE.EQ.10 ) THEN
00876 *
00877 *              Positive definite tridiagonal, eigenvalues specified.
00878 *
00879                CALL DLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND,
00880      $                      ANORM, 1, 1, 'N', A, LDA, WORK( N+1 ),
00881      $                      IINFO )
00882                DO 90 I = 2, N
00883                   TEMP1 = ABS( A( I-1, I ) ) /
00884      $                    SQRT( ABS( A( I-1, I-1 )*A( I, I ) ) )
00885                   IF( TEMP1.GT.HALF ) THEN
00886                      A( I-1, I ) = HALF*SQRT( ABS( A( I-1, I-1 )*A( I,
00887      $                             I ) ) )
00888                      A( I, I-1 ) = A( I-1, I )
00889                   END IF
00890    90          CONTINUE
00891 *
00892             ELSE
00893 *
00894                IINFO = 1
00895             END IF
00896 *
00897             IF( IINFO.NE.0 ) THEN
00898                WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
00899      $            IOLDSD
00900                INFO = ABS( IINFO )
00901                RETURN
00902             END IF
00903 *
00904   100       CONTINUE
00905 *
00906 *           Call DSYTRD and DORGTR to compute S and U from
00907 *           upper triangle.
00908 *
00909             CALL DLACPY( 'U', N, N, A, LDA, V, LDU )
00910 *
00911             NTEST = 1
00912             CALL DSYTRD( 'U', N, V, LDU, SD, SE, TAU, WORK, LWORK,
00913      $                   IINFO )
00914 *
00915             IF( IINFO.NE.0 ) THEN
00916                WRITE( NOUNIT, FMT = 9999 )'DSYTRD(U)', IINFO, N, JTYPE,
00917      $            IOLDSD
00918                INFO = ABS( IINFO )
00919                IF( IINFO.LT.0 ) THEN
00920                   RETURN
00921                ELSE
00922                   RESULT( 1 ) = ULPINV
00923                   GO TO 280
00924                END IF
00925             END IF
00926 *
00927             CALL DLACPY( 'U', N, N, V, LDU, U, LDU )
00928 *
00929             NTEST = 2
00930             CALL DORGTR( 'U', N, U, LDU, TAU, WORK, LWORK, IINFO )
00931             IF( IINFO.NE.0 ) THEN
00932                WRITE( NOUNIT, FMT = 9999 )'DORGTR(U)', IINFO, N, JTYPE,
00933      $            IOLDSD
00934                INFO = ABS( IINFO )
00935                IF( IINFO.LT.0 ) THEN
00936                   RETURN
00937                ELSE
00938                   RESULT( 2 ) = ULPINV
00939                   GO TO 280
00940                END IF
00941             END IF
00942 *
00943 *           Do tests 1 and 2
00944 *
00945             CALL DSYT21( 2, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V,
00946      $                   LDU, TAU, WORK, RESULT( 1 ) )
00947             CALL DSYT21( 3, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V,
00948      $                   LDU, TAU, WORK, RESULT( 2 ) )
00949 *
00950 *           Call DSYTRD and DORGTR to compute S and U from
00951 *           lower triangle, do tests.
00952 *
00953             CALL DLACPY( 'L', N, N, A, LDA, V, LDU )
00954 *
00955             NTEST = 3
00956             CALL DSYTRD( 'L', N, V, LDU, SD, SE, TAU, WORK, LWORK,
00957      $                   IINFO )
00958 *
00959             IF( IINFO.NE.0 ) THEN
00960                WRITE( NOUNIT, FMT = 9999 )'DSYTRD(L)', IINFO, N, JTYPE,
00961      $            IOLDSD
00962                INFO = ABS( IINFO )
00963                IF( IINFO.LT.0 ) THEN
00964                   RETURN
00965                ELSE
00966                   RESULT( 3 ) = ULPINV
00967                   GO TO 280
00968                END IF
00969             END IF
00970 *
00971             CALL DLACPY( 'L', N, N, V, LDU, U, LDU )
00972 *
00973             NTEST = 4
00974             CALL DORGTR( 'L', N, U, LDU, TAU, WORK, LWORK, IINFO )
00975             IF( IINFO.NE.0 ) THEN
00976                WRITE( NOUNIT, FMT = 9999 )'DORGTR(L)', IINFO, N, JTYPE,
00977      $            IOLDSD
00978                INFO = ABS( IINFO )
00979                IF( IINFO.LT.0 ) THEN
00980                   RETURN
00981                ELSE
00982                   RESULT( 4 ) = ULPINV
00983                   GO TO 280
00984                END IF
00985             END IF
00986 *
00987             CALL DSYT21( 2, 'Lower', N, 1, A, LDA, SD, SE, U, LDU, V,
00988      $                   LDU, TAU, WORK, RESULT( 3 ) )
00989             CALL DSYT21( 3, 'Lower', N, 1, A, LDA, SD, SE, U, LDU, V,
00990      $                   LDU, TAU, WORK, RESULT( 4 ) )
00991 *
00992 *           Store the upper triangle of A in AP
00993 *
00994             I = 0
00995             DO 120 JC = 1, N
00996                DO 110 JR = 1, JC
00997                   I = I + 1
00998                   AP( I ) = A( JR, JC )
00999   110          CONTINUE
01000   120       CONTINUE
01001 *
01002 *           Call DSPTRD and DOPGTR to compute S and U from AP
01003 *
01004             CALL DCOPY( NAP, AP, 1, VP, 1 )
01005 *
01006             NTEST = 5
01007             CALL DSPTRD( 'U', N, VP, SD, SE, TAU, IINFO )
01008 *
01009             IF( IINFO.NE.0 ) THEN
01010                WRITE( NOUNIT, FMT = 9999 )'DSPTRD(U)', IINFO, N, JTYPE,
01011      $            IOLDSD
01012                INFO = ABS( IINFO )
01013                IF( IINFO.LT.0 ) THEN
01014                   RETURN
01015                ELSE
01016                   RESULT( 5 ) = ULPINV
01017                   GO TO 280
01018                END IF
01019             END IF
01020 *
01021             NTEST = 6
01022             CALL DOPGTR( 'U', N, VP, TAU, U, LDU, WORK, IINFO )
01023             IF( IINFO.NE.0 ) THEN
01024                WRITE( NOUNIT, FMT = 9999 )'DOPGTR(U)', IINFO, N, JTYPE,
01025      $            IOLDSD
01026                INFO = ABS( IINFO )
01027                IF( IINFO.LT.0 ) THEN
01028                   RETURN
01029                ELSE
01030                   RESULT( 6 ) = ULPINV
01031                   GO TO 280
01032                END IF
01033             END IF
01034 *
01035 *           Do tests 5 and 6
01036 *
01037             CALL DSPT21( 2, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU,
01038      $                   WORK, RESULT( 5 ) )
01039             CALL DSPT21( 3, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU,
01040      $                   WORK, RESULT( 6 ) )
01041 *
01042 *           Store the lower triangle of A in AP
01043 *
01044             I = 0
01045             DO 140 JC = 1, N
01046                DO 130 JR = JC, N
01047                   I = I + 1
01048                   AP( I ) = A( JR, JC )
01049   130          CONTINUE
01050   140       CONTINUE
01051 *
01052 *           Call DSPTRD and DOPGTR to compute S and U from AP
01053 *
01054             CALL DCOPY( NAP, AP, 1, VP, 1 )
01055 *
01056             NTEST = 7
01057             CALL DSPTRD( 'L', N, VP, SD, SE, TAU, IINFO )
01058 *
01059             IF( IINFO.NE.0 ) THEN
01060                WRITE( NOUNIT, FMT = 9999 )'DSPTRD(L)', IINFO, N, JTYPE,
01061      $            IOLDSD
01062                INFO = ABS( IINFO )
01063                IF( IINFO.LT.0 ) THEN
01064                   RETURN
01065                ELSE
01066                   RESULT( 7 ) = ULPINV
01067                   GO TO 280
01068                END IF
01069             END IF
01070 *
01071             NTEST = 8
01072             CALL DOPGTR( 'L', N, VP, TAU, U, LDU, WORK, IINFO )
01073             IF( IINFO.NE.0 ) THEN
01074                WRITE( NOUNIT, FMT = 9999 )'DOPGTR(L)', IINFO, N, JTYPE,
01075      $            IOLDSD
01076                INFO = ABS( IINFO )
01077                IF( IINFO.LT.0 ) THEN
01078                   RETURN
01079                ELSE
01080                   RESULT( 8 ) = ULPINV
01081                   GO TO 280
01082                END IF
01083             END IF
01084 *
01085             CALL DSPT21( 2, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU,
01086      $                   WORK, RESULT( 7 ) )
01087             CALL DSPT21( 3, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU,
01088      $                   WORK, RESULT( 8 ) )
01089 *
01090 *           Call DSTEQR to compute D1, D2, and Z, do tests.
01091 *
01092 *           Compute D1 and Z
01093 *
01094             CALL DCOPY( N, SD, 1, D1, 1 )
01095             IF( N.GT.0 )
01096      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
01097             CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
01098 *
01099             NTEST = 9
01100             CALL DSTEQR( 'V', N, D1, WORK, Z, LDU, WORK( N+1 ), IINFO )
01101             IF( IINFO.NE.0 ) THEN
01102                WRITE( NOUNIT, FMT = 9999 )'DSTEQR(V)', IINFO, N, JTYPE,
01103      $            IOLDSD
01104                INFO = ABS( IINFO )
01105                IF( IINFO.LT.0 ) THEN
01106                   RETURN
01107                ELSE
01108                   RESULT( 9 ) = ULPINV
01109                   GO TO 280
01110                END IF
01111             END IF
01112 *
01113 *           Compute D2
01114 *
01115             CALL DCOPY( N, SD, 1, D2, 1 )
01116             IF( N.GT.0 )
01117      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
01118 *
01119             NTEST = 11
01120             CALL DSTEQR( 'N', N, D2, WORK, WORK( N+1 ), LDU,
01121      $                   WORK( N+1 ), IINFO )
01122             IF( IINFO.NE.0 ) THEN
01123                WRITE( NOUNIT, FMT = 9999 )'DSTEQR(N)', IINFO, N, JTYPE,
01124      $            IOLDSD
01125                INFO = ABS( IINFO )
01126                IF( IINFO.LT.0 ) THEN
01127                   RETURN
01128                ELSE
01129                   RESULT( 11 ) = ULPINV
01130                   GO TO 280
01131                END IF
01132             END IF
01133 *
01134 *           Compute D3 (using PWK method)
01135 *
01136             CALL DCOPY( N, SD, 1, D3, 1 )
01137             IF( N.GT.0 )
01138      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
01139 *
01140             NTEST = 12
01141             CALL DSTERF( N, D3, WORK, IINFO )
01142             IF( IINFO.NE.0 ) THEN
01143                WRITE( NOUNIT, FMT = 9999 )'DSTERF', IINFO, N, JTYPE,
01144      $            IOLDSD
01145                INFO = ABS( IINFO )
01146                IF( IINFO.LT.0 ) THEN
01147                   RETURN
01148                ELSE
01149                   RESULT( 12 ) = ULPINV
01150                   GO TO 280
01151                END IF
01152             END IF
01153 *
01154 *           Do Tests 9 and 10
01155 *
01156             CALL DSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
01157      $                   RESULT( 9 ) )
01158 *
01159 *           Do Tests 11 and 12
01160 *
01161             TEMP1 = ZERO
01162             TEMP2 = ZERO
01163             TEMP3 = ZERO
01164             TEMP4 = ZERO
01165 *
01166             DO 150 J = 1, N
01167                TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
01168                TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
01169                TEMP3 = MAX( TEMP3, ABS( D1( J ) ), ABS( D3( J ) ) )
01170                TEMP4 = MAX( TEMP4, ABS( D1( J )-D3( J ) ) )
01171   150       CONTINUE
01172 *
01173             RESULT( 11 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
01174             RESULT( 12 ) = TEMP4 / MAX( UNFL, ULP*MAX( TEMP3, TEMP4 ) )
01175 *
01176 *           Do Test 13 -- Sturm Sequence Test of Eigenvalues
01177 *                         Go up by factors of two until it succeeds
01178 *
01179             NTEST = 13
01180             TEMP1 = THRESH*( HALF-ULP )
01181 *
01182             DO 160 J = 0, LOG2UI
01183                CALL DSTECH( N, SD, SE, D1, TEMP1, WORK, IINFO )
01184                IF( IINFO.EQ.0 )
01185      $            GO TO 170
01186                TEMP1 = TEMP1*TWO
01187   160       CONTINUE
01188 *
01189   170       CONTINUE
01190             RESULT( 13 ) = TEMP1
01191 *
01192 *           For positive definite matrices ( JTYPE.GT.15 ) call DPTEQR
01193 *           and do tests 14, 15, and 16 .
01194 *
01195             IF( JTYPE.GT.15 ) THEN
01196 *
01197 *              Compute D4 and Z4
01198 *
01199                CALL DCOPY( N, SD, 1, D4, 1 )
01200                IF( N.GT.0 )
01201      $            CALL DCOPY( N-1, SE, 1, WORK, 1 )
01202                CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
01203 *
01204                NTEST = 14
01205                CALL DPTEQR( 'V', N, D4, WORK, Z, LDU, WORK( N+1 ),
01206      $                      IINFO )
01207                IF( IINFO.NE.0 ) THEN
01208                   WRITE( NOUNIT, FMT = 9999 )'DPTEQR(V)', IINFO, N,
01209      $               JTYPE, IOLDSD
01210                   INFO = ABS( IINFO )
01211                   IF( IINFO.LT.0 ) THEN
01212                      RETURN
01213                   ELSE
01214                      RESULT( 14 ) = ULPINV
01215                      GO TO 280
01216                   END IF
01217                END IF
01218 *
01219 *              Do Tests 14 and 15
01220 *
01221                CALL DSTT21( N, 0, SD, SE, D4, DUMMA, Z, LDU, WORK,
01222      $                      RESULT( 14 ) )
01223 *
01224 *              Compute D5
01225 *
01226                CALL DCOPY( N, SD, 1, D5, 1 )
01227                IF( N.GT.0 )
01228      $            CALL DCOPY( N-1, SE, 1, WORK, 1 )
01229 *
01230                NTEST = 16
01231                CALL DPTEQR( 'N', N, D5, WORK, Z, LDU, WORK( N+1 ),
01232      $                      IINFO )
01233                IF( IINFO.NE.0 ) THEN
01234                   WRITE( NOUNIT, FMT = 9999 )'DPTEQR(N)', IINFO, N,
01235      $               JTYPE, IOLDSD
01236                   INFO = ABS( IINFO )
01237                   IF( IINFO.LT.0 ) THEN
01238                      RETURN
01239                   ELSE
01240                      RESULT( 16 ) = ULPINV
01241                      GO TO 280
01242                   END IF
01243                END IF
01244 *
01245 *              Do Test 16
01246 *
01247                TEMP1 = ZERO
01248                TEMP2 = ZERO
01249                DO 180 J = 1, N
01250                   TEMP1 = MAX( TEMP1, ABS( D4( J ) ), ABS( D5( J ) ) )
01251                   TEMP2 = MAX( TEMP2, ABS( D4( J )-D5( J ) ) )
01252   180          CONTINUE
01253 *
01254                RESULT( 16 ) = TEMP2 / MAX( UNFL,
01255      $                        HUN*ULP*MAX( TEMP1, TEMP2 ) )
01256             ELSE
01257                RESULT( 14 ) = ZERO
01258                RESULT( 15 ) = ZERO
01259                RESULT( 16 ) = ZERO
01260             END IF
01261 *
01262 *           Call DSTEBZ with different options and do tests 17-18.
01263 *
01264 *              If S is positive definite and diagonally dominant,
01265 *              ask for all eigenvalues with high relative accuracy.
01266 *
01267             VL = ZERO
01268             VU = ZERO
01269             IL = 0
01270             IU = 0
01271             IF( JTYPE.EQ.21 ) THEN
01272                NTEST = 17
01273                ABSTOL = UNFL + UNFL
01274                CALL DSTEBZ( 'A', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE,
01275      $                      M, NSPLIT, WR, IWORK( 1 ), IWORK( N+1 ),
01276      $                      WORK, IWORK( 2*N+1 ), IINFO )
01277                IF( IINFO.NE.0 ) THEN
01278                   WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(A,rel)', IINFO, N,
01279      $               JTYPE, IOLDSD
01280                   INFO = ABS( IINFO )
01281                   IF( IINFO.LT.0 ) THEN
01282                      RETURN
01283                   ELSE
01284                      RESULT( 17 ) = ULPINV
01285                      GO TO 280
01286                   END IF
01287                END IF
01288 *
01289 *              Do test 17
01290 *
01291                TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) /
01292      $                 ( ONE-HALF )**4
01293 *
01294                TEMP1 = ZERO
01295                DO 190 J = 1, N
01296                   TEMP1 = MAX( TEMP1, ABS( D4( J )-WR( N-J+1 ) ) /
01297      $                    ( ABSTOL+ABS( D4( J ) ) ) )
01298   190          CONTINUE
01299 *
01300                RESULT( 17 ) = TEMP1 / TEMP2
01301             ELSE
01302                RESULT( 17 ) = ZERO
01303             END IF
01304 *
01305 *           Now ask for all eigenvalues with high absolute accuracy.
01306 *
01307             NTEST = 18
01308             ABSTOL = UNFL + UNFL
01309             CALL DSTEBZ( 'A', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE, M,
01310      $                   NSPLIT, WA1, IWORK( 1 ), IWORK( N+1 ), WORK,
01311      $                   IWORK( 2*N+1 ), IINFO )
01312             IF( IINFO.NE.0 ) THEN
01313                WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(A)', IINFO, N, JTYPE,
01314      $            IOLDSD
01315                INFO = ABS( IINFO )
01316                IF( IINFO.LT.0 ) THEN
01317                   RETURN
01318                ELSE
01319                   RESULT( 18 ) = ULPINV
01320                   GO TO 280
01321                END IF
01322             END IF
01323 *
01324 *           Do test 18
01325 *
01326             TEMP1 = ZERO
01327             TEMP2 = ZERO
01328             DO 200 J = 1, N
01329                TEMP1 = MAX( TEMP1, ABS( D3( J ) ), ABS( WA1( J ) ) )
01330                TEMP2 = MAX( TEMP2, ABS( D3( J )-WA1( J ) ) )
01331   200       CONTINUE
01332 *
01333             RESULT( 18 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
01334 *
01335 *           Choose random values for IL and IU, and ask for the
01336 *           IL-th through IU-th eigenvalues.
01337 *
01338             NTEST = 19
01339             IF( N.LE.1 ) THEN
01340                IL = 1
01341                IU = N
01342             ELSE
01343                IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
01344                IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
01345                IF( IU.LT.IL ) THEN
01346                   ITEMP = IU
01347                   IU = IL
01348                   IL = ITEMP
01349                END IF
01350             END IF
01351 *
01352             CALL DSTEBZ( 'I', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE,
01353      $                   M2, NSPLIT, WA2, IWORK( 1 ), IWORK( N+1 ),
01354      $                   WORK, IWORK( 2*N+1 ), IINFO )
01355             IF( IINFO.NE.0 ) THEN
01356                WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(I)', IINFO, N, JTYPE,
01357      $            IOLDSD
01358                INFO = ABS( IINFO )
01359                IF( IINFO.LT.0 ) THEN
01360                   RETURN
01361                ELSE
01362                   RESULT( 19 ) = ULPINV
01363                   GO TO 280
01364                END IF
01365             END IF
01366 *
01367 *           Determine the values VL and VU of the IL-th and IU-th
01368 *           eigenvalues and ask for all eigenvalues in this range.
01369 *
01370             IF( N.GT.0 ) THEN
01371                IF( IL.NE.1 ) THEN
01372                   VL = WA1( IL ) - MAX( HALF*( WA1( IL )-WA1( IL-1 ) ),
01373      $                 ULP*ANORM, TWO*RTUNFL )
01374                ELSE
01375                   VL = WA1( 1 ) - MAX( HALF*( WA1( N )-WA1( 1 ) ),
01376      $                 ULP*ANORM, TWO*RTUNFL )
01377                END IF
01378                IF( IU.NE.N ) THEN
01379                   VU = WA1( IU ) + MAX( HALF*( WA1( IU+1 )-WA1( IU ) ),
01380      $                 ULP*ANORM, TWO*RTUNFL )
01381                ELSE
01382                   VU = WA1( N ) + MAX( HALF*( WA1( N )-WA1( 1 ) ),
01383      $                 ULP*ANORM, TWO*RTUNFL )
01384                END IF
01385             ELSE
01386                VL = ZERO
01387                VU = ONE
01388             END IF
01389 *
01390             CALL DSTEBZ( 'V', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE,
01391      $                   M3, NSPLIT, WA3, IWORK( 1 ), IWORK( N+1 ),
01392      $                   WORK, IWORK( 2*N+1 ), IINFO )
01393             IF( IINFO.NE.0 ) THEN
01394                WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(V)', IINFO, N, JTYPE,
01395      $            IOLDSD
01396                INFO = ABS( IINFO )
01397                IF( IINFO.LT.0 ) THEN
01398                   RETURN
01399                ELSE
01400                   RESULT( 19 ) = ULPINV
01401                   GO TO 280
01402                END IF
01403             END IF
01404 *
01405             IF( M3.EQ.0 .AND. N.NE.0 ) THEN
01406                RESULT( 19 ) = ULPINV
01407                GO TO 280
01408             END IF
01409 *
01410 *           Do test 19
01411 *
01412             TEMP1 = DSXT1( 1, WA2, M2, WA3, M3, ABSTOL, ULP, UNFL )
01413             TEMP2 = DSXT1( 1, WA3, M3, WA2, M2, ABSTOL, ULP, UNFL )
01414             IF( N.GT.0 ) THEN
01415                TEMP3 = MAX( ABS( WA1( N ) ), ABS( WA1( 1 ) ) )
01416             ELSE
01417                TEMP3 = ZERO
01418             END IF
01419 *
01420             RESULT( 19 ) = ( TEMP1+TEMP2 ) / MAX( UNFL, TEMP3*ULP )
01421 *
01422 *           Call DSTEIN to compute eigenvectors corresponding to
01423 *           eigenvalues in WA1.  (First call DSTEBZ again, to make sure
01424 *           it returns these eigenvalues in the correct order.)
01425 *
01426             NTEST = 21
01427             CALL DSTEBZ( 'A', 'B', N, VL, VU, IL, IU, ABSTOL, SD, SE, M,
01428      $                   NSPLIT, WA1, IWORK( 1 ), IWORK( N+1 ), WORK,
01429      $                   IWORK( 2*N+1 ), IINFO )
01430             IF( IINFO.NE.0 ) THEN
01431                WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(A,B)', IINFO, N,
01432      $            JTYPE, IOLDSD
01433                INFO = ABS( IINFO )
01434                IF( IINFO.LT.0 ) THEN
01435                   RETURN
01436                ELSE
01437                   RESULT( 20 ) = ULPINV
01438                   RESULT( 21 ) = ULPINV
01439                   GO TO 280
01440                END IF
01441             END IF
01442 *
01443             CALL DSTEIN( N, SD, SE, M, WA1, IWORK( 1 ), IWORK( N+1 ), Z,
01444      $                   LDU, WORK, IWORK( 2*N+1 ), IWORK( 3*N+1 ),
01445      $                   IINFO )
01446             IF( IINFO.NE.0 ) THEN
01447                WRITE( NOUNIT, FMT = 9999 )'DSTEIN', IINFO, N, JTYPE,
01448      $            IOLDSD
01449                INFO = ABS( IINFO )
01450                IF( IINFO.LT.0 ) THEN
01451                   RETURN
01452                ELSE
01453                   RESULT( 20 ) = ULPINV
01454                   RESULT( 21 ) = ULPINV
01455                   GO TO 280
01456                END IF
01457             END IF
01458 *
01459 *           Do tests 20 and 21
01460 *
01461             CALL DSTT21( N, 0, SD, SE, WA1, DUMMA, Z, LDU, WORK,
01462      $                   RESULT( 20 ) )
01463 *
01464 *           Call DSTEDC(I) to compute D1 and Z, do tests.
01465 *
01466 *           Compute D1 and Z
01467 *
01468             CALL DCOPY( N, SD, 1, D1, 1 )
01469             IF( N.GT.0 )
01470      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
01471             CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
01472 *
01473             NTEST = 22
01474             CALL DSTEDC( 'I', N, D1, WORK, Z, LDU, WORK( N+1 ), LWEDC-N,
01475      $                   IWORK, LIWEDC, IINFO )
01476             IF( IINFO.NE.0 ) THEN
01477                WRITE( NOUNIT, FMT = 9999 )'DSTEDC(I)', IINFO, N, JTYPE,
01478      $            IOLDSD
01479                INFO = ABS( IINFO )
01480                IF( IINFO.LT.0 ) THEN
01481                   RETURN
01482                ELSE
01483                   RESULT( 22 ) = ULPINV
01484                   GO TO 280
01485                END IF
01486             END IF
01487 *
01488 *           Do Tests 22 and 23
01489 *
01490             CALL DSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
01491      $                   RESULT( 22 ) )
01492 *
01493 *           Call DSTEDC(V) to compute D1 and Z, do tests.
01494 *
01495 *           Compute D1 and Z
01496 *
01497             CALL DCOPY( N, SD, 1, D1, 1 )
01498             IF( N.GT.0 )
01499      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
01500             CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
01501 *
01502             NTEST = 24
01503             CALL DSTEDC( 'V', N, D1, WORK, Z, LDU, WORK( N+1 ), LWEDC-N,
01504      $                   IWORK, LIWEDC, IINFO )
01505             IF( IINFO.NE.0 ) THEN
01506                WRITE( NOUNIT, FMT = 9999 )'DSTEDC(V)', IINFO, N, JTYPE,
01507      $            IOLDSD
01508                INFO = ABS( IINFO )
01509                IF( IINFO.LT.0 ) THEN
01510                   RETURN
01511                ELSE
01512                   RESULT( 24 ) = ULPINV
01513                   GO TO 280
01514                END IF
01515             END IF
01516 *
01517 *           Do Tests 24 and 25
01518 *
01519             CALL DSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
01520      $                   RESULT( 24 ) )
01521 *
01522 *           Call DSTEDC(N) to compute D2, do tests.
01523 *
01524 *           Compute D2
01525 *
01526             CALL DCOPY( N, SD, 1, D2, 1 )
01527             IF( N.GT.0 )
01528      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
01529             CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
01530 *
01531             NTEST = 26
01532             CALL DSTEDC( 'N', N, D2, WORK, Z, LDU, WORK( N+1 ), LWEDC-N,
01533      $                   IWORK, LIWEDC, IINFO )
01534             IF( IINFO.NE.0 ) THEN
01535                WRITE( NOUNIT, FMT = 9999 )'DSTEDC(N)', IINFO, N, JTYPE,
01536      $            IOLDSD
01537                INFO = ABS( IINFO )
01538                IF( IINFO.LT.0 ) THEN
01539                   RETURN
01540                ELSE
01541                   RESULT( 26 ) = ULPINV
01542                   GO TO 280
01543                END IF
01544             END IF
01545 *
01546 *           Do Test 26
01547 *
01548             TEMP1 = ZERO
01549             TEMP2 = ZERO
01550 *
01551             DO 210 J = 1, N
01552                TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
01553                TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
01554   210       CONTINUE
01555 *
01556             RESULT( 26 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
01557 *
01558 *           Only test DSTEMR if IEEE compliant
01559 *
01560             IF( ILAENV( 10, 'DSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND.
01561      $          ILAENV( 11, 'DSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN
01562 *
01563 *           Call DSTEMR, do test 27 (relative eigenvalue accuracy)
01564 *
01565 *              If S is positive definite and diagonally dominant,
01566 *              ask for all eigenvalues with high relative accuracy.
01567 *
01568                VL = ZERO
01569                VU = ZERO
01570                IL = 0
01571                IU = 0
01572                IF( JTYPE.EQ.21 .AND. SREL ) THEN
01573                   NTEST = 27
01574                   ABSTOL = UNFL + UNFL
01575                   CALL DSTEMR( 'V', 'A', N, SD, SE, VL, VU, IL, IU,
01576      $                         M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC,
01577      $                         WORK, LWORK, IWORK( 2*N+1 ), LWORK-2*N,
01578      $                         IINFO )
01579                   IF( IINFO.NE.0 ) THEN
01580                      WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,A,rel)',
01581      $                  IINFO, N, JTYPE, IOLDSD
01582                      INFO = ABS( IINFO )
01583                      IF( IINFO.LT.0 ) THEN
01584                         RETURN
01585                      ELSE
01586                         RESULT( 27 ) = ULPINV
01587                         GO TO 270
01588                      END IF
01589                   END IF
01590 *
01591 *              Do test 27
01592 *
01593                   TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) /
01594      $                    ( ONE-HALF )**4
01595 *
01596                   TEMP1 = ZERO
01597                   DO 220 J = 1, N
01598                      TEMP1 = MAX( TEMP1, ABS( D4( J )-WR( N-J+1 ) ) /
01599      $                       ( ABSTOL+ABS( D4( J ) ) ) )
01600   220             CONTINUE
01601 *
01602                   RESULT( 27 ) = TEMP1 / TEMP2
01603 *
01604                   IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
01605                   IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
01606                   IF( IU.LT.IL ) THEN
01607                      ITEMP = IU
01608                      IU = IL
01609                      IL = ITEMP
01610                   END IF
01611 *
01612                   IF( SRANGE ) THEN
01613                      NTEST = 28
01614                      ABSTOL = UNFL + UNFL
01615                      CALL DSTEMR( 'V', 'I', N, SD, SE, VL, VU, IL, IU,
01616      $                            M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC,
01617      $                            WORK, LWORK, IWORK( 2*N+1 ),
01618      $                            LWORK-2*N, IINFO )
01619 *
01620                      IF( IINFO.NE.0 ) THEN
01621                         WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,I,rel)',
01622      $                     IINFO, N, JTYPE, IOLDSD
01623                         INFO = ABS( IINFO )
01624                         IF( IINFO.LT.0 ) THEN
01625                            RETURN
01626                         ELSE
01627                            RESULT( 28 ) = ULPINV
01628                            GO TO 270
01629                         END IF
01630                      END IF
01631 *
01632 *
01633 *                 Do test 28
01634 *
01635                      TEMP2 = TWO*( TWO*N-ONE )*ULP*
01636      $                       ( ONE+EIGHT*HALF**2 ) / ( ONE-HALF )**4
01637 *
01638                      TEMP1 = ZERO
01639                      DO 230 J = IL, IU
01640                         TEMP1 = MAX( TEMP1, ABS( WR( J-IL+1 )-D4( N-J+
01641      $                          1 ) ) / ( ABSTOL+ABS( WR( J-IL+1 ) ) ) )
01642   230                CONTINUE
01643 *
01644                      RESULT( 28 ) = TEMP1 / TEMP2
01645                   ELSE
01646                      RESULT( 28 ) = ZERO
01647                   END IF
01648                ELSE
01649                   RESULT( 27 ) = ZERO
01650                   RESULT( 28 ) = ZERO
01651                END IF
01652 *
01653 *           Call DSTEMR(V,I) to compute D1 and Z, do tests.
01654 *
01655 *           Compute D1 and Z
01656 *
01657                CALL DCOPY( N, SD, 1, D5, 1 )
01658                IF( N.GT.0 )
01659      $            CALL DCOPY( N-1, SE, 1, WORK, 1 )
01660                CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
01661 *
01662                IF( SRANGE ) THEN
01663                   NTEST = 29
01664                   IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
01665                   IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
01666                   IF( IU.LT.IL ) THEN
01667                      ITEMP = IU
01668                      IU = IL
01669                      IL = ITEMP
01670                   END IF
01671                   CALL DSTEMR( 'V', 'I', N, D5, WORK, VL, VU, IL, IU,
01672      $                         M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC,
01673      $                         WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
01674      $                         LIWORK-2*N, IINFO )
01675                   IF( IINFO.NE.0 ) THEN
01676                      WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,I)', IINFO,
01677      $                  N, JTYPE, IOLDSD
01678                      INFO = ABS( IINFO )
01679                      IF( IINFO.LT.0 ) THEN
01680                         RETURN
01681                      ELSE
01682                         RESULT( 29 ) = ULPINV
01683                         GO TO 280
01684                      END IF
01685                   END IF
01686 *
01687 *           Do Tests 29 and 30
01688 *
01689                   CALL DSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
01690      $                         M, RESULT( 29 ) )
01691 *
01692 *           Call DSTEMR to compute D2, do tests.
01693 *
01694 *           Compute D2
01695 *
01696                   CALL DCOPY( N, SD, 1, D5, 1 )
01697                   IF( N.GT.0 )
01698      $               CALL DCOPY( N-1, SE, 1, WORK, 1 )
01699 *
01700                   NTEST = 31
01701                   CALL DSTEMR( 'N', 'I', N, D5, WORK, VL, VU, IL, IU,
01702      $                         M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC,
01703      $                         WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
01704      $                         LIWORK-2*N, IINFO )
01705                   IF( IINFO.NE.0 ) THEN
01706                      WRITE( NOUNIT, FMT = 9999 )'DSTEMR(N,I)', IINFO,
01707      $                  N, JTYPE, IOLDSD
01708                      INFO = ABS( IINFO )
01709                      IF( IINFO.LT.0 ) THEN
01710                         RETURN
01711                      ELSE
01712                         RESULT( 31 ) = ULPINV
01713                         GO TO 280
01714                      END IF
01715                   END IF
01716 *
01717 *           Do Test 31
01718 *
01719                   TEMP1 = ZERO
01720                   TEMP2 = ZERO
01721 *
01722                   DO 240 J = 1, IU - IL + 1
01723                      TEMP1 = MAX( TEMP1, ABS( D1( J ) ),
01724      $                       ABS( D2( J ) ) )
01725                      TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
01726   240             CONTINUE
01727 *
01728                   RESULT( 31 ) = TEMP2 / MAX( UNFL,
01729      $                           ULP*MAX( TEMP1, TEMP2 ) )
01730 *
01731 *
01732 *           Call DSTEMR(V,V) to compute D1 and Z, do tests.
01733 *
01734 *           Compute D1 and Z
01735 *
01736                   CALL DCOPY( N, SD, 1, D5, 1 )
01737                   IF( N.GT.0 )
01738      $               CALL DCOPY( N-1, SE, 1, WORK, 1 )
01739                   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
01740 *
01741                   NTEST = 32
01742 *
01743                   IF( N.GT.0 ) THEN
01744                      IF( IL.NE.1 ) THEN
01745                         VL = D2( IL ) - MAX( HALF*
01746      $                       ( D2( IL )-D2( IL-1 ) ), ULP*ANORM,
01747      $                       TWO*RTUNFL )
01748                      ELSE
01749                         VL = D2( 1 ) - MAX( HALF*( D2( N )-D2( 1 ) ),
01750      $                       ULP*ANORM, TWO*RTUNFL )
01751                      END IF
01752                      IF( IU.NE.N ) THEN
01753                         VU = D2( IU ) + MAX( HALF*
01754      $                       ( D2( IU+1 )-D2( IU ) ), ULP*ANORM,
01755      $                       TWO*RTUNFL )
01756                      ELSE
01757                         VU = D2( N ) + MAX( HALF*( D2( N )-D2( 1 ) ),
01758      $                       ULP*ANORM, TWO*RTUNFL )
01759                      END IF
01760                   ELSE
01761                      VL = ZERO
01762                      VU = ONE
01763                   END IF
01764 *
01765                   CALL DSTEMR( 'V', 'V', N, D5, WORK, VL, VU, IL, IU,
01766      $                         M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC,
01767      $                         WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
01768      $                         LIWORK-2*N, IINFO )
01769                   IF( IINFO.NE.0 ) THEN
01770                      WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,V)', IINFO,
01771      $                  N, JTYPE, IOLDSD
01772                      INFO = ABS( IINFO )
01773                      IF( IINFO.LT.0 ) THEN
01774                         RETURN
01775                      ELSE
01776                         RESULT( 32 ) = ULPINV
01777                         GO TO 280
01778                      END IF
01779                   END IF
01780 *
01781 *           Do Tests 32 and 33
01782 *
01783                   CALL DSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
01784      $                         M, RESULT( 32 ) )
01785 *
01786 *           Call DSTEMR to compute D2, do tests.
01787 *
01788 *           Compute D2
01789 *
01790                   CALL DCOPY( N, SD, 1, D5, 1 )
01791                   IF( N.GT.0 )
01792      $               CALL DCOPY( N-1, SE, 1, WORK, 1 )
01793 *
01794                   NTEST = 34
01795                   CALL DSTEMR( 'N', 'V', N, D5, WORK, VL, VU, IL, IU,
01796      $                         M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC,
01797      $                         WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
01798      $                         LIWORK-2*N, IINFO )
01799                   IF( IINFO.NE.0 ) THEN
01800                      WRITE( NOUNIT, FMT = 9999 )'DSTEMR(N,V)', IINFO,
01801      $                  N, JTYPE, IOLDSD
01802                      INFO = ABS( IINFO )
01803                      IF( IINFO.LT.0 ) THEN
01804                         RETURN
01805                      ELSE
01806                         RESULT( 34 ) = ULPINV
01807                         GO TO 280
01808                      END IF
01809                   END IF
01810 *
01811 *           Do Test 34
01812 *
01813                   TEMP1 = ZERO
01814                   TEMP2 = ZERO
01815 *
01816                   DO 250 J = 1, IU - IL + 1
01817                      TEMP1 = MAX( TEMP1, ABS( D1( J ) ),
01818      $                       ABS( D2( J ) ) )
01819                      TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
01820   250             CONTINUE
01821 *
01822                   RESULT( 34 ) = TEMP2 / MAX( UNFL,
01823      $                           ULP*MAX( TEMP1, TEMP2 ) )
01824                ELSE
01825                   RESULT( 29 ) = ZERO
01826                   RESULT( 30 ) = ZERO
01827                   RESULT( 31 ) = ZERO
01828                   RESULT( 32 ) = ZERO
01829                   RESULT( 33 ) = ZERO
01830                   RESULT( 34 ) = ZERO
01831                END IF
01832 *
01833 *
01834 *           Call DSTEMR(V,A) to compute D1 and Z, do tests.
01835 *
01836 *           Compute D1 and Z
01837 *
01838                CALL DCOPY( N, SD, 1, D5, 1 )
01839                IF( N.GT.0 )
01840      $            CALL DCOPY( N-1, SE, 1, WORK, 1 )
01841 *
01842                NTEST = 35
01843 *
01844                CALL DSTEMR( 'V', 'A', N, D5, WORK, VL, VU, IL, IU,
01845      $                      M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC,
01846      $                      WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
01847      $                      LIWORK-2*N, IINFO )
01848                IF( IINFO.NE.0 ) THEN
01849                   WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,A)', IINFO, N,
01850      $               JTYPE, IOLDSD
01851                   INFO = ABS( IINFO )
01852                   IF( IINFO.LT.0 ) THEN
01853                      RETURN
01854                   ELSE
01855                      RESULT( 35 ) = ULPINV
01856                      GO TO 280
01857                   END IF
01858                END IF
01859 *
01860 *           Do Tests 35 and 36
01861 *
01862                CALL DSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, M,
01863      $                      RESULT( 35 ) )
01864 *
01865 *           Call DSTEMR to compute D2, do tests.
01866 *
01867 *           Compute D2
01868 *
01869                CALL DCOPY( N, SD, 1, D5, 1 )
01870                IF( N.GT.0 )
01871      $            CALL DCOPY( N-1, SE, 1, WORK, 1 )
01872 *
01873                NTEST = 37
01874                CALL DSTEMR( 'N', 'A', N, D5, WORK, VL, VU, IL, IU,
01875      $                      M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC,
01876      $                      WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
01877      $                      LIWORK-2*N, IINFO )
01878                IF( IINFO.NE.0 ) THEN
01879                   WRITE( NOUNIT, FMT = 9999 )'DSTEMR(N,A)', IINFO, N,
01880      $               JTYPE, IOLDSD
01881                   INFO = ABS( IINFO )
01882                   IF( IINFO.LT.0 ) THEN
01883                      RETURN
01884                   ELSE
01885                      RESULT( 37 ) = ULPINV
01886                      GO TO 280
01887                   END IF
01888                END IF
01889 *
01890 *           Do Test 34
01891 *
01892                TEMP1 = ZERO
01893                TEMP2 = ZERO
01894 *
01895                DO 260 J = 1, N
01896                   TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
01897                   TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
01898   260          CONTINUE
01899 *
01900                RESULT( 37 ) = TEMP2 / MAX( UNFL,
01901      $                        ULP*MAX( TEMP1, TEMP2 ) )
01902             END IF
01903   270       CONTINUE
01904   280       CONTINUE
01905             NTESTT = NTESTT + NTEST
01906 *
01907 *           End of Loop -- Check for RESULT(j) > THRESH
01908 *
01909 *
01910 *           Print out tests which fail.
01911 *
01912             DO 290 JR = 1, NTEST
01913                IF( RESULT( JR ).GE.THRESH ) THEN
01914 *
01915 *                 If this is the first test to fail,
01916 *                 print a header to the data file.
01917 *
01918                   IF( NERRS.EQ.0 ) THEN
01919                      WRITE( NOUNIT, FMT = 9998 )'DST'
01920                      WRITE( NOUNIT, FMT = 9997 )
01921                      WRITE( NOUNIT, FMT = 9996 )
01922                      WRITE( NOUNIT, FMT = 9995 )'Symmetric'
01923                      WRITE( NOUNIT, FMT = 9994 )
01924 *
01925 *                    Tests performed
01926 *
01927                      WRITE( NOUNIT, FMT = 9988 )
01928                   END IF
01929                   NERRS = NERRS + 1
01930                   WRITE( NOUNIT, FMT = 9990 )N, IOLDSD, JTYPE, JR,
01931      $               RESULT( JR )
01932                END IF
01933   290       CONTINUE
01934   300    CONTINUE
01935   310 CONTINUE
01936 *
01937 *     Summary
01938 *
01939       CALL DLASUM( 'DST', NOUNIT, NERRS, NTESTT )
01940       RETURN
01941 *
01942  9999 FORMAT( ' DCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
01943      $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
01944 *
01945  9998 FORMAT( / 1X, A3, ' -- Real Symmetric eigenvalue problem' )
01946  9997 FORMAT( ' Matrix types (see DCHKST for details): ' )
01947 *
01948  9996 FORMAT( / ' Special Matrices:',
01949      $      / '  1=Zero matrix.                        ',
01950      $      '  5=Diagonal: clustered entries.',
01951      $      / '  2=Identity matrix.                    ',
01952      $      '  6=Diagonal: large, evenly spaced.',
01953      $      / '  3=Diagonal: evenly spaced entries.    ',
01954      $      '  7=Diagonal: small, evenly spaced.',
01955      $      / '  4=Diagonal: geometr. spaced entries.' )
01956  9995 FORMAT( ' Dense ', A, ' Matrices:',
01957      $      / '  8=Evenly spaced eigenvals.            ',
01958      $      ' 12=Small, evenly spaced eigenvals.',
01959      $      / '  9=Geometrically spaced eigenvals.     ',
01960      $      ' 13=Matrix with random O(1) entries.',
01961      $      / ' 10=Clustered eigenvalues.              ',
01962      $      ' 14=Matrix with large random entries.',
01963      $      / ' 11=Large, evenly spaced eigenvals.     ',
01964      $      ' 15=Matrix with small random entries.' )
01965  9994 FORMAT( ' 16=Positive definite, evenly spaced eigenvalues',
01966      $      / ' 17=Positive definite, geometrically spaced eigenvlaues',
01967      $      / ' 18=Positive definite, clustered eigenvalues',
01968      $      / ' 19=Positive definite, small evenly spaced eigenvalues',
01969      $      / ' 20=Positive definite, large evenly spaced eigenvalues',
01970      $      / ' 21=Diagonally dominant tridiagonal, geometrically',
01971      $      ' spaced eigenvalues' )
01972 *
01973  9990 FORMAT( ' N=', I5, ', seed=', 4( I4, ',' ), ' type ', I2,
01974      $      ', test(', I2, ')=', G10.3 )
01975 *
01976  9988 FORMAT( / 'Test performed:  see DCHKST for details.', / )
01977 *     End of DCHKST
01978 *
01979       END
 All Files Functions