lapack-netlib/TESTING/EIG/dchkst2stg.f

   1 *> \brief \b DCHKST2STG
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE DCHKST2STG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
  12 *                          NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5,
  13 *                          WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK,
  14 *                          LWORK, IWORK, LIWORK, RESULT, INFO )
  15 *
  16 *       .. Scalar Arguments ..
  17 *       INTEGER            INFO, LDA, LDU, LIWORK, LWORK, NOUNIT, NSIZES,
  18 *      $                   NTYPES
  19 *       DOUBLE PRECISION   THRESH
  20 *       ..
  21 *       .. Array Arguments ..
  22 *       LOGICAL            DOTYPE( * )
  23 *       INTEGER            ISEED( 4 ), IWORK( * ), NN( * )
  24 *       DOUBLE PRECISION   A( LDA, * ), AP( * ), D1( * ), D2( * ),
  25 *      $                   D3( * ), D4( * ), D5( * ), RESULT( * ),
  26 *      $                   SD( * ), SE( * ), TAU( * ), U( LDU, * ),
  27 *      $                   V( LDU, * ), VP( * ), WA1( * ), WA2( * ),
  28 *      $                   WA3( * ), WORK( * ), WR( * ), Z( LDU, * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DCHKST2STG  checks the symmetric eigenvalue problem routines
  38 *> using the 2-stage reduction techniques. Since the generation
  39 *> of Q or the vectors is not available in this release, we only
  40 *> compare the eigenvalue resulting when using the 2-stage to the
  41 *> one considered as reference using the standard 1-stage reduction
  42 *> DSYTRD. For that, we call the standard DSYTRD and compute D1 using
  43 *> DSTEQR, then we call the 2-stage DSYTRD_2STAGE with Upper and Lower
  44 *> and we compute D2 and D3 using DSTEQR and then we replaced tests
  45 *> 3 and 4 by tests 11 and 12. test 1 and 2 remain to verify that
  46 *> the 1-stage results are OK and can be trusted.
  47 *> This testing routine will converge to the DCHKST in the next
  48 *> release when vectors and generation of Q will be implemented.
  49 *>
  50 *>    DSYTRD factors A as  U S U' , where ' means transpose,
  51 *>    S is symmetric tridiagonal, and U is orthogonal.
  52 *>    DSYTRD can use either just the lower or just the upper triangle
  53 *>    of A; DCHKST2STG checks both cases.
  54 *>    U is represented as a product of Householder
  55 *>    transformations, whose vectors are stored in the first
  56 *>    n-1 columns of V, and whose scale factors are in TAU.
  57 *>
  58 *>    DSPTRD does the same as DSYTRD, except that A and V are stored
  59 *>    in "packed" format.
  60 *>
  61 *>    DORGTR constructs the matrix U from the contents of V and TAU.
  62 *>
  63 *>    DOPGTR constructs the matrix U from the contents of VP and TAU.
  64 *>
  65 *>    DSTEQR factors S as  Z D1 Z' , where Z is the orthogonal
  66 *>    matrix of eigenvectors and D1 is a diagonal matrix with
  67 *>    the eigenvalues on the diagonal.  D2 is the matrix of
  68 *>    eigenvalues computed when Z is not computed.
  69 *>
  70 *>    DSTERF computes D3, the matrix of eigenvalues, by the
  71 *>    PWK method, which does not yield eigenvectors.
  72 *>
  73 *>    DPTEQR factors S as  Z4 D4 Z4' , for a
  74 *>    symmetric positive definite tridiagonal matrix.
  75 *>    D5 is the matrix of eigenvalues computed when Z is not
  76 *>    computed.
  77 *>
  78 *>    DSTEBZ computes selected eigenvalues.  WA1, WA2, and
  79 *>    WA3 will denote eigenvalues computed to high
  80 *>    absolute accuracy, with different range options.
  81 *>    WR will denote eigenvalues computed to high relative
  82 *>    accuracy.
  83 *>
  84 *>    DSTEIN computes Y, the eigenvectors of S, given the
  85 *>    eigenvalues.
  86 *>
  87 *>    DSTEDC factors S as Z D1 Z' , where Z is the orthogonal
  88 *>    matrix of eigenvectors and D1 is a diagonal matrix with
  89 *>    the eigenvalues on the diagonal ('I' option). It may also
  90 *>    update an input orthogonal matrix, usually the output
  91 *>    from DSYTRD/DORGTR or DSPTRD/DOPGTR ('V' option). It may
  92 *>    also just compute eigenvalues ('N' option).
  93 *>
  94 *>    DSTEMR factors S as Z D1 Z' , where Z is the orthogonal
  95 *>    matrix of eigenvectors and D1 is a diagonal matrix with
  96 *>    the eigenvalues on the diagonal ('I' option).  DSTEMR
  97 *>    uses the Relatively Robust Representation whenever possible.
  98 *>
  99 *> When DCHKST2STG is called, a number of matrix "sizes" ("n's") and a
 100 *> number of matrix "types" are specified.  For each size ("n")
 101 *> and each type of matrix, one matrix will be generated and used
 102 *> to test the symmetric eigenroutines.  For each matrix, a number
 103 *> of tests will be performed:
 104 *>
 105 *> (1)     | A - V S V' | / ( |A| n ulp ) DSYTRD( UPLO='U', ... )
 106 *>
 107 *> (2)     | I - UV' | / ( n ulp )        DORGTR( UPLO='U', ... )
 108 *>
 109 *> (3)     | A - V S V' | / ( |A| n ulp ) DSYTRD( UPLO='L', ... )
 110 *>         replaced by | D1 - D2 | / ( |D1| ulp ) where D1 is the
 111 *>         eigenvalue matrix computed using S and D2 is the
 112 *>         eigenvalue matrix computed using S_2stage the output of
 113 *>         DSYTRD_2STAGE("N", "U",....). D1 and D2 are computed
 114 *>         via DSTEQR('N',...)
 115 *>
 116 *> (4)     | I - UV' | / ( n ulp )        DORGTR( UPLO='L', ... )
 117 *>         replaced by | D1 - D3 | / ( |D1| ulp ) where D1 is the
 118 *>         eigenvalue matrix computed using S and D3 is the
 119 *>         eigenvalue matrix computed using S_2stage the output of
 120 *>         DSYTRD_2STAGE("N", "L",....). D1 and D3 are computed
 121 *>         via DSTEQR('N',...)
 122 *>
 123 *> (5-8)   Same as 1-4, but for DSPTRD and DOPGTR.
 124 *>
 125 *> (9)     | S - Z D Z' | / ( |S| n ulp ) DSTEQR('V',...)
 126 *>
 127 *> (10)    | I - ZZ' | / ( n ulp )        DSTEQR('V',...)
 128 *>
 129 *> (11)    | D1 - D2 | / ( |D1| ulp )        DSTEQR('N',...)
 130 *>
 131 *> (12)    | D1 - D3 | / ( |D1| ulp )        DSTERF
 132 *>
 133 *> (13)    0 if the true eigenvalues (computed by sturm count)
 134 *>         of S are within THRESH of
 135 *>         those in D1.  2*THRESH if they are not.  (Tested using
 136 *>         DSTECH)
 137 *>
 138 *> For S positive definite,
 139 *>
 140 *> (14)    | S - Z4 D4 Z4' | / ( |S| n ulp ) DPTEQR('V',...)
 141 *>
 142 *> (15)    | I - Z4 Z4' | / ( n ulp )        DPTEQR('V',...)
 143 *>
 144 *> (16)    | D4 - D5 | / ( 100 |D4| ulp )       DPTEQR('N',...)
 145 *>
 146 *> When S is also diagonally dominant by the factor gamma < 1,
 147 *>
 148 *> (17)    max | D4(i) - WR(i) | / ( |D4(i)| omega ) ,
 149 *>          i
 150 *>         omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
 151 *>                                              DSTEBZ( 'A', 'E', ...)
 152 *>
 153 *> (18)    | WA1 - D3 | / ( |D3| ulp )          DSTEBZ( 'A', 'E', ...)
 154 *>
 155 *> (19)    ( max { min | WA2(i)-WA3(j) | } +
 156 *>            i     j
 157 *>           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
 158 *>            i     j
 159 *>                                              DSTEBZ( 'I', 'E', ...)
 160 *>
 161 *> (20)    | S - Y WA1 Y' | / ( |S| n ulp )  DSTEBZ, SSTEIN
 162 *>
 163 *> (21)    | I - Y Y' | / ( n ulp )          DSTEBZ, SSTEIN
 164 *>
 165 *> (22)    | S - Z D Z' | / ( |S| n ulp )    DSTEDC('I')
 166 *>
 167 *> (23)    | I - ZZ' | / ( n ulp )           DSTEDC('I')
 168 *>
 169 *> (24)    | S - Z D Z' | / ( |S| n ulp )    DSTEDC('V')
 170 *>
 171 *> (25)    | I - ZZ' | / ( n ulp )           DSTEDC('V')
 172 *>
 173 *> (26)    | D1 - D2 | / ( |D1| ulp )           DSTEDC('V') and
 174 *>                                              DSTEDC('N')
 175 *>
 176 *> Test 27 is disabled at the moment because DSTEMR does not
 177 *> guarantee high relatvie accuracy.
 178 *>
 179 *> (27)    max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
 180 *>          i
 181 *>         omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
 182 *>                                              DSTEMR('V', 'A')
 183 *>
 184 *> (28)    max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
 185 *>          i
 186 *>         omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
 187 *>                                              DSTEMR('V', 'I')
 188 *>
 189 *> Tests 29 through 34 are disable at present because DSTEMR
 190 *> does not handle partial spectrum requests.
 191 *>
 192 *> (29)    | S - Z D Z' | / ( |S| n ulp )    DSTEMR('V', 'I')
 193 *>
 194 *> (30)    | I - ZZ' | / ( n ulp )           DSTEMR('V', 'I')
 195 *>
 196 *> (31)    ( max { min | WA2(i)-WA3(j) | } +
 197 *>            i     j
 198 *>           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
 199 *>            i     j
 200 *>         DSTEMR('N', 'I') vs. SSTEMR('V', 'I')
 201 *>
 202 *> (32)    | S - Z D Z' | / ( |S| n ulp )    DSTEMR('V', 'V')
 203 *>
 204 *> (33)    | I - ZZ' | / ( n ulp )           DSTEMR('V', 'V')
 205 *>
 206 *> (34)    ( max { min | WA2(i)-WA3(j) | } +
 207 *>            i     j
 208 *>           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
 209 *>            i     j
 210 *>         DSTEMR('N', 'V') vs. SSTEMR('V', 'V')
 211 *>
 212 *> (35)    | S - Z D Z' | / ( |S| n ulp )    DSTEMR('V', 'A')
 213 *>
 214 *> (36)    | I - ZZ' | / ( n ulp )           DSTEMR('V', 'A')
 215 *>
 216 *> (37)    ( max { min | WA2(i)-WA3(j) | } +
 217 *>            i     j
 218 *>           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
 219 *>            i     j
 220 *>         DSTEMR('N', 'A') vs. SSTEMR('V', 'A')
 221 *>
 222 *> The "sizes" are specified by an array NN(1:NSIZES); the value of
 223 *> each element NN(j) specifies one size.
 224 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
 225 *> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
 226 *> Currently, the list of possible types is:
 227 *>
 228 *> (1)  The zero matrix.
 229 *> (2)  The identity matrix.
 230 *>
 231 *> (3)  A diagonal matrix with evenly spaced entries
 232 *>      1, ..., ULP  and random signs.
 233 *>      (ULP = (first number larger than 1) - 1 )
 234 *> (4)  A diagonal matrix with geometrically spaced entries
 235 *>      1, ..., ULP  and random signs.
 236 *> (5)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
 237 *>      and random signs.
 238 *>
 239 *> (6)  Same as (4), but multiplied by SQRT( overflow threshold )
 240 *> (7)  Same as (4), but multiplied by SQRT( underflow threshold )
 241 *>
 242 *> (8)  A matrix of the form  U' D U, where U is orthogonal and
 243 *>      D has evenly spaced entries 1, ..., ULP with random signs
 244 *>      on the diagonal.
 245 *>
 246 *> (9)  A matrix of the form  U' D U, where U is orthogonal and
 247 *>      D has geometrically spaced entries 1, ..., ULP with random
 248 *>      signs on the diagonal.
 249 *>
 250 *> (10) A matrix of the form  U' D U, where U is orthogonal and
 251 *>      D has "clustered" entries 1, ULP,..., ULP with random
 252 *>      signs on the diagonal.
 253 *>
 254 *> (11) Same as (8), but multiplied by SQRT( overflow threshold )
 255 *> (12) Same as (8), but multiplied by SQRT( underflow threshold )
 256 *>
 257 *> (13) Symmetric matrix with random entries chosen from (-1,1).
 258 *> (14) Same as (13), but multiplied by SQRT( overflow threshold )
 259 *> (15) Same as (13), but multiplied by SQRT( underflow threshold )
 260 *> (16) Same as (8), but diagonal elements are all positive.
 261 *> (17) Same as (9), but diagonal elements are all positive.
 262 *> (18) Same as (10), but diagonal elements are all positive.
 263 *> (19) Same as (16), but multiplied by SQRT( overflow threshold )
 264 *> (20) Same as (16), but multiplied by SQRT( underflow threshold )
 265 *> (21) A diagonally dominant tridiagonal matrix with geometrically
 266 *>      spaced diagonal entries 1, ..., ULP.
 267 *> \endverbatim
 268 *
 269 *  Arguments:
 270 *  ==========
 271 *
 272 *> \param[in] NSIZES
 273 *> \verbatim
 274 *>          NSIZES is INTEGER
 275 *>          The number of sizes of matrices to use.  If it is zero,
 276 *>          DCHKST2STG does nothing.  It must be at least zero.
 277 *> \endverbatim
 278 *>
 279 *> \param[in] NN
 280 *> \verbatim
 281 *>          NN is INTEGER array, dimension (NSIZES)
 282 *>          An array containing the sizes to be used for the matrices.
 283 *>          Zero values will be skipped.  The values must be at least
 284 *>          zero.
 285 *> \endverbatim
 286 *>
 287 *> \param[in] NTYPES
 288 *> \verbatim
 289 *>          NTYPES is INTEGER
 290 *>          The number of elements in DOTYPE.   If it is zero, DCHKST2STG
 291 *>          does nothing.  It must be at least zero.  If it is MAXTYP+1
 292 *>          and NSIZES is 1, then an additional type, MAXTYP+1 is
 293 *>          defined, which is to use whatever matrix is in A.  This
 294 *>          is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
 295 *>          DOTYPE(MAXTYP+1) is .TRUE. .
 296 *> \endverbatim
 297 *>
 298 *> \param[in] DOTYPE
 299 *> \verbatim
 300 *>          DOTYPE is LOGICAL array, dimension (NTYPES)
 301 *>          If DOTYPE(j) is .TRUE., then for each size in NN a
 302 *>          matrix of that size and of type j will be generated.
 303 *>          If NTYPES is smaller than the maximum number of types
 304 *>          defined (PARAMETER MAXTYP), then types NTYPES+1 through
 305 *>          MAXTYP will not be generated.  If NTYPES is larger
 306 *>          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
 307 *>          will be ignored.
 308 *> \endverbatim
 309 *>
 310 *> \param[in,out] ISEED
 311 *> \verbatim
 312 *>          ISEED is INTEGER array, dimension (4)
 313 *>          On entry ISEED specifies the seed of the random number
 314 *>          generator. The array elements should be between 0 and 4095;
 315 *>          if not they will be reduced mod 4096.  Also, ISEED(4) must
 316 *>          be odd.  The random number generator uses a linear
 317 *>          congruential sequence limited to small integers, and so
 318 *>          should produce machine independent random numbers. The
 319 *>          values of ISEED are changed on exit, and can be used in the
 320 *>          next call to DCHKST2STG to continue the same random number
 321 *>          sequence.
 322 *> \endverbatim
 323 *>
 324 *> \param[in] THRESH
 325 *> \verbatim
 326 *>          THRESH is DOUBLE PRECISION
 327 *>          A test will count as "failed" if the "error", computed as
 328 *>          described above, exceeds THRESH.  Note that the error
 329 *>          is scaled to be O(1), so THRESH should be a reasonably
 330 *>          small multiple of 1, e.g., 10 or 100.  In particular,
 331 *>          it should not depend on the precision (single vs. double)
 332 *>          or the size of the matrix.  It must be at least zero.
 333 *> \endverbatim
 334 *>
 335 *> \param[in] NOUNIT
 336 *> \verbatim
 337 *>          NOUNIT is INTEGER
 338 *>          The FORTRAN unit number for printing out error messages
 339 *>          (e.g., if a routine returns IINFO not equal to 0.)
 340 *> \endverbatim
 341 *>
 342 *> \param[in,out] A
 343 *> \verbatim
 344 *>          A is DOUBLE PRECISION array of
 345 *>                                  dimension ( LDA , max(NN) )
 346 *>          Used to hold the matrix whose eigenvalues are to be
 347 *>          computed.  On exit, A contains the last matrix actually
 348 *>          used.
 349 *> \endverbatim
 350 *>
 351 *> \param[in] LDA
 352 *> \verbatim
 353 *>          LDA is INTEGER
 354 *>          The leading dimension of A.  It must be at
 355 *>          least 1 and at least max( NN ).
 356 *> \endverbatim
 357 *>
 358 *> \param[out] AP
 359 *> \verbatim
 360 *>          AP is DOUBLE PRECISION array of
 361 *>                      dimension( max(NN)*max(NN+1)/2 )
 362 *>          The matrix A stored in packed format.
 363 *> \endverbatim
 364 *>
 365 *> \param[out] SD
 366 *> \verbatim
 367 *>          SD is DOUBLE PRECISION array of
 368 *>                             dimension( max(NN) )
 369 *>          The diagonal of the tridiagonal matrix computed by DSYTRD.
 370 *>          On exit, SD and SE contain the tridiagonal form of the
 371 *>          matrix in A.
 372 *> \endverbatim
 373 *>
 374 *> \param[out] SE
 375 *> \verbatim
 376 *>          SE is DOUBLE PRECISION array of
 377 *>                             dimension( max(NN) )
 378 *>          The off-diagonal of the tridiagonal matrix computed by
 379 *>          DSYTRD.  On exit, SD and SE contain the tridiagonal form of
 380 *>          the matrix in A.
 381 *> \endverbatim
 382 *>
 383 *> \param[out] D1
 384 *> \verbatim
 385 *>          D1 is DOUBLE PRECISION array of
 386 *>                             dimension( max(NN) )
 387 *>          The eigenvalues of A, as computed by DSTEQR simlutaneously
 388 *>          with Z.  On exit, the eigenvalues in D1 correspond with the
 389 *>          matrix in A.
 390 *> \endverbatim
 391 *>
 392 *> \param[out] D2
 393 *> \verbatim
 394 *>          D2 is DOUBLE PRECISION array of
 395 *>                             dimension( max(NN) )
 396 *>          The eigenvalues of A, as computed by DSTEQR if Z is not
 397 *>          computed.  On exit, the eigenvalues in D2 correspond with
 398 *>          the matrix in A.
 399 *> \endverbatim
 400 *>
 401 *> \param[out] D3
 402 *> \verbatim
 403 *>          D3 is DOUBLE PRECISION array of
 404 *>                             dimension( max(NN) )
 405 *>          The eigenvalues of A, as computed by DSTERF.  On exit, the
 406 *>          eigenvalues in D3 correspond with the matrix in A.
 407 *> \endverbatim
 408 *>
 409 *> \param[out] D4
 410 *> \verbatim
 411 *>          D4 is DOUBLE PRECISION array of
 412 *>                             dimension( max(NN) )
 413 *>          The eigenvalues of A, as computed by DPTEQR(V).
 414 *>          DPTEQR factors S as  Z4 D4 Z4*
 415 *>          On exit, the eigenvalues in D4 correspond with the matrix in A.
 416 *> \endverbatim
 417 *>
 418 *> \param[out] D5
 419 *> \verbatim
 420 *>          D5 is DOUBLE PRECISION array of
 421 *>                             dimension( max(NN) )
 422 *>          The eigenvalues of A, as computed by DPTEQR(N)
 423 *>          when Z is not computed. On exit, the
 424 *>          eigenvalues in D4 correspond with the matrix in A.
 425 *> \endverbatim
 426 *>
 427 *> \param[out] WA1
 428 *> \verbatim
 429 *>          WA1 is DOUBLE PRECISION array of
 430 *>                             dimension( max(NN) )
 431 *>          All eigenvalues of A, computed to high
 432 *>          absolute accuracy, with different range options.
 433 *>          as computed by DSTEBZ.
 434 *> \endverbatim
 435 *>
 436 *> \param[out] WA2
 437 *> \verbatim
 438 *>          WA2 is DOUBLE PRECISION array of
 439 *>                             dimension( max(NN) )
 440 *>          Selected eigenvalues of A, computed to high
 441 *>          absolute accuracy, with different range options.
 442 *>          as computed by DSTEBZ.
 443 *>          Choose random values for IL and IU, and ask for the
 444 *>          IL-th through IU-th eigenvalues.
 445 *> \endverbatim
 446 *>
 447 *> \param[out] WA3
 448 *> \verbatim
 449 *>          WA3 is DOUBLE PRECISION array of
 450 *>                             dimension( max(NN) )
 451 *>          Selected eigenvalues of A, computed to high
 452 *>          absolute accuracy, with different range options.
 453 *>          as computed by DSTEBZ.
 454 *>          Determine the values VL and VU of the IL-th and IU-th
 455 *>          eigenvalues and ask for all eigenvalues in this range.
 456 *> \endverbatim
 457 *>
 458 *> \param[out] WR
 459 *> \verbatim
 460 *>          WR is DOUBLE PRECISION array of
 461 *>                             dimension( max(NN) )
 462 *>          All eigenvalues of A, computed to high
 463 *>          absolute accuracy, with different options.
 464 *>          as computed by DSTEBZ.
 465 *> \endverbatim
 466 *>
 467 *> \param[out] U
 468 *> \verbatim
 469 *>          U is DOUBLE PRECISION array of
 470 *>                             dimension( LDU, max(NN) ).
 471 *>          The orthogonal matrix computed by DSYTRD + DORGTR.
 472 *> \endverbatim
 473 *>
 474 *> \param[in] LDU
 475 *> \verbatim
 476 *>          LDU is INTEGER
 477 *>          The leading dimension of U, Z, and V.  It must be at least 1
 478 *>          and at least max( NN ).
 479 *> \endverbatim
 480 *>
 481 *> \param[out] V
 482 *> \verbatim
 483 *>          V is DOUBLE PRECISION array of
 484 *>                             dimension( LDU, max(NN) ).
 485 *>          The Housholder vectors computed by DSYTRD in reducing A to
 486 *>          tridiagonal form.  The vectors computed with UPLO='U' are
 487 *>          in the upper triangle, and the vectors computed with UPLO='L'
 488 *>          are in the lower triangle.  (As described in DSYTRD, the
 489 *>          sub- and superdiagonal are not set to 1, although the
 490 *>          true Householder vector has a 1 in that position.  The
 491 *>          routines that use V, such as DORGTR, set those entries to
 492 *>          1 before using them, and then restore them later.)
 493 *> \endverbatim
 494 *>
 495 *> \param[out] VP
 496 *> \verbatim
 497 *>          VP is DOUBLE PRECISION array of
 498 *>                      dimension( max(NN)*max(NN+1)/2 )
 499 *>          The matrix V stored in packed format.
 500 *> \endverbatim
 501 *>
 502 *> \param[out] TAU
 503 *> \verbatim
 504 *>          TAU is DOUBLE PRECISION array of
 505 *>                             dimension( max(NN) )
 506 *>          The Householder factors computed by DSYTRD in reducing A
 507 *>          to tridiagonal form.
 508 *> \endverbatim
 509 *>
 510 *> \param[out] Z
 511 *> \verbatim
 512 *>          Z is DOUBLE PRECISION array of
 513 *>                             dimension( LDU, max(NN) ).
 514 *>          The orthogonal matrix of eigenvectors computed by DSTEQR,
 515 *>          DPTEQR, and DSTEIN.
 516 *> \endverbatim
 517 *>
 518 *> \param[out] WORK
 519 *> \verbatim
 520 *>          WORK is DOUBLE PRECISION array of
 521 *>                      dimension( LWORK )
 522 *> \endverbatim
 523 *>
 524 *> \param[in] LWORK
 525 *> \verbatim
 526 *>          LWORK is INTEGER
 527 *>          The number of entries in WORK.  This must be at least
 528 *>          1 + 4 * Nmax + 2 * Nmax * lg Nmax + 3 * Nmax**2
 529 *>          where Nmax = max( NN(j), 2 ) and lg = log base 2.
 530 *> \endverbatim
 531 *>
 532 *> \param[out] IWORK
 533 *> \verbatim
 534 *>          IWORK is INTEGER array,
 535 *>          Workspace.
 536 *> \endverbatim
 537 *>
 538 *> \param[out] LIWORK
 539 *> \verbatim
 540 *>          LIWORK is INTEGER
 541 *>          The number of entries in IWORK.  This must be at least
 542 *>                  6 + 6*Nmax + 5 * Nmax * lg Nmax
 543 *>          where Nmax = max( NN(j), 2 ) and lg = log base 2.
 544 *> \endverbatim
 545 *>
 546 *> \param[out] RESULT
 547 *> \verbatim
 548 *>          RESULT is DOUBLE PRECISION array, dimension (26)
 549 *>          The values computed by the tests described above.
 550 *>          The values are currently limited to 1/ulp, to avoid
 551 *>          overflow.
 552 *> \endverbatim
 553 *>
 554 *> \param[out] INFO
 555 *> \verbatim
 556 *>          INFO is INTEGER
 557 *>          If 0, then everything ran OK.
 558 *>           -1: NSIZES < 0
 559 *>           -2: Some NN(j) < 0
 560 *>           -3: NTYPES < 0
 561 *>           -5: THRESH < 0
 562 *>           -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
 563 *>          -23: LDU < 1 or LDU < NMAX.
 564 *>          -29: LWORK too small.
 565 *>          If  DLATMR, SLATMS, DSYTRD, DORGTR, DSTEQR, SSTERF,
 566 *>              or DORMC2 returns an error code, the
 567 *>              absolute value of it is returned.
 568 *>
 569 *>-----------------------------------------------------------------------
 570 *>
 571 *>       Some Local Variables and Parameters:
 572 *>       ---- ----- --------- --- ----------
 573 *>       ZERO, ONE       Real 0 and 1.
 574 *>       MAXTYP          The number of types defined.
 575 *>       NTEST           The number of tests performed, or which can
 576 *>                       be performed so far, for the current matrix.
 577 *>       NTESTT          The total number of tests performed so far.
 578 *>       NBLOCK          Blocksize as returned by ENVIR.
 579 *>       NMAX            Largest value in NN.
 580 *>       NMATS           The number of matrices generated so far.
 581 *>       NERRS           The number of tests which have exceeded THRESH
 582 *>                       so far.
 583 *>       COND, IMODE     Values to be passed to the matrix generators.
 584 *>       ANORM           Norm of A; passed to matrix generators.
 585 *>
 586 *>       OVFL, UNFL      Overflow and underflow thresholds.
 587 *>       ULP, ULPINV     Finest relative precision and its inverse.
 588 *>       RTOVFL, RTUNFL  Square roots of the previous 2 values.
 589 *>               The following four arrays decode JTYPE:
 590 *>       KTYPE(j)        The general type (1-10) for type "j".
 591 *>       KMODE(j)        The MODE value to be passed to the matrix
 592 *>                       generator for type "j".
 593 *>       KMAGN(j)        The order of magnitude ( O(1),
 594 *>                       O(overflow^(1/2) ), O(underflow^(1/2) )
 595 *> \endverbatim
 596 *
 597 *  Authors:
 598 *  ========
 599 *
 600 *> \author Univ. of Tennessee
 601 *> \author Univ. of California Berkeley
 602 *> \author Univ. of Colorado Denver
 603 *> \author NAG Ltd.
 604 *
 605 *> \date December 2016
 606 *
 607 *> \ingroup double_eig
 608 *
 609 *  =====================================================================
 610       SUBROUTINE DCHKST2STG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
 611      $                   NOUNIT, A, LDA, AP, SD, SE, D1, D2, D3, D4, D5,
 612      $                   WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK,
 613      $                   LWORK, IWORK, LIWORK, RESULT, INFO )
 614 *
 615 *  -- LAPACK test routine (version 3.7.0) --
 616 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 617 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 618 *     December 2016
 619 *
 620 *     .. Scalar Arguments ..
 621       INTEGER            INFO, LDA, LDU, LIWORK, LWORK, NOUNIT, NSIZES,
 622      $                   NTYPES
 623       DOUBLE PRECISION   THRESH
 624 *     ..
 625 *     .. Array Arguments ..
 626       LOGICAL            DOTYPE( * )
 627       INTEGER            ISEED( 4 ), IWORK( * ), NN( * )
 628       DOUBLE PRECISION   A( LDA, * ), AP( * ), D1( * ), D2( * ),
 629      $                   D3( * ), D4( * ), D5( * ), RESULT( * ),
 630      $                   SD( * ), SE( * ), TAU( * ), U( LDU, * ),
 631      $                   V( LDU, * ), VP( * ), WA1( * ), WA2( * ),
 632      $                   WA3( * ), WORK( * ), WR( * ), Z( LDU, * )
 633 *     ..
 634 *
 635 *  =====================================================================
 636 *
 637 *     .. Parameters ..
 638       DOUBLE PRECISION   ZERO, ONE, TWO, EIGHT, TEN, HUN
 639       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
 640      $                   EIGHT = 8.0D0, TEN = 10.0D0, HUN = 100.0D0 )
 641       DOUBLE PRECISION   HALF
 642       PARAMETER          ( HALF = ONE / TWO )
 643       INTEGER            MAXTYP
 644       PARAMETER          ( MAXTYP = 21 )
 645       LOGICAL            SRANGE
 646       PARAMETER          ( SRANGE = .FALSE. )
 647       LOGICAL            SREL
 648       PARAMETER          ( SREL = .FALSE. )
 649 *     ..
 650 *     .. Local Scalars ..
 651       LOGICAL            BADNN, TRYRAC
 652       INTEGER            I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, J, JC,
 653      $                   JR, JSIZE, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC,
 654      $                   M, M2, M3, MTYPES, N, NAP, NBLOCK, NERRS,
 655      $                   NMATS, NMAX, NSPLIT, NTEST, NTESTT, LH, LW
 656       DOUBLE PRECISION   ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL,
 657      $                   RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP,
 658      $                   ULPINV, UNFL, VL, VU
 659 *     ..
 660 *     .. Local Arrays ..
 661       INTEGER            IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
 662      $                   KMAGN( MAXTYP ), KMODE( MAXTYP ),
 663      $                   KTYPE( MAXTYP )
 664       DOUBLE PRECISION   DUMMA( 1 )
 665 *     ..
 666 *     .. External Functions ..
 667       INTEGER            ILAENV
 668       DOUBLE PRECISION   DLAMCH, DLARND, DSXT1
 669       EXTERNAL           ILAENV, DLAMCH, DLARND, DSXT1
 670 *     ..
 671 *     .. External Subroutines ..
 672       EXTERNAL           DCOPY, DLABAD, DLACPY, DLASET, DLASUM, DLATMR,
 673      $                   DLATMS, DOPGTR, DORGTR, DPTEQR, DSPT21, DSPTRD,
 674      $                   DSTEBZ, DSTECH, DSTEDC, DSTEMR, DSTEIN, DSTEQR,
 675      $                   DSTERF, DSTT21, DSTT22, DSYT21, DSYTRD, XERBLA,
 676      $                   DSYTRD_2STAGE
 677 *     ..
 678 *     .. Intrinsic Functions ..
 679       INTRINSIC          ABS, DBLE, INT, LOG, MAX, MIN, SQRT
 680 *     ..
 681 *     .. Data statements ..
 682       DATA               KTYPE / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8,
 683      $                   8, 8, 9, 9, 9, 9, 9, 10 /
 684       DATA               KMAGN / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
 685      $                   2, 3, 1, 1, 1, 2, 3, 1 /
 686       DATA               KMODE / 0, 0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
 687      $                   0, 0, 4, 3, 1, 4, 4, 3 /
 688 *     ..
 689 *     .. Executable Statements ..
 690 *
 691 *     Keep ftnchek happy
 692       IDUMMA( 1 ) = 1
 693 *
 694 *     Check for errors
 695 *
 696       NTESTT = 0
 697       INFO = 0
 698 *
 699 *     Important constants
 700 *
 701       BADNN = .FALSE.
 702       TRYRAC = .TRUE.
 703       NMAX = 1
 704       DO 10 J = 1, NSIZES
 705          NMAX = MAX( NMAX, NN( J ) )
 706          IF( NN( J ).LT.0 )
 707      $      BADNN = .TRUE.
 708    10 CONTINUE
 709 *
 710       NBLOCK = ILAENV( 1, 'DSYTRD', 'L', NMAX, -1, -1, -1 )
 711       NBLOCK = MIN( NMAX, MAX( 1, NBLOCK ) )
 712 *
 713 *     Check for errors
 714 *
 715       IF( NSIZES.LT.0 ) THEN
 716          INFO = -1
 717       ELSE IF( BADNN ) THEN
 718          INFO = -2
 719       ELSE IF( NTYPES.LT.0 ) THEN
 720          INFO = -3
 721       ELSE IF( LDA.LT.NMAX ) THEN
 722          INFO = -9
 723       ELSE IF( LDU.LT.NMAX ) THEN
 724          INFO = -23
 725       ELSE IF( 2*MAX( 2, NMAX )**2.GT.LWORK ) THEN
 726          INFO = -29
 727       END IF
 728 *
 729       IF( INFO.NE.0 ) THEN
 730          CALL XERBLA( 'DCHKST2STG', -INFO )
 731          RETURN
 732       END IF
 733 *
 734 *     Quick return if possible
 735 *
 736       IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
 737      $   RETURN
 738 *
 739 *     More Important constants
 740 *
 741       UNFL = DLAMCH( 'Safe minimum' )
 742       OVFL = ONE / UNFL
 743       CALL DLABAD( UNFL, OVFL )
 744       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
 745       ULPINV = ONE / ULP
 746       LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) )
 747       RTUNFL = SQRT( UNFL )
 748       RTOVFL = SQRT( OVFL )
 749 *
 750 *     Loop over sizes, types
 751 *
 752       DO 20 I = 1, 4
 753          ISEED2( I ) = ISEED( I )
 754    20 CONTINUE
 755       NERRS = 0
 756       NMATS = 0
 757 *
 758       DO 310 JSIZE = 1, NSIZES
 759          N = NN( JSIZE )
 760          IF( N.GT.0 ) THEN
 761             LGN = INT( LOG( DBLE( N ) ) / LOG( TWO ) )
 762             IF( 2**LGN.LT.N )
 763      $         LGN = LGN + 1
 764             IF( 2**LGN.LT.N )
 765      $         LGN = LGN + 1
 766             LWEDC = 1 + 4*N + 2*N*LGN + 4*N**2
 767             LIWEDC = 6 + 6*N + 5*N*LGN
 768          ELSE
 769             LWEDC = 8
 770             LIWEDC = 12
 771          END IF
 772          NAP = ( N*( N+1 ) ) / 2
 773          ANINV = ONE / DBLE( MAX( 1, N ) )
 774 *
 775          IF( NSIZES.NE.1 ) THEN
 776             MTYPES = MIN( MAXTYP, NTYPES )
 777          ELSE
 778             MTYPES = MIN( MAXTYP+1, NTYPES )
 779          END IF
 780 *
 781          DO 300 JTYPE = 1, MTYPES
 782             IF( .NOT.DOTYPE( JTYPE ) )
 783      $         GO TO 300
 784             NMATS = NMATS + 1
 785             NTEST = 0
 786 *
 787             DO 30 J = 1, 4
 788                IOLDSD( J ) = ISEED( J )
 789    30       CONTINUE
 790 *
 791 *           Compute "A"
 792 *
 793 *           Control parameters:
 794 *
 795 *               KMAGN  KMODE        KTYPE
 796 *           =1  O(1)   clustered 1  zero
 797 *           =2  large  clustered 2  identity
 798 *           =3  small  exponential  (none)
 799 *           =4         arithmetic   diagonal, (w/ eigenvalues)
 800 *           =5         random log   symmetric, w/ eigenvalues
 801 *           =6         random       (none)
 802 *           =7                      random diagonal
 803 *           =8                      random symmetric
 804 *           =9                      positive definite
 805 *           =10                     diagonally dominant tridiagonal
 806 *
 807             IF( MTYPES.GT.MAXTYP )
 808      $         GO TO 100
 809 *
 810             ITYPE = KTYPE( JTYPE )
 811             IMODE = KMODE( JTYPE )
 812 *
 813 *           Compute norm
 814 *
 815             GO TO ( 40, 50, 60 )KMAGN( JTYPE )
 816 *
 817    40       CONTINUE
 818             ANORM = ONE
 819             GO TO 70
 820 *
 821    50       CONTINUE
 822             ANORM = ( RTOVFL*ULP )*ANINV
 823             GO TO 70
 824 *
 825    60       CONTINUE
 826             ANORM = RTUNFL*N*ULPINV
 827             GO TO 70
 828 *
 829    70       CONTINUE
 830 *
 831             CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
 832             IINFO = 0
 833             IF( JTYPE.LE.15 ) THEN
 834                COND = ULPINV
 835             ELSE
 836                COND = ULPINV*ANINV / TEN
 837             END IF
 838 *
 839 *           Special Matrices -- Identity & Jordan block
 840 *
 841 *              Zero
 842 *
 843             IF( ITYPE.EQ.1 ) THEN
 844                IINFO = 0
 845 *
 846             ELSE IF( ITYPE.EQ.2 ) THEN
 847 *
 848 *              Identity
 849 *
 850                DO 80 JC = 1, N
 851                   A( JC, JC ) = ANORM
 852    80          CONTINUE
 853 *
 854             ELSE IF( ITYPE.EQ.4 ) THEN
 855 *
 856 *              Diagonal Matrix, [Eigen]values Specified
 857 *
 858                CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
 859      $                      ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
 860      $                      IINFO )
 861 *
 862 *
 863             ELSE IF( ITYPE.EQ.5 ) THEN
 864 *
 865 *              Symmetric, eigenvalues specified
 866 *
 867                CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
 868      $                      ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
 869      $                      IINFO )
 870 *
 871             ELSE IF( ITYPE.EQ.7 ) THEN
 872 *
 873 *              Diagonal, random eigenvalues
 874 *
 875                CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
 876      $                      'T', 'N', WORK( N+1 ), 1, ONE,
 877      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
 878      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
 879 *
 880             ELSE IF( ITYPE.EQ.8 ) THEN
 881 *
 882 *              Symmetric, random eigenvalues
 883 *
 884                CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
 885      $                      'T', 'N', WORK( N+1 ), 1, ONE,
 886      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
 887      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
 888 *
 889             ELSE IF( ITYPE.EQ.9 ) THEN
 890 *
 891 *              Positive definite, eigenvalues specified.
 892 *
 893                CALL DLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND,
 894      $                      ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
 895      $                      IINFO )
 896 *
 897             ELSE IF( ITYPE.EQ.10 ) THEN
 898 *
 899 *              Positive definite tridiagonal, eigenvalues specified.
 900 *
 901                CALL DLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND,
 902      $                      ANORM, 1, 1, 'N', A, LDA, WORK( N+1 ),
 903      $                      IINFO )
 904                DO 90 I = 2, N
 905                   TEMP1 = ABS( A( I-1, I ) ) /
 906      $                    SQRT( ABS( A( I-1, I-1 )*A( I, I ) ) )
 907                   IF( TEMP1.GT.HALF ) THEN
 908                      A( I-1, I ) = HALF*SQRT( ABS( A( I-1, I-1 )*A( I,
 909      $                             I ) ) )
 910                      A( I, I-1 ) = A( I-1, I )
 911                   END IF
 912    90          CONTINUE
 913 *
 914             ELSE
 915 *
 916                IINFO = 1
 917             END IF
 918 *
 919             IF( IINFO.NE.0 ) THEN
 920                WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
 921      $            IOLDSD
 922                INFO = ABS( IINFO )
 923                RETURN
 924             END IF
 925 *
 926   100       CONTINUE
 927 *
 928 *           Call DSYTRD and DORGTR to compute S and U from
 929 *           upper triangle.
 930 *
 931             CALL DLACPY( 'U', N, N, A, LDA, V, LDU )
 932 *
 933             NTEST = 1
 934             CALL DSYTRD( 'U', N, V, LDU, SD, SE, TAU, WORK, LWORK,
 935      $                   IINFO )
 936 *
 937             IF( IINFO.NE.0 ) THEN
 938                WRITE( NOUNIT, FMT = 9999 )'DSYTRD(U)', IINFO, N, JTYPE,
 939      $            IOLDSD
 940                INFO = ABS( IINFO )
 941                IF( IINFO.LT.0 ) THEN
 942                   RETURN
 943                ELSE
 944                   RESULT( 1 ) = ULPINV
 945                   GO TO 280
 946                END IF
 947             END IF
 948 *
 949             CALL DLACPY( 'U', N, N, V, LDU, U, LDU )
 950 *
 951             NTEST = 2
 952             CALL DORGTR( 'U', N, U, LDU, TAU, WORK, LWORK, IINFO )
 953             IF( IINFO.NE.0 ) THEN
 954                WRITE( NOUNIT, FMT = 9999 )'DORGTR(U)', IINFO, N, JTYPE,
 955      $            IOLDSD
 956                INFO = ABS( IINFO )
 957                IF( IINFO.LT.0 ) THEN
 958                   RETURN
 959                ELSE
 960                   RESULT( 2 ) = ULPINV
 961                   GO TO 280
 962                END IF
 963             END IF
 964 *
 965 *           Do tests 1 and 2
 966 *
 967             CALL DSYT21( 2, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V,
 968      $                   LDU, TAU, WORK, RESULT( 1 ) )
 969             CALL DSYT21( 3, 'Upper', N, 1, A, LDA, SD, SE, U, LDU, V,
 970      $                   LDU, TAU, WORK, RESULT( 2 ) )
 971 *
 972 *           Compute D1 the eigenvalues resulting from the tridiagonal
 973 *           form using the standard 1-stage algorithm and use it as a
 974 *           reference to compare with the 2-stage technique
 975 *
 976 *           Compute D1 from the 1-stage and used as reference for the
 977 *           2-stage
 978 *
 979             CALL DCOPY( N, SD, 1, D1, 1 )
 980             IF( N.GT.0 )
 981      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
 982 *
 983             CALL DSTEQR( 'N', N, D1, WORK, WORK( N+1 ), LDU,
 984      $                   WORK( N+1 ), IINFO )
 985             IF( IINFO.NE.0 ) THEN
 986                WRITE( NOUNIT, FMT = 9999 )'DSTEQR(N)', IINFO, N, JTYPE,
 987      $            IOLDSD
 988                INFO = ABS( IINFO )
 989                IF( IINFO.LT.0 ) THEN
 990                   RETURN
 991                ELSE
 992                   RESULT( 3 ) = ULPINV
 993                   GO TO 280
 994                END IF
 995             END IF
 996 *
 997 *           2-STAGE TRD Upper case is used to compute D2.
 998 *           Note to set SD and SE to zero to be sure not reusing
 999 *           the one from above. Compare it with D1 computed
1000 *           using the 1-stage.
1001 *
1002             CALL DLASET( 'Full', N, 1, ZERO, ZERO, SD, N )
1003             CALL DLASET( 'Full', N, 1, ZERO, ZERO, SE, N )
1004             CALL DLACPY( "U", N, N, A, LDA, V, LDU )
1005             LH = MAX(1, 4*N)
1006             LW = LWORK - LH
1007             CALL DSYTRD_2STAGE( 'N', "U", N, V, LDU, SD, SE, TAU,
1008      $                   WORK, LH, WORK( LH+1 ), LW, IINFO )
1009 *
1010 *           Compute D2 from the 2-stage Upper case
1011 *
1012             CALL DCOPY( N, SD, 1, D2, 1 )
1013             IF( N.GT.0 )
1014      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
1015 *
1016             CALL DSTEQR( 'N', N, D2, WORK, WORK( N+1 ), LDU,
1017      $                   WORK( N+1 ), IINFO )
1018             IF( IINFO.NE.0 ) THEN
1019                WRITE( NOUNIT, FMT = 9999 )'DSTEQR(N)', IINFO, N, JTYPE,
1020      $            IOLDSD
1021                INFO = ABS( IINFO )
1022                IF( IINFO.LT.0 ) THEN
1023                   RETURN
1024                ELSE
1025                   RESULT( 3 ) = ULPINV
1026                   GO TO 280
1027                END IF
1028             END IF
1029 *
1030 *           2-STAGE TRD Lower case is used to compute D3.
1031 *           Note to set SD and SE to zero to be sure not reusing
1032 *           the one from above. Compare it with D1 computed
1033 *           using the 1-stage.
1034 *
1035             CALL DLASET( 'Full', N, 1, ZERO, ZERO, SD, N )
1036             CALL DLASET( 'Full', N, 1, ZERO, ZERO, SE, N )
1037             CALL DLACPY( "L", N, N, A, LDA, V, LDU )
1038             CALL DSYTRD_2STAGE( 'N', "L", N, V, LDU, SD, SE, TAU,
1039      $                   WORK, LH, WORK( LH+1 ), LW, IINFO )
1040 *
1041 *           Compute D3 from the 2-stage Upper case
1042 *
1043             CALL DCOPY( N, SD, 1, D3, 1 )
1044             IF( N.GT.0 )
1045      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
1046 *
1047             CALL DSTEQR( 'N', N, D3, WORK, WORK( N+1 ), LDU,
1048      $                   WORK( N+1 ), IINFO )
1049             IF( IINFO.NE.0 ) THEN
1050                WRITE( NOUNIT, FMT = 9999 )'DSTEQR(N)', IINFO, N, JTYPE,
1051      $            IOLDSD
1052                INFO = ABS( IINFO )
1053                IF( IINFO.LT.0 ) THEN
1054                   RETURN
1055                ELSE
1056                   RESULT( 4 ) = ULPINV
1057                   GO TO 280
1058                END IF
1059             END IF
1060 *
1061 *
1062 *           Do Tests 3 and 4 which are similar to 11 and 12 but with the
1063 *           D1 computed using the standard 1-stage reduction as reference
1064 *
1065             NTEST = 4
1066             TEMP1 = ZERO
1067             TEMP2 = ZERO
1068             TEMP3 = ZERO
1069             TEMP4 = ZERO
1070 *
1071             DO 151 J = 1, N
1072                TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
1073                TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
1074                TEMP3 = MAX( TEMP3, ABS( D1( J ) ), ABS( D3( J ) ) )
1075                TEMP4 = MAX( TEMP4, ABS( D1( J )-D3( J ) ) )
1076   151       CONTINUE
1077 *
1078             RESULT( 3 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
1079             RESULT( 4 ) = TEMP4 / MAX( UNFL, ULP*MAX( TEMP3, TEMP4 ) )
1080 *
1081 *           Store the upper triangle of A in AP
1082 *
1083             I = 0
1084             DO 120 JC = 1, N
1085                DO 110 JR = 1, JC
1086                   I = I + 1
1087                   AP( I ) = A( JR, JC )
1088   110          CONTINUE
1089   120       CONTINUE
1090 *
1091 *           Call DSPTRD and DOPGTR to compute S and U from AP
1092 *
1093             CALL DCOPY( NAP, AP, 1, VP, 1 )
1094 *
1095             NTEST = 5
1096             CALL DSPTRD( 'U', N, VP, SD, SE, TAU, IINFO )
1097 *
1098             IF( IINFO.NE.0 ) THEN
1099                WRITE( NOUNIT, FMT = 9999 )'DSPTRD(U)', IINFO, N, JTYPE,
1100      $            IOLDSD
1101                INFO = ABS( IINFO )
1102                IF( IINFO.LT.0 ) THEN
1103                   RETURN
1104                ELSE
1105                   RESULT( 5 ) = ULPINV
1106                   GO TO 280
1107                END IF
1108             END IF
1109 *
1110             NTEST = 6
1111             CALL DOPGTR( 'U', N, VP, TAU, U, LDU, WORK, IINFO )
1112             IF( IINFO.NE.0 ) THEN
1113                WRITE( NOUNIT, FMT = 9999 )'DOPGTR(U)', IINFO, N, JTYPE,
1114      $            IOLDSD
1115                INFO = ABS( IINFO )
1116                IF( IINFO.LT.0 ) THEN
1117                   RETURN
1118                ELSE
1119                   RESULT( 6 ) = ULPINV
1120                   GO TO 280
1121                END IF
1122             END IF
1123 *
1124 *           Do tests 5 and 6
1125 *
1126             CALL DSPT21( 2, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU,
1127      $                   WORK, RESULT( 5 ) )
1128             CALL DSPT21( 3, 'Upper', N, 1, AP, SD, SE, U, LDU, VP, TAU,
1129      $                   WORK, RESULT( 6 ) )
1130 *
1131 *           Store the lower triangle of A in AP
1132 *
1133             I = 0
1134             DO 140 JC = 1, N
1135                DO 130 JR = JC, N
1136                   I = I + 1
1137                   AP( I ) = A( JR, JC )
1138   130          CONTINUE
1139   140       CONTINUE
1140 *
1141 *           Call DSPTRD and DOPGTR to compute S and U from AP
1142 *
1143             CALL DCOPY( NAP, AP, 1, VP, 1 )
1144 *
1145             NTEST = 7
1146             CALL DSPTRD( 'L', N, VP, SD, SE, TAU, IINFO )
1147 *
1148             IF( IINFO.NE.0 ) THEN
1149                WRITE( NOUNIT, FMT = 9999 )'DSPTRD(L)', IINFO, N, JTYPE,
1150      $            IOLDSD
1151                INFO = ABS( IINFO )
1152                IF( IINFO.LT.0 ) THEN
1153                   RETURN
1154                ELSE
1155                   RESULT( 7 ) = ULPINV
1156                   GO TO 280
1157                END IF
1158             END IF
1159 *
1160             NTEST = 8
1161             CALL DOPGTR( 'L', N, VP, TAU, U, LDU, WORK, IINFO )
1162             IF( IINFO.NE.0 ) THEN
1163                WRITE( NOUNIT, FMT = 9999 )'DOPGTR(L)', IINFO, N, JTYPE,
1164      $            IOLDSD
1165                INFO = ABS( IINFO )
1166                IF( IINFO.LT.0 ) THEN
1167                   RETURN
1168                ELSE
1169                   RESULT( 8 ) = ULPINV
1170                   GO TO 280
1171                END IF
1172             END IF
1173 *
1174             CALL DSPT21( 2, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU,
1175      $                   WORK, RESULT( 7 ) )
1176             CALL DSPT21( 3, 'Lower', N, 1, AP, SD, SE, U, LDU, VP, TAU,
1177      $                   WORK, RESULT( 8 ) )
1178 *
1179 *           Call DSTEQR to compute D1, D2, and Z, do tests.
1180 *
1181 *           Compute D1 and Z
1182 *
1183             CALL DCOPY( N, SD, 1, D1, 1 )
1184             IF( N.GT.0 )
1185      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
1186             CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
1187 *
1188             NTEST = 9
1189             CALL DSTEQR( 'V', N, D1, WORK, Z, LDU, WORK( N+1 ), IINFO )
1190             IF( IINFO.NE.0 ) THEN
1191                WRITE( NOUNIT, FMT = 9999 )'DSTEQR(V)', IINFO, N, JTYPE,
1192      $            IOLDSD
1193                INFO = ABS( IINFO )
1194                IF( IINFO.LT.0 ) THEN
1195                   RETURN
1196                ELSE
1197                   RESULT( 9 ) = ULPINV
1198                   GO TO 280
1199                END IF
1200             END IF
1201 *
1202 *           Compute D2
1203 *
1204             CALL DCOPY( N, SD, 1, D2, 1 )
1205             IF( N.GT.0 )
1206      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
1207 *
1208             NTEST = 11
1209             CALL DSTEQR( 'N', N, D2, WORK, WORK( N+1 ), LDU,
1210      $                   WORK( N+1 ), IINFO )
1211             IF( IINFO.NE.0 ) THEN
1212                WRITE( NOUNIT, FMT = 9999 )'DSTEQR(N)', IINFO, N, JTYPE,
1213      $            IOLDSD
1214                INFO = ABS( IINFO )
1215                IF( IINFO.LT.0 ) THEN
1216                   RETURN
1217                ELSE
1218                   RESULT( 11 ) = ULPINV
1219                   GO TO 280
1220                END IF
1221             END IF
1222 *
1223 *           Compute D3 (using PWK method)
1224 *
1225             CALL DCOPY( N, SD, 1, D3, 1 )
1226             IF( N.GT.0 )
1227      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
1228 *
1229             NTEST = 12
1230             CALL DSTERF( N, D3, WORK, IINFO )
1231             IF( IINFO.NE.0 ) THEN
1232                WRITE( NOUNIT, FMT = 9999 )'DSTERF', IINFO, N, JTYPE,
1233      $            IOLDSD
1234                INFO = ABS( IINFO )
1235                IF( IINFO.LT.0 ) THEN
1236                   RETURN
1237                ELSE
1238                   RESULT( 12 ) = ULPINV
1239                   GO TO 280
1240                END IF
1241             END IF
1242 *
1243 *           Do Tests 9 and 10
1244 *
1245             CALL DSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
1246      $                   RESULT( 9 ) )
1247 *
1248 *           Do Tests 11 and 12
1249 *
1250             TEMP1 = ZERO
1251             TEMP2 = ZERO
1252             TEMP3 = ZERO
1253             TEMP4 = ZERO
1254 *
1255             DO 150 J = 1, N
1256                TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
1257                TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
1258                TEMP3 = MAX( TEMP3, ABS( D1( J ) ), ABS( D3( J ) ) )
1259                TEMP4 = MAX( TEMP4, ABS( D1( J )-D3( J ) ) )
1260   150       CONTINUE
1261 *
1262             RESULT( 11 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
1263             RESULT( 12 ) = TEMP4 / MAX( UNFL, ULP*MAX( TEMP3, TEMP4 ) )
1264 *
1265 *           Do Test 13 -- Sturm Sequence Test of Eigenvalues
1266 *                         Go up by factors of two until it succeeds
1267 *
1268             NTEST = 13
1269             TEMP1 = THRESH*( HALF-ULP )
1270 *
1271             DO 160 J = 0, LOG2UI
1272                CALL DSTECH( N, SD, SE, D1, TEMP1, WORK, IINFO )
1273                IF( IINFO.EQ.0 )
1274      $            GO TO 170
1275                TEMP1 = TEMP1*TWO
1276   160       CONTINUE
1277 *
1278   170       CONTINUE
1279             RESULT( 13 ) = TEMP1
1280 *
1281 *           For positive definite matrices ( JTYPE.GT.15 ) call DPTEQR
1282 *           and do tests 14, 15, and 16 .
1283 *
1284             IF( JTYPE.GT.15 ) THEN
1285 *
1286 *              Compute D4 and Z4
1287 *
1288                CALL DCOPY( N, SD, 1, D4, 1 )
1289                IF( N.GT.0 )
1290      $            CALL DCOPY( N-1, SE, 1, WORK, 1 )
1291                CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
1292 *
1293                NTEST = 14
1294                CALL DPTEQR( 'V', N, D4, WORK, Z, LDU, WORK( N+1 ),
1295      $                      IINFO )
1296                IF( IINFO.NE.0 ) THEN
1297                   WRITE( NOUNIT, FMT = 9999 )'DPTEQR(V)', IINFO, N,
1298      $               JTYPE, IOLDSD
1299                   INFO = ABS( IINFO )
1300                   IF( IINFO.LT.0 ) THEN
1301                      RETURN
1302                   ELSE
1303                      RESULT( 14 ) = ULPINV
1304                      GO TO 280
1305                   END IF
1306                END IF
1307 *
1308 *              Do Tests 14 and 15
1309 *
1310                CALL DSTT21( N, 0, SD, SE, D4, DUMMA, Z, LDU, WORK,
1311      $                      RESULT( 14 ) )
1312 *
1313 *              Compute D5
1314 *
1315                CALL DCOPY( N, SD, 1, D5, 1 )
1316                IF( N.GT.0 )
1317      $            CALL DCOPY( N-1, SE, 1, WORK, 1 )
1318 *
1319                NTEST = 16
1320                CALL DPTEQR( 'N', N, D5, WORK, Z, LDU, WORK( N+1 ),
1321      $                      IINFO )
1322                IF( IINFO.NE.0 ) THEN
1323                   WRITE( NOUNIT, FMT = 9999 )'DPTEQR(N)', IINFO, N,
1324      $               JTYPE, IOLDSD
1325                   INFO = ABS( IINFO )
1326                   IF( IINFO.LT.0 ) THEN
1327                      RETURN
1328                   ELSE
1329                      RESULT( 16 ) = ULPINV
1330                      GO TO 280
1331                   END IF
1332                END IF
1333 *
1334 *              Do Test 16
1335 *
1336                TEMP1 = ZERO
1337                TEMP2 = ZERO
1338                DO 180 J = 1, N
1339                   TEMP1 = MAX( TEMP1, ABS( D4( J ) ), ABS( D5( J ) ) )
1340                   TEMP2 = MAX( TEMP2, ABS( D4( J )-D5( J ) ) )
1341   180          CONTINUE
1342 *
1343                RESULT( 16 ) = TEMP2 / MAX( UNFL,
1344      $                        HUN*ULP*MAX( TEMP1, TEMP2 ) )
1345             ELSE
1346                RESULT( 14 ) = ZERO
1347                RESULT( 15 ) = ZERO
1348                RESULT( 16 ) = ZERO
1349             END IF
1350 *
1351 *           Call DSTEBZ with different options and do tests 17-18.
1352 *
1353 *              If S is positive definite and diagonally dominant,
1354 *              ask for all eigenvalues with high relative accuracy.
1355 *
1356             VL = ZERO
1357             VU = ZERO
1358             IL = 0
1359             IU = 0
1360             IF( JTYPE.EQ.21 ) THEN
1361                NTEST = 17
1362                ABSTOL = UNFL + UNFL
1363                CALL DSTEBZ( 'A', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE,
1364      $                      M, NSPLIT, WR, IWORK( 1 ), IWORK( N+1 ),
1365      $                      WORK, IWORK( 2*N+1 ), IINFO )
1366                IF( IINFO.NE.0 ) THEN
1367                   WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(A,rel)', IINFO, N,
1368      $               JTYPE, IOLDSD
1369                   INFO = ABS( IINFO )
1370                   IF( IINFO.LT.0 ) THEN
1371                      RETURN
1372                   ELSE
1373                      RESULT( 17 ) = ULPINV
1374                      GO TO 280
1375                   END IF
1376                END IF
1377 *
1378 *              Do test 17
1379 *
1380                TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) /
1381      $                 ( ONE-HALF )**4
1382 *
1383                TEMP1 = ZERO
1384                DO 190 J = 1, N
1385                   TEMP1 = MAX( TEMP1, ABS( D4( J )-WR( N-J+1 ) ) /
1386      $                    ( ABSTOL+ABS( D4( J ) ) ) )
1387   190          CONTINUE
1388 *
1389                RESULT( 17 ) = TEMP1 / TEMP2
1390             ELSE
1391                RESULT( 17 ) = ZERO
1392             END IF
1393 *
1394 *           Now ask for all eigenvalues with high absolute accuracy.
1395 *
1396             NTEST = 18
1397             ABSTOL = UNFL + UNFL
1398             CALL DSTEBZ( 'A', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE, M,
1399      $                   NSPLIT, WA1, IWORK( 1 ), IWORK( N+1 ), WORK,
1400      $                   IWORK( 2*N+1 ), IINFO )
1401             IF( IINFO.NE.0 ) THEN
1402                WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(A)', IINFO, N, JTYPE,
1403      $            IOLDSD
1404                INFO = ABS( IINFO )
1405                IF( IINFO.LT.0 ) THEN
1406                   RETURN
1407                ELSE
1408                   RESULT( 18 ) = ULPINV
1409                   GO TO 280
1410                END IF
1411             END IF
1412 *
1413 *           Do test 18
1414 *
1415             TEMP1 = ZERO
1416             TEMP2 = ZERO
1417             DO 200 J = 1, N
1418                TEMP1 = MAX( TEMP1, ABS( D3( J ) ), ABS( WA1( J ) ) )
1419                TEMP2 = MAX( TEMP2, ABS( D3( J )-WA1( J ) ) )
1420   200       CONTINUE
1421 *
1422             RESULT( 18 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
1423 *
1424 *           Choose random values for IL and IU, and ask for the
1425 *           IL-th through IU-th eigenvalues.
1426 *
1427             NTEST = 19
1428             IF( N.LE.1 ) THEN
1429                IL = 1
1430                IU = N
1431             ELSE
1432                IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
1433                IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
1434                IF( IU.LT.IL ) THEN
1435                   ITEMP = IU
1436                   IU = IL
1437                   IL = ITEMP
1438                END IF
1439             END IF
1440 *
1441             CALL DSTEBZ( 'I', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE,
1442      $                   M2, NSPLIT, WA2, IWORK( 1 ), IWORK( N+1 ),
1443      $                   WORK, IWORK( 2*N+1 ), IINFO )
1444             IF( IINFO.NE.0 ) THEN
1445                WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(I)', IINFO, N, JTYPE,
1446      $            IOLDSD
1447                INFO = ABS( IINFO )
1448                IF( IINFO.LT.0 ) THEN
1449                   RETURN
1450                ELSE
1451                   RESULT( 19 ) = ULPINV
1452                   GO TO 280
1453                END IF
1454             END IF
1455 *
1456 *           Determine the values VL and VU of the IL-th and IU-th
1457 *           eigenvalues and ask for all eigenvalues in this range.
1458 *
1459             IF( N.GT.0 ) THEN
1460                IF( IL.NE.1 ) THEN
1461                   VL = WA1( IL ) - MAX( HALF*( WA1( IL )-WA1( IL-1 ) ),
1462      $                 ULP*ANORM, TWO*RTUNFL )
1463                ELSE
1464                   VL = WA1( 1 ) - MAX( HALF*( WA1( N )-WA1( 1 ) ),
1465      $                 ULP*ANORM, TWO*RTUNFL )
1466                END IF
1467                IF( IU.NE.N ) THEN
1468                   VU = WA1( IU ) + MAX( HALF*( WA1( IU+1 )-WA1( IU ) ),
1469      $                 ULP*ANORM, TWO*RTUNFL )
1470                ELSE
1471                   VU = WA1( N ) + MAX( HALF*( WA1( N )-WA1( 1 ) ),
1472      $                 ULP*ANORM, TWO*RTUNFL )
1473                END IF
1474             ELSE
1475                VL = ZERO
1476                VU = ONE
1477             END IF
1478 *
1479             CALL DSTEBZ( 'V', 'E', N, VL, VU, IL, IU, ABSTOL, SD, SE,
1480      $                   M3, NSPLIT, WA3, IWORK( 1 ), IWORK( N+1 ),
1481      $                   WORK, IWORK( 2*N+1 ), IINFO )
1482             IF( IINFO.NE.0 ) THEN
1483                WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(V)', IINFO, N, JTYPE,
1484      $            IOLDSD
1485                INFO = ABS( IINFO )
1486                IF( IINFO.LT.0 ) THEN
1487                   RETURN
1488                ELSE
1489                   RESULT( 19 ) = ULPINV
1490                   GO TO 280
1491                END IF
1492             END IF
1493 *
1494             IF( M3.EQ.0 .AND. N.NE.0 ) THEN
1495                RESULT( 19 ) = ULPINV
1496                GO TO 280
1497             END IF
1498 *
1499 *           Do test 19
1500 *
1501             TEMP1 = DSXT1( 1, WA2, M2, WA3, M3, ABSTOL, ULP, UNFL )
1502             TEMP2 = DSXT1( 1, WA3, M3, WA2, M2, ABSTOL, ULP, UNFL )
1503             IF( N.GT.0 ) THEN
1504                TEMP3 = MAX( ABS( WA1( N ) ), ABS( WA1( 1 ) ) )
1505             ELSE
1506                TEMP3 = ZERO
1507             END IF
1508 *
1509             RESULT( 19 ) = ( TEMP1+TEMP2 ) / MAX( UNFL, TEMP3*ULP )
1510 *
1511 *           Call DSTEIN to compute eigenvectors corresponding to
1512 *           eigenvalues in WA1.  (First call DSTEBZ again, to make sure
1513 *           it returns these eigenvalues in the correct order.)
1514 *
1515             NTEST = 21
1516             CALL DSTEBZ( 'A', 'B', N, VL, VU, IL, IU, ABSTOL, SD, SE, M,
1517      $                   NSPLIT, WA1, IWORK( 1 ), IWORK( N+1 ), WORK,
1518      $                   IWORK( 2*N+1 ), IINFO )
1519             IF( IINFO.NE.0 ) THEN
1520                WRITE( NOUNIT, FMT = 9999 )'DSTEBZ(A,B)', IINFO, N,
1521      $            JTYPE, IOLDSD
1522                INFO = ABS( IINFO )
1523                IF( IINFO.LT.0 ) THEN
1524                   RETURN
1525                ELSE
1526                   RESULT( 20 ) = ULPINV
1527                   RESULT( 21 ) = ULPINV
1528                   GO TO 280
1529                END IF
1530             END IF
1531 *
1532             CALL DSTEIN( N, SD, SE, M, WA1, IWORK( 1 ), IWORK( N+1 ), Z,
1533      $                   LDU, WORK, IWORK( 2*N+1 ), IWORK( 3*N+1 ),
1534      $                   IINFO )
1535             IF( IINFO.NE.0 ) THEN
1536                WRITE( NOUNIT, FMT = 9999 )'DSTEIN', IINFO, N, JTYPE,
1537      $            IOLDSD
1538                INFO = ABS( IINFO )
1539                IF( IINFO.LT.0 ) THEN
1540                   RETURN
1541                ELSE
1542                   RESULT( 20 ) = ULPINV
1543                   RESULT( 21 ) = ULPINV
1544                   GO TO 280
1545                END IF
1546             END IF
1547 *
1548 *           Do tests 20 and 21
1549 *
1550             CALL DSTT21( N, 0, SD, SE, WA1, DUMMA, Z, LDU, WORK,
1551      $                   RESULT( 20 ) )
1552 *
1553 *           Call DSTEDC(I) to compute D1 and Z, do tests.
1554 *
1555 *           Compute D1 and Z
1556 *
1557             CALL DCOPY( N, SD, 1, D1, 1 )
1558             IF( N.GT.0 )
1559      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
1560             CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
1561 *
1562             NTEST = 22
1563             CALL DSTEDC( 'I', N, D1, WORK, Z, LDU, WORK( N+1 ), LWEDC-N,
1564      $                   IWORK, LIWEDC, IINFO )
1565             IF( IINFO.NE.0 ) THEN
1566                WRITE( NOUNIT, FMT = 9999 )'DSTEDC(I)', IINFO, N, JTYPE,
1567      $            IOLDSD
1568                INFO = ABS( IINFO )
1569                IF( IINFO.LT.0 ) THEN
1570                   RETURN
1571                ELSE
1572                   RESULT( 22 ) = ULPINV
1573                   GO TO 280
1574                END IF
1575             END IF
1576 *
1577 *           Do Tests 22 and 23
1578 *
1579             CALL DSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
1580      $                   RESULT( 22 ) )
1581 *
1582 *           Call DSTEDC(V) to compute D1 and Z, do tests.
1583 *
1584 *           Compute D1 and Z
1585 *
1586             CALL DCOPY( N, SD, 1, D1, 1 )
1587             IF( N.GT.0 )
1588      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
1589             CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
1590 *
1591             NTEST = 24
1592             CALL DSTEDC( 'V', N, D1, WORK, Z, LDU, WORK( N+1 ), LWEDC-N,
1593      $                   IWORK, LIWEDC, IINFO )
1594             IF( IINFO.NE.0 ) THEN
1595                WRITE( NOUNIT, FMT = 9999 )'DSTEDC(V)', IINFO, N, JTYPE,
1596      $            IOLDSD
1597                INFO = ABS( IINFO )
1598                IF( IINFO.LT.0 ) THEN
1599                   RETURN
1600                ELSE
1601                   RESULT( 24 ) = ULPINV
1602                   GO TO 280
1603                END IF
1604             END IF
1605 *
1606 *           Do Tests 24 and 25
1607 *
1608             CALL DSTT21( N, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
1609      $                   RESULT( 24 ) )
1610 *
1611 *           Call DSTEDC(N) to compute D2, do tests.
1612 *
1613 *           Compute D2
1614 *
1615             CALL DCOPY( N, SD, 1, D2, 1 )
1616             IF( N.GT.0 )
1617      $         CALL DCOPY( N-1, SE, 1, WORK, 1 )
1618             CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
1619 *
1620             NTEST = 26
1621             CALL DSTEDC( 'N', N, D2, WORK, Z, LDU, WORK( N+1 ), LWEDC-N,
1622      $                   IWORK, LIWEDC, IINFO )
1623             IF( IINFO.NE.0 ) THEN
1624                WRITE( NOUNIT, FMT = 9999 )'DSTEDC(N)', IINFO, N, JTYPE,
1625      $            IOLDSD
1626                INFO = ABS( IINFO )
1627                IF( IINFO.LT.0 ) THEN
1628                   RETURN
1629                ELSE
1630                   RESULT( 26 ) = ULPINV
1631                   GO TO 280
1632                END IF
1633             END IF
1634 *
1635 *           Do Test 26
1636 *
1637             TEMP1 = ZERO
1638             TEMP2 = ZERO
1639 *
1640             DO 210 J = 1, N
1641                TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
1642                TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
1643   210       CONTINUE
1644 *
1645             RESULT( 26 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
1646 *
1647 *           Only test DSTEMR if IEEE compliant
1648 *
1649             IF( ILAENV( 10, 'DSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND.
1650      $          ILAENV( 11, 'DSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN
1651 *
1652 *           Call DSTEMR, do test 27 (relative eigenvalue accuracy)
1653 *
1654 *              If S is positive definite and diagonally dominant,
1655 *              ask for all eigenvalues with high relative accuracy.
1656 *
1657                VL = ZERO
1658                VU = ZERO
1659                IL = 0
1660                IU = 0
1661                IF( JTYPE.EQ.21 .AND. SREL ) THEN
1662                   NTEST = 27
1663                   ABSTOL = UNFL + UNFL
1664                   CALL DSTEMR( 'V', 'A', N, SD, SE, VL, VU, IL, IU,
1665      $                         M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC,
1666      $                         WORK, LWORK, IWORK( 2*N+1 ), LWORK-2*N,
1667      $                         IINFO )
1668                   IF( IINFO.NE.0 ) THEN
1669                      WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,A,rel)',
1670      $                  IINFO, N, JTYPE, IOLDSD
1671                      INFO = ABS( IINFO )
1672                      IF( IINFO.LT.0 ) THEN
1673                         RETURN
1674                      ELSE
1675                         RESULT( 27 ) = ULPINV
1676                         GO TO 270
1677                      END IF
1678                   END IF
1679 *
1680 *              Do test 27
1681 *
1682                   TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) /
1683      $                    ( ONE-HALF )**4
1684 *
1685                   TEMP1 = ZERO
1686                   DO 220 J = 1, N
1687                      TEMP1 = MAX( TEMP1, ABS( D4( J )-WR( N-J+1 ) ) /
1688      $                       ( ABSTOL+ABS( D4( J ) ) ) )
1689   220             CONTINUE
1690 *
1691                   RESULT( 27 ) = TEMP1 / TEMP2
1692 *
1693                   IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
1694                   IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
1695                   IF( IU.LT.IL ) THEN
1696                      ITEMP = IU
1697                      IU = IL
1698                      IL = ITEMP
1699                   END IF
1700 *
1701                   IF( SRANGE ) THEN
1702                      NTEST = 28
1703                      ABSTOL = UNFL + UNFL
1704                      CALL DSTEMR( 'V', 'I', N, SD, SE, VL, VU, IL, IU,
1705      $                            M, WR, Z, LDU, N, IWORK( 1 ), TRYRAC,
1706      $                            WORK, LWORK, IWORK( 2*N+1 ),
1707      $                            LWORK-2*N, IINFO )
1708 *
1709                      IF( IINFO.NE.0 ) THEN
1710                         WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,I,rel)',
1711      $                     IINFO, N, JTYPE, IOLDSD
1712                         INFO = ABS( IINFO )
1713                         IF( IINFO.LT.0 ) THEN
1714                            RETURN
1715                         ELSE
1716                            RESULT( 28 ) = ULPINV
1717                            GO TO 270
1718                         END IF
1719                      END IF
1720 *
1721 *
1722 *                 Do test 28
1723 *
1724                      TEMP2 = TWO*( TWO*N-ONE )*ULP*
1725      $                       ( ONE+EIGHT*HALF**2 ) / ( ONE-HALF )**4
1726 *
1727                      TEMP1 = ZERO
1728                      DO 230 J = IL, IU
1729                         TEMP1 = MAX( TEMP1, ABS( WR( J-IL+1 )-D4( N-J+
1730      $                          1 ) ) / ( ABSTOL+ABS( WR( J-IL+1 ) ) ) )
1731   230                CONTINUE
1732 *
1733                      RESULT( 28 ) = TEMP1 / TEMP2
1734                   ELSE
1735                      RESULT( 28 ) = ZERO
1736                   END IF
1737                ELSE
1738                   RESULT( 27 ) = ZERO
1739                   RESULT( 28 ) = ZERO
1740                END IF
1741 *
1742 *           Call DSTEMR(V,I) to compute D1 and Z, do tests.
1743 *
1744 *           Compute D1 and Z
1745 *
1746                CALL DCOPY( N, SD, 1, D5, 1 )
1747                IF( N.GT.0 )
1748      $            CALL DCOPY( N-1, SE, 1, WORK, 1 )
1749                CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
1750 *
1751                IF( SRANGE ) THEN
1752                   NTEST = 29
1753                   IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
1754                   IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) )
1755                   IF( IU.LT.IL ) THEN
1756                      ITEMP = IU
1757                      IU = IL
1758                      IL = ITEMP
1759                   END IF
1760                   CALL DSTEMR( 'V', 'I', N, D5, WORK, VL, VU, IL, IU,
1761      $                         M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC,
1762      $                         WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
1763      $                         LIWORK-2*N, IINFO )
1764                   IF( IINFO.NE.0 ) THEN
1765                      WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,I)', IINFO,
1766      $                  N, JTYPE, IOLDSD
1767                      INFO = ABS( IINFO )
1768                      IF( IINFO.LT.0 ) THEN
1769                         RETURN
1770                      ELSE
1771                         RESULT( 29 ) = ULPINV
1772                         GO TO 280
1773                      END IF
1774                   END IF
1775 *
1776 *           Do Tests 29 and 30
1777 *
1778                   CALL DSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
1779      $                         M, RESULT( 29 ) )
1780 *
1781 *           Call DSTEMR to compute D2, do tests.
1782 *
1783 *           Compute D2
1784 *
1785                   CALL DCOPY( N, SD, 1, D5, 1 )
1786                   IF( N.GT.0 )
1787      $               CALL DCOPY( N-1, SE, 1, WORK, 1 )
1788 *
1789                   NTEST = 31
1790                   CALL DSTEMR( 'N', 'I', N, D5, WORK, VL, VU, IL, IU,
1791      $                         M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC,
1792      $                         WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
1793      $                         LIWORK-2*N, IINFO )
1794                   IF( IINFO.NE.0 ) THEN
1795                      WRITE( NOUNIT, FMT = 9999 )'DSTEMR(N,I)', IINFO,
1796      $                  N, JTYPE, IOLDSD
1797                      INFO = ABS( IINFO )
1798                      IF( IINFO.LT.0 ) THEN
1799                         RETURN
1800                      ELSE
1801                         RESULT( 31 ) = ULPINV
1802                         GO TO 280
1803                      END IF
1804                   END IF
1805 *
1806 *           Do Test 31
1807 *
1808                   TEMP1 = ZERO
1809                   TEMP2 = ZERO
1810 *
1811                   DO 240 J = 1, IU - IL + 1
1812                      TEMP1 = MAX( TEMP1, ABS( D1( J ) ),
1813      $                       ABS( D2( J ) ) )
1814                      TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
1815   240             CONTINUE
1816 *
1817                   RESULT( 31 ) = TEMP2 / MAX( UNFL,
1818      $                           ULP*MAX( TEMP1, TEMP2 ) )
1819 *
1820 *
1821 *           Call DSTEMR(V,V) to compute D1 and Z, do tests.
1822 *
1823 *           Compute D1 and Z
1824 *
1825                   CALL DCOPY( N, SD, 1, D5, 1 )
1826                   IF( N.GT.0 )
1827      $               CALL DCOPY( N-1, SE, 1, WORK, 1 )
1828                   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU )
1829 *
1830                   NTEST = 32
1831 *
1832                   IF( N.GT.0 ) THEN
1833                      IF( IL.NE.1 ) THEN
1834                         VL = D2( IL ) - MAX( HALF*
1835      $                       ( D2( IL )-D2( IL-1 ) ), ULP*ANORM,
1836      $                       TWO*RTUNFL )
1837                      ELSE
1838                         VL = D2( 1 ) - MAX( HALF*( D2( N )-D2( 1 ) ),
1839      $                       ULP*ANORM, TWO*RTUNFL )
1840                      END IF
1841                      IF( IU.NE.N ) THEN
1842                         VU = D2( IU ) + MAX( HALF*
1843      $                       ( D2( IU+1 )-D2( IU ) ), ULP*ANORM,
1844      $                       TWO*RTUNFL )
1845                      ELSE
1846                         VU = D2( N ) + MAX( HALF*( D2( N )-D2( 1 ) ),
1847      $                       ULP*ANORM, TWO*RTUNFL )
1848                      END IF
1849                   ELSE
1850                      VL = ZERO
1851                      VU = ONE
1852                   END IF
1853 *
1854                   CALL DSTEMR( 'V', 'V', N, D5, WORK, VL, VU, IL, IU,
1855      $                         M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC,
1856      $                         WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
1857      $                         LIWORK-2*N, IINFO )
1858                   IF( IINFO.NE.0 ) THEN
1859                      WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,V)', IINFO,
1860      $                  N, JTYPE, IOLDSD
1861                      INFO = ABS( IINFO )
1862                      IF( IINFO.LT.0 ) THEN
1863                         RETURN
1864                      ELSE
1865                         RESULT( 32 ) = ULPINV
1866                         GO TO 280
1867                      END IF
1868                   END IF
1869 *
1870 *           Do Tests 32 and 33
1871 *
1872                   CALL DSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK,
1873      $                         M, RESULT( 32 ) )
1874 *
1875 *           Call DSTEMR to compute D2, do tests.
1876 *
1877 *           Compute D2
1878 *
1879                   CALL DCOPY( N, SD, 1, D5, 1 )
1880                   IF( N.GT.0 )
1881      $               CALL DCOPY( N-1, SE, 1, WORK, 1 )
1882 *
1883                   NTEST = 34
1884                   CALL DSTEMR( 'N', 'V', N, D5, WORK, VL, VU, IL, IU,
1885      $                         M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC,
1886      $                         WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
1887      $                         LIWORK-2*N, IINFO )
1888                   IF( IINFO.NE.0 ) THEN
1889                      WRITE( NOUNIT, FMT = 9999 )'DSTEMR(N,V)', IINFO,
1890      $                  N, JTYPE, IOLDSD
1891                      INFO = ABS( IINFO )
1892                      IF( IINFO.LT.0 ) THEN
1893                         RETURN
1894                      ELSE
1895                         RESULT( 34 ) = ULPINV
1896                         GO TO 280
1897                      END IF
1898                   END IF
1899 *
1900 *           Do Test 34
1901 *
1902                   TEMP1 = ZERO
1903                   TEMP2 = ZERO
1904 *
1905                   DO 250 J = 1, IU - IL + 1
1906                      TEMP1 = MAX( TEMP1, ABS( D1( J ) ),
1907      $                       ABS( D2( J ) ) )
1908                      TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
1909   250             CONTINUE
1910 *
1911                   RESULT( 34 ) = TEMP2 / MAX( UNFL,
1912      $                           ULP*MAX( TEMP1, TEMP2 ) )
1913                ELSE
1914                   RESULT( 29 ) = ZERO
1915                   RESULT( 30 ) = ZERO
1916                   RESULT( 31 ) = ZERO
1917                   RESULT( 32 ) = ZERO
1918                   RESULT( 33 ) = ZERO
1919                   RESULT( 34 ) = ZERO
1920                END IF
1921 *
1922 *
1923 *           Call DSTEMR(V,A) to compute D1 and Z, do tests.
1924 *
1925 *           Compute D1 and Z
1926 *
1927                CALL DCOPY( N, SD, 1, D5, 1 )
1928                IF( N.GT.0 )
1929      $            CALL DCOPY( N-1, SE, 1, WORK, 1 )
1930 *
1931                NTEST = 35
1932 *
1933                CALL DSTEMR( 'V', 'A', N, D5, WORK, VL, VU, IL, IU,
1934      $                      M, D1, Z, LDU, N, IWORK( 1 ), TRYRAC,
1935      $                      WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
1936      $                      LIWORK-2*N, IINFO )
1937                IF( IINFO.NE.0 ) THEN
1938                   WRITE( NOUNIT, FMT = 9999 )'DSTEMR(V,A)', IINFO, N,
1939      $               JTYPE, IOLDSD
1940                   INFO = ABS( IINFO )
1941                   IF( IINFO.LT.0 ) THEN
1942                      RETURN
1943                   ELSE
1944                      RESULT( 35 ) = ULPINV
1945                      GO TO 280
1946                   END IF
1947                END IF
1948 *
1949 *           Do Tests 35 and 36
1950 *
1951                CALL DSTT22( N, M, 0, SD, SE, D1, DUMMA, Z, LDU, WORK, M,
1952      $                      RESULT( 35 ) )
1953 *
1954 *           Call DSTEMR to compute D2, do tests.
1955 *
1956 *           Compute D2
1957 *
1958                CALL DCOPY( N, SD, 1, D5, 1 )
1959                IF( N.GT.0 )
1960      $            CALL DCOPY( N-1, SE, 1, WORK, 1 )
1961 *
1962                NTEST = 37
1963                CALL DSTEMR( 'N', 'A', N, D5, WORK, VL, VU, IL, IU,
1964      $                      M, D2, Z, LDU, N, IWORK( 1 ), TRYRAC,
1965      $                      WORK( N+1 ), LWORK-N, IWORK( 2*N+1 ),
1966      $                      LIWORK-2*N, IINFO )
1967                IF( IINFO.NE.0 ) THEN
1968                   WRITE( NOUNIT, FMT = 9999 )'DSTEMR(N,A)', IINFO, N,
1969      $               JTYPE, IOLDSD
1970                   INFO = ABS( IINFO )
1971                   IF( IINFO.LT.0 ) THEN
1972                      RETURN
1973                   ELSE
1974                      RESULT( 37 ) = ULPINV
1975                      GO TO 280
1976                   END IF
1977                END IF
1978 *
1979 *           Do Test 34
1980 *
1981                TEMP1 = ZERO
1982                TEMP2 = ZERO
1983 *
1984                DO 260 J = 1, N
1985                   TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D2( J ) ) )
1986                   TEMP2 = MAX( TEMP2, ABS( D1( J )-D2( J ) ) )
1987   260          CONTINUE
1988 *
1989                RESULT( 37 ) = TEMP2 / MAX( UNFL,
1990      $                        ULP*MAX( TEMP1, TEMP2 ) )
1991             END IF
1992   270       CONTINUE
1993   280       CONTINUE
1994             NTESTT = NTESTT + NTEST
1995 *
1996 *           End of Loop -- Check for RESULT(j) > THRESH
1997 *
1998 *
1999 *           Print out tests which fail.
2000 *
2001             DO 290 JR = 1, NTEST
2002                IF( RESULT( JR ).GE.THRESH ) THEN
2003 *
2004 *                 If this is the first test to fail,
2005 *                 print a header to the data file.
2006 *
2007                   IF( NERRS.EQ.0 ) THEN
2008                      WRITE( NOUNIT, FMT = 9998 )'DST'
2009                      WRITE( NOUNIT, FMT = 9997 )
2010                      WRITE( NOUNIT, FMT = 9996 )
2011                      WRITE( NOUNIT, FMT = 9995 )'Symmetric'
2012                      WRITE( NOUNIT, FMT = 9994 )
2013 *
2014 *                    Tests performed
2015 *
2016                      WRITE( NOUNIT, FMT = 9988 )
2017                   END IF
2018                   NERRS = NERRS + 1
2019                   WRITE( NOUNIT, FMT = 9990 )N, IOLDSD, JTYPE, JR,
2020      $               RESULT( JR )
2021                END IF
2022   290       CONTINUE
2023   300    CONTINUE
2024   310 CONTINUE
2025 *
2026 *     Summary
2027 *
2028       CALL DLASUM( 'DST', NOUNIT, NERRS, NTESTT )
2029       RETURN
2030 *
2031  9999 FORMAT( ' DCHKST2STG: ', A, ' returned INFO=', I6, '.', / 9X,
2032      $  'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
2033 *
2034  9998 FORMAT( / 1X, A3, ' -- Real Symmetric eigenvalue problem' )
2035  9997 FORMAT( ' Matrix types (see DCHKST2STG for details): ' )
2036 *
2037  9996 FORMAT( / ' Special Matrices:',
2038      $      / '  1=Zero matrix.                        ',
2039      $      '  5=Diagonal: clustered entries.',
2040      $      / '  2=Identity matrix.                    ',
2041      $      '  6=Diagonal: large, evenly spaced.',
2042      $      / '  3=Diagonal: evenly spaced entries.    ',
2043      $      '  7=Diagonal: small, evenly spaced.',
2044      $      / '  4=Diagonal: geometr. spaced entries.' )
2045  9995 FORMAT( ' Dense ', A, ' Matrices:',
2046      $      / '  8=Evenly spaced eigenvals.            ',
2047      $      ' 12=Small, evenly spaced eigenvals.',
2048      $      / '  9=Geometrically spaced eigenvals.     ',
2049      $      ' 13=Matrix with random O(1) entries.',
2050      $      / ' 10=Clustered eigenvalues.              ',
2051      $      ' 14=Matrix with large random entries.',
2052      $      / ' 11=Large, evenly spaced eigenvals.     ',
2053      $      ' 15=Matrix with small random entries.' )
2054  9994 FORMAT( ' 16=Positive definite, evenly spaced eigenvalues',
2055      $      / ' 17=Positive definite, geometrically spaced eigenvlaues',
2056      $      / ' 18=Positive definite, clustered eigenvalues',
2057      $      / ' 19=Positive definite, small evenly spaced eigenvalues',
2058      $      / ' 20=Positive definite, large evenly spaced eigenvalues',
2059      $      / ' 21=Diagonally dominant tridiagonal, geometrically',
2060      $      ' spaced eigenvalues' )
2061 *
2062  9990 FORMAT( ' N=', I5, ', seed=', 4( I4, ',' ), ' type ', I2,
2063      $      ', test(', I2, ')=', G10.3 )
2064 *
2065  9988 FORMAT( / 'Test performed:  see DCHKST2STG for details.', / )
2066 *     End of DCHKST2STG
2067 *
2068       END