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