LAPACK 3.3.0

dchkst.f

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