LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dchkst2stg()

subroutine dchkst2stg ( integer  NSIZES,
integer, dimension( * )  NN,
integer  NTYPES,
logical, dimension( * )  DOTYPE,
integer, dimension( 4 )  ISEED,
double precision  THRESH,
integer  NOUNIT,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  AP,
double precision, dimension( * )  SD,
double precision, dimension( * )  SE,
double precision, dimension( * )  D1,
double precision, dimension( * )  D2,
double precision, dimension( * )  D3,
double precision, dimension( * )  D4,
double precision, dimension( * )  D5,
double precision, dimension( * )  WA1,
double precision, dimension( * )  WA2,
double precision, dimension( * )  WA3,
double precision, dimension( * )  WR,
double precision, dimension( ldu, * )  U,
integer  LDU,
double precision, dimension( ldu, * )  V,
double precision, dimension( * )  VP,
double precision, dimension( * )  TAU,
double precision, dimension( ldu, * )  Z,
double precision, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
double precision, dimension( * )  RESULT,
integer  INFO 
)

DCHKST2STG

Purpose:
 DCHKST2STG  checks the symmetric eigenvalue problem routines
 using the 2-stage reduction techniques. Since the generation
 of Q or the vectors is not available in this release, we only 
 compare the eigenvalue resulting when using the 2-stage to the 
 one considered as reference using the standard 1-stage reduction
 DSYTRD. For that, we call the standard DSYTRD and compute D1 using 
 DSTEQR, then we call the 2-stage DSYTRD_2STAGE with Upper and Lower
 and we compute D2 and D3 using DSTEQR and then we replaced tests
 3 and 4 by tests 11 and 12. test 1 and 2 remain to verify that 
 the 1-stage results are OK and can be trusted.
 This testing routine will converge to the DCHKST in the next 
 release when vectors and generation of Q will be implemented.

    DSYTRD factors A as  U S U' , where ' means transpose,
    S is symmetric tridiagonal, and U is orthogonal.
    DSYTRD can use either just the lower or just the upper triangle
    of A; DCHKST2STG checks both cases.
    U is represented as a product of Householder
    transformations, whose vectors are stored in the first
    n-1 columns of V, and whose scale factors are in TAU.

    DSPTRD does the same as DSYTRD, except that A and V are stored
    in "packed" format.

    DORGTR constructs the matrix U from the contents of V and TAU.

    DOPGTR constructs the matrix U from the contents of VP and TAU.

    DSTEQR factors S as  Z D1 Z' , where Z is the orthogonal
    matrix of eigenvectors and D1 is a diagonal matrix with
    the eigenvalues on the diagonal.  D2 is the matrix of
    eigenvalues computed when Z is not computed.

    DSTERF computes D3, the matrix of eigenvalues, by the
    PWK method, which does not yield eigenvectors.

    DPTEQR factors S as  Z4 D4 Z4' , for a
    symmetric positive definite tridiagonal matrix.
    D5 is the matrix of eigenvalues computed when Z is not
    computed.

    DSTEBZ computes selected eigenvalues.  WA1, WA2, and
    WA3 will denote eigenvalues computed to high
    absolute accuracy, with different range options.
    WR will denote eigenvalues computed to high relative
    accuracy.

    DSTEIN computes Y, the eigenvectors of S, given the
    eigenvalues.

    DSTEDC factors S as Z D1 Z' , where Z is the orthogonal
    matrix of eigenvectors and D1 is a diagonal matrix with
    the eigenvalues on the diagonal ('I' option). It may also
    update an input orthogonal matrix, usually the output
    from DSYTRD/DORGTR or DSPTRD/DOPGTR ('V' option). It may
    also just compute eigenvalues ('N' option).

    DSTEMR factors S as Z D1 Z' , where Z is the orthogonal
    matrix of eigenvectors and D1 is a diagonal matrix with
    the eigenvalues on the diagonal ('I' option).  DSTEMR
    uses the Relatively Robust Representation whenever possible.

 When DCHKST2STG is called, a number of matrix "sizes" ("n's") and a
 number of matrix "types" are specified.  For each size ("n")
 and each type of matrix, one matrix will be generated and used
 to test the symmetric eigenroutines.  For each matrix, a number
 of tests will be performed:

 (1)     | A - V S V' | / ( |A| n ulp ) DSYTRD( UPLO='U', ... )

 (2)     | I - UV' | / ( n ulp )        DORGTR( UPLO='U', ... )

 (3)     | A - V S V' | / ( |A| n ulp ) DSYTRD( UPLO='L', ... )
         replaced by | D1 - D2 | / ( |D1| ulp ) where D1 is the 
         eigenvalue matrix computed using S and D2 is the 
         eigenvalue matrix computed using S_2stage the output of
         DSYTRD_2STAGE("N", "U",....). D1 and D2 are computed 
         via DSTEQR('N',...)  

 (4)     | I - UV' | / ( n ulp )        DORGTR( UPLO='L', ... )
         replaced by | D1 - D3 | / ( |D1| ulp ) where D1 is the 
         eigenvalue matrix computed using S and D3 is the 
         eigenvalue matrix computed using S_2stage the output of
         DSYTRD_2STAGE("N", "L",....). D1 and D3 are computed 
         via DSTEQR('N',...)  

 (5-8)   Same as 1-4, but for DSPTRD and DOPGTR.

 (9)     | S - Z D Z' | / ( |S| n ulp ) DSTEQR('V',...)

 (10)    | I - ZZ' | / ( n ulp )        DSTEQR('V',...)

 (11)    | D1 - D2 | / ( |D1| ulp )        DSTEQR('N',...)

 (12)    | D1 - D3 | / ( |D1| ulp )        DSTERF

 (13)    0 if the true eigenvalues (computed by sturm count)
         of S are within THRESH of
         those in D1.  2*THRESH if they are not.  (Tested using
         DSTECH)

 For S positive definite,

 (14)    | S - Z4 D4 Z4' | / ( |S| n ulp ) DPTEQR('V',...)

 (15)    | I - Z4 Z4' | / ( n ulp )        DPTEQR('V',...)

 (16)    | D4 - D5 | / ( 100 |D4| ulp )       DPTEQR('N',...)

 When S is also diagonally dominant by the factor gamma < 1,

 (17)    max | D4(i) - WR(i) | / ( |D4(i)| omega ) ,
          i
         omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
                                              DSTEBZ( 'A', 'E', ...)

 (18)    | WA1 - D3 | / ( |D3| ulp )          DSTEBZ( 'A', 'E', ...)

 (19)    ( max { min | WA2(i)-WA3(j) | } +
            i     j
           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
            i     j
                                              DSTEBZ( 'I', 'E', ...)

 (20)    | S - Y WA1 Y' | / ( |S| n ulp )  DSTEBZ, SSTEIN

 (21)    | I - Y Y' | / ( n ulp )          DSTEBZ, SSTEIN

 (22)    | S - Z D Z' | / ( |S| n ulp )    DSTEDC('I')

 (23)    | I - ZZ' | / ( n ulp )           DSTEDC('I')

 (24)    | S - Z D Z' | / ( |S| n ulp )    DSTEDC('V')

 (25)    | I - ZZ' | / ( n ulp )           DSTEDC('V')

 (26)    | D1 - D2 | / ( |D1| ulp )           DSTEDC('V') and
                                              DSTEDC('N')

 Test 27 is disabled at the moment because DSTEMR does not
 guarantee high relatvie accuracy.

 (27)    max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
          i
         omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
                                              DSTEMR('V', 'A')

 (28)    max | D6(i) - WR(i) | / ( |D6(i)| omega ) ,
          i
         omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4
                                              DSTEMR('V', 'I')

 Tests 29 through 34 are disable at present because DSTEMR
 does not handle partial spectrum requests.

 (29)    | S - Z D Z' | / ( |S| n ulp )    DSTEMR('V', 'I')

 (30)    | I - ZZ' | / ( n ulp )           DSTEMR('V', 'I')

 (31)    ( max { min | WA2(i)-WA3(j) | } +
            i     j
           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
            i     j
         DSTEMR('N', 'I') vs. SSTEMR('V', 'I')

 (32)    | S - Z D Z' | / ( |S| n ulp )    DSTEMR('V', 'V')

 (33)    | I - ZZ' | / ( n ulp )           DSTEMR('V', 'V')

 (34)    ( max { min | WA2(i)-WA3(j) | } +
            i     j
           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
            i     j
         DSTEMR('N', 'V') vs. SSTEMR('V', 'V')

 (35)    | S - Z D Z' | / ( |S| n ulp )    DSTEMR('V', 'A')

 (36)    | I - ZZ' | / ( n ulp )           DSTEMR('V', 'A')

 (37)    ( max { min | WA2(i)-WA3(j) | } +
            i     j
           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
            i     j
         DSTEMR('N', 'A') vs. SSTEMR('V', 'A')

 The "sizes" are specified by an array NN(1:NSIZES); the value of
 each element NN(j) specifies one size.
 The "types" are specified by a logical array DOTYPE( 1:NTYPES );
 if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
 Currently, the list of possible types is:

 (1)  The zero matrix.
 (2)  The identity matrix.

 (3)  A diagonal matrix with evenly spaced entries
      1, ..., ULP  and random signs.
      (ULP = (first number larger than 1) - 1 )
 (4)  A diagonal matrix with geometrically spaced entries
      1, ..., ULP  and random signs.
 (5)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
      and random signs.

 (6)  Same as (4), but multiplied by SQRT( overflow threshold )
 (7)  Same as (4), but multiplied by SQRT( underflow threshold )

 (8)  A matrix of the form  U' D U, where U is orthogonal and
      D has evenly spaced entries 1, ..., ULP with random signs
      on the diagonal.

 (9)  A matrix of the form  U' D U, where U is orthogonal and
      D has geometrically spaced entries 1, ..., ULP with random
      signs on the diagonal.

 (10) A matrix of the form  U' D U, where U is orthogonal and
      D has "clustered" entries 1, ULP,..., ULP with random
      signs on the diagonal.

 (11) Same as (8), but multiplied by SQRT( overflow threshold )
 (12) Same as (8), but multiplied by SQRT( underflow threshold )

 (13) Symmetric matrix with random entries chosen from (-1,1).
 (14) Same as (13), but multiplied by SQRT( overflow threshold )
 (15) Same as (13), but multiplied by SQRT( underflow threshold )
 (16) Same as (8), but diagonal elements are all positive.
 (17) Same as (9), but diagonal elements are all positive.
 (18) Same as (10), but diagonal elements are all positive.
 (19) Same as (16), but multiplied by SQRT( overflow threshold )
 (20) Same as (16), but multiplied by SQRT( underflow threshold )
 (21) A diagonally dominant tridiagonal matrix with geometrically
      spaced diagonal entries 1, ..., ULP.
Parameters
[in]NSIZES
          NSIZES is INTEGER
          The number of sizes of matrices to use.  If it is zero,
          DCHKST2STG does nothing.  It must be at least zero.
[in]NN
          NN is INTEGER array, dimension (NSIZES)
          An array containing the sizes to be used for the matrices.
          Zero values will be skipped.  The values must be at least
          zero.
[in]NTYPES
          NTYPES is INTEGER
          The number of elements in DOTYPE.   If it is zero, DCHKST2STG
          does nothing.  It must be at least zero.  If it is MAXTYP+1
          and NSIZES is 1, then an additional type, MAXTYP+1 is
          defined, which is to use whatever matrix is in A.  This
          is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
          DOTYPE(MAXTYP+1) is .TRUE. .
[in]DOTYPE
          DOTYPE is LOGICAL array, dimension (NTYPES)
          If DOTYPE(j) is .TRUE., then for each size in NN a
          matrix of that size and of type j will be generated.
          If NTYPES is smaller than the maximum number of types
          defined (PARAMETER MAXTYP), then types NTYPES+1 through
          MAXTYP will not be generated.  If NTYPES is larger
          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
          will be ignored.
[in,out]ISEED
          ISEED is INTEGER array, dimension (4)
          On entry ISEED specifies the seed of the random number
          generator. The array elements should be between 0 and 4095;
          if not they will be reduced mod 4096.  Also, ISEED(4) must
          be odd.  The random number generator uses a linear
          congruential sequence limited to small integers, and so
          should produce machine independent random numbers. The
          values of ISEED are changed on exit, and can be used in the
          next call to DCHKST2STG to continue the same random number
          sequence.
[in]THRESH
          THRESH is DOUBLE PRECISION
          A test will count as "failed" if the "error", computed as
          described above, exceeds THRESH.  Note that the error
          is scaled to be O(1), so THRESH should be a reasonably
          small multiple of 1, e.g., 10 or 100.  In particular,
          it should not depend on the precision (single vs. double)
          or the size of the matrix.  It must be at least zero.
[in]NOUNIT
          NOUNIT is INTEGER
          The FORTRAN unit number for printing out error messages
          (e.g., if a routine returns IINFO not equal to 0.)
[in,out]A
          A is DOUBLE PRECISION array of
                                  dimension ( LDA , max(NN) )
          Used to hold the matrix whose eigenvalues are to be
          computed.  On exit, A contains the last matrix actually
          used.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  It must be at
          least 1 and at least max( NN ).
[out]AP
          AP is DOUBLE PRECISION array of
                      dimension( max(NN)*max(NN+1)/2 )
          The matrix A stored in packed format.
[out]SD
          SD is DOUBLE PRECISION array of
                             dimension( max(NN) )
          The diagonal of the tridiagonal matrix computed by DSYTRD.
          On exit, SD and SE contain the tridiagonal form of the
          matrix in A.
[out]SE
          SE is DOUBLE PRECISION array of
                             dimension( max(NN) )
          The off-diagonal of the tridiagonal matrix computed by
          DSYTRD.  On exit, SD and SE contain the tridiagonal form of
          the matrix in A.
[out]D1
          D1 is DOUBLE PRECISION array of
                             dimension( max(NN) )
          The eigenvalues of A, as computed by DSTEQR simlutaneously
          with Z.  On exit, the eigenvalues in D1 correspond with the
          matrix in A.
[out]D2
          D2 is DOUBLE PRECISION array of
                             dimension( max(NN) )
          The eigenvalues of A, as computed by DSTEQR if Z is not
          computed.  On exit, the eigenvalues in D2 correspond with
          the matrix in A.
[out]D3
          D3 is DOUBLE PRECISION array of
                             dimension( max(NN) )
          The eigenvalues of A, as computed by DSTERF.  On exit, the
          eigenvalues in D3 correspond with the matrix in A.
[out]D4
          D4 is DOUBLE PRECISION array of
                             dimension( max(NN) )
          The eigenvalues of A, as computed by DPTEQR(V).
          DPTEQR factors S as  Z4 D4 Z4*
          On exit, the eigenvalues in D4 correspond with the matrix in A.
[out]D5
          D5 is DOUBLE PRECISION array of
                             dimension( max(NN) )
          The eigenvalues of A, as computed by DPTEQR(N)
          when Z is not computed. On exit, the
          eigenvalues in D4 correspond with the matrix in A.
[out]WA1
          WA1 is DOUBLE PRECISION array of
                             dimension( max(NN) )
          All eigenvalues of A, computed to high
          absolute accuracy, with different range options.
          as computed by DSTEBZ.
[out]WA2
          WA2 is DOUBLE PRECISION array of
                             dimension( max(NN) )
          Selected eigenvalues of A, computed to high
          absolute accuracy, with different range options.
          as computed by DSTEBZ.
          Choose random values for IL and IU, and ask for the
          IL-th through IU-th eigenvalues.
[out]WA3
          WA3 is DOUBLE PRECISION array of
                             dimension( max(NN) )
          Selected eigenvalues of A, computed to high
          absolute accuracy, with different range options.
          as computed by DSTEBZ.
          Determine the values VL and VU of the IL-th and IU-th
          eigenvalues and ask for all eigenvalues in this range.
[out]WR
          WR is DOUBLE PRECISION array of
                             dimension( max(NN) )
          All eigenvalues of A, computed to high
          absolute accuracy, with different options.
          as computed by DSTEBZ.
[out]U
          U is DOUBLE PRECISION array of
                             dimension( LDU, max(NN) ).
          The orthogonal matrix computed by DSYTRD + DORGTR.
[in]LDU
          LDU is INTEGER
          The leading dimension of U, Z, and V.  It must be at least 1
          and at least max( NN ).
[out]V
          V is DOUBLE PRECISION array of
                             dimension( LDU, max(NN) ).
          The Housholder vectors computed by DSYTRD in reducing A to
          tridiagonal form.  The vectors computed with UPLO='U' are
          in the upper triangle, and the vectors computed with UPLO='L'
          are in the lower triangle.  (As described in DSYTRD, the
          sub- and superdiagonal are not set to 1, although the
          true Householder vector has a 1 in that position.  The
          routines that use V, such as DORGTR, set those entries to
          1 before using them, and then restore them later.)
[out]VP
          VP is DOUBLE PRECISION array of
                      dimension( max(NN)*max(NN+1)/2 )
          The matrix V stored in packed format.
[out]TAU
          TAU is DOUBLE PRECISION array of
                             dimension( max(NN) )
          The Householder factors computed by DSYTRD in reducing A
          to tridiagonal form.
[out]Z
          Z is DOUBLE PRECISION array of
                             dimension( LDU, max(NN) ).
          The orthogonal matrix of eigenvectors computed by DSTEQR,
          DPTEQR, and DSTEIN.
[out]WORK
          WORK is DOUBLE PRECISION array of
                      dimension( LWORK )
[in]LWORK
          LWORK is INTEGER
          The number of entries in WORK.  This must be at least
          1 + 4 * Nmax + 2 * Nmax * lg Nmax + 3 * Nmax**2
          where Nmax = max( NN(j), 2 ) and lg = log base 2.
[out]IWORK
          IWORK is INTEGER array,
          Workspace.
[out]LIWORK
          LIWORK is INTEGER
          The number of entries in IWORK.  This must be at least
                  6 + 6*Nmax + 5 * Nmax * lg Nmax
          where Nmax = max( NN(j), 2 ) and lg = log base 2.
[out]RESULT
          RESULT is DOUBLE PRECISION array, dimension (26)
          The values computed by the tests described above.
          The values are currently limited to 1/ulp, to avoid
          overflow.
[out]INFO
          INFO is INTEGER
          If 0, then everything ran OK.
           -1: NSIZES < 0
           -2: Some NN(j) < 0
           -3: NTYPES < 0
           -5: THRESH < 0
           -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
          -23: LDU < 1 or LDU < NMAX.
          -29: LWORK too small.
          If  DLATMR, SLATMS, DSYTRD, DORGTR, DSTEQR, SSTERF,
              or DORMC2 returns an error code, the
              absolute value of it is returned.

-----------------------------------------------------------------------

       Some Local Variables and Parameters:
       ---- ----- --------- --- ----------
       ZERO, ONE       Real 0 and 1.
       MAXTYP          The number of types defined.
       NTEST           The number of tests performed, or which can
                       be performed so far, for the current matrix.
       NTESTT          The total number of tests performed so far.
       NBLOCK          Blocksize as returned by ENVIR.
       NMAX            Largest value in NN.
       NMATS           The number of matrices generated so far.
       NERRS           The number of tests which have exceeded THRESH
                       so far.
       COND, IMODE     Values to be passed to the matrix generators.
       ANORM           Norm of A; passed to matrix generators.

       OVFL, UNFL      Overflow and underflow thresholds.
       ULP, ULPINV     Finest relative precision and its inverse.
       RTOVFL, RTUNFL  Square roots of the previous 2 values.
               The following four arrays decode JTYPE:
       KTYPE(j)        The general type (1-10) for type "j".
       KMODE(j)        The MODE value to be passed to the matrix
                       generator for type "j".
       KMAGN(j)        The order of magnitude ( O(1),
                       O(overflow^(1/2) ), O(underflow^(1/2) )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 608 of file dchkst2stg.f.

612 *
613 * -- LAPACK test routine --
614 * -- LAPACK is a software package provided by Univ. of Tennessee, --
615 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
616 *
617 * .. Scalar Arguments ..
618  INTEGER INFO, LDA, LDU, LIWORK, LWORK, NOUNIT, NSIZES,
619  $ NTYPES
620  DOUBLE PRECISION THRESH
621 * ..
622 * .. Array Arguments ..
623  LOGICAL DOTYPE( * )
624  INTEGER ISEED( 4 ), IWORK( * ), NN( * )
625  DOUBLE PRECISION A( LDA, * ), AP( * ), D1( * ), D2( * ),
626  $ D3( * ), D4( * ), D5( * ), RESULT( * ),
627  $ SD( * ), SE( * ), TAU( * ), U( LDU, * ),
628  $ V( LDU, * ), VP( * ), WA1( * ), WA2( * ),
629  $ WA3( * ), WORK( * ), WR( * ), Z( LDU, * )
630 * ..
631 *
632 * =====================================================================
633 *
634 * .. Parameters ..
635  DOUBLE PRECISION ZERO, ONE, TWO, EIGHT, TEN, HUN
636  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
637  $ eight = 8.0d0, ten = 10.0d0, hun = 100.0d0 )
638  DOUBLE PRECISION HALF
639  parameter( half = one / two )
640  INTEGER MAXTYP
641  parameter( maxtyp = 21 )
642  LOGICAL SRANGE
643  parameter( srange = .false. )
644  LOGICAL SREL
645  parameter( srel = .false. )
646 * ..
647 * .. Local Scalars ..
648  LOGICAL BADNN, TRYRAC
649  INTEGER I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, J, JC,
650  $ JR, JSIZE, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC,
651  $ M, M2, M3, MTYPES, N, NAP, NBLOCK, NERRS,
652  $ NMATS, NMAX, NSPLIT, NTEST, NTESTT, LH, LW
653  DOUBLE PRECISION ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL,
654  $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP,
655  $ ULPINV, UNFL, VL, VU
656 * ..
657 * .. Local Arrays ..
658  INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
659  $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
660  $ KTYPE( MAXTYP )
661  DOUBLE PRECISION DUMMA( 1 )
662 * ..
663 * .. External Functions ..
664  INTEGER ILAENV
665  DOUBLE PRECISION DLAMCH, DLARND, DSXT1
666  EXTERNAL ilaenv, dlamch, dlarnd, dsxt1
667 * ..
668 * .. External Subroutines ..
669  EXTERNAL dcopy, dlabad, dlacpy, dlaset, dlasum, dlatmr,
673  $ dsytrd_2stage
674 * ..
675 * .. Intrinsic Functions ..
676  INTRINSIC abs, dble, int, log, max, min, sqrt
677 * ..
678 * .. Data statements ..
679  DATA ktype / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8,
680  $ 8, 8, 9, 9, 9, 9, 9, 10 /
681  DATA kmagn / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
682  $ 2, 3, 1, 1, 1, 2, 3, 1 /
683  DATA kmode / 0, 0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
684  $ 0, 0, 4, 3, 1, 4, 4, 3 /
685 * ..
686 * .. Executable Statements ..
687 *
688 * Keep ftnchek happy
689  idumma( 1 ) = 1
690 *
691 * Check for errors
692 *
693  ntestt = 0
694  info = 0
695 *
696 * Important constants
697 *
698  badnn = .false.
699  tryrac = .true.
700  nmax = 1
701  DO 10 j = 1, nsizes
702  nmax = max( nmax, nn( j ) )
703  IF( nn( j ).LT.0 )
704  $ badnn = .true.
705  10 CONTINUE
706 *
707  nblock = ilaenv( 1, 'DSYTRD', 'L', nmax, -1, -1, -1 )
708  nblock = min( nmax, max( 1, nblock ) )
709 *
710 * Check for errors
711 *
712  IF( nsizes.LT.0 ) THEN
713  info = -1
714  ELSE IF( badnn ) THEN
715  info = -2
716  ELSE IF( ntypes.LT.0 ) THEN
717  info = -3
718  ELSE IF( lda.LT.nmax ) THEN
719  info = -9
720  ELSE IF( ldu.LT.nmax ) THEN
721  info = -23
722  ELSE IF( 2*max( 2, nmax )**2.GT.lwork ) THEN
723  info = -29
724  END IF
725 *
726  IF( info.NE.0 ) THEN
727  CALL xerbla( 'DCHKST2STG', -info )
728  RETURN
729  END IF
730 *
731 * Quick return if possible
732 *
733  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
734  $ RETURN
735 *
736 * More Important constants
737 *
738  unfl = dlamch( 'Safe minimum' )
739  ovfl = one / unfl
740  CALL dlabad( unfl, ovfl )
741  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
742  ulpinv = one / ulp
743  log2ui = int( log( ulpinv ) / log( two ) )
744  rtunfl = sqrt( unfl )
745  rtovfl = sqrt( ovfl )
746 *
747 * Loop over sizes, types
748 *
749  DO 20 i = 1, 4
750  iseed2( i ) = iseed( i )
751  20 CONTINUE
752  nerrs = 0
753  nmats = 0
754 *
755  DO 310 jsize = 1, nsizes
756  n = nn( jsize )
757  IF( n.GT.0 ) THEN
758  lgn = int( log( dble( n ) ) / log( two ) )
759  IF( 2**lgn.LT.n )
760  $ lgn = lgn + 1
761  IF( 2**lgn.LT.n )
762  $ lgn = lgn + 1
763  lwedc = 1 + 4*n + 2*n*lgn + 4*n**2
764  liwedc = 6 + 6*n + 5*n*lgn
765  ELSE
766  lwedc = 8
767  liwedc = 12
768  END IF
769  nap = ( n*( n+1 ) ) / 2
770  aninv = one / dble( max( 1, n ) )
771 *
772  IF( nsizes.NE.1 ) THEN
773  mtypes = min( maxtyp, ntypes )
774  ELSE
775  mtypes = min( maxtyp+1, ntypes )
776  END IF
777 *
778  DO 300 jtype = 1, mtypes
779  IF( .NOT.dotype( jtype ) )
780  $ GO TO 300
781  nmats = nmats + 1
782  ntest = 0
783 *
784  DO 30 j = 1, 4
785  ioldsd( j ) = iseed( j )
786  30 CONTINUE
787 *
788 * Compute "A"
789 *
790 * Control parameters:
791 *
792 * KMAGN KMODE KTYPE
793 * =1 O(1) clustered 1 zero
794 * =2 large clustered 2 identity
795 * =3 small exponential (none)
796 * =4 arithmetic diagonal, (w/ eigenvalues)
797 * =5 random log symmetric, w/ eigenvalues
798 * =6 random (none)
799 * =7 random diagonal
800 * =8 random symmetric
801 * =9 positive definite
802 * =10 diagonally dominant tridiagonal
803 *
804  IF( mtypes.GT.maxtyp )
805  $ GO TO 100
806 *
807  itype = ktype( jtype )
808  imode = kmode( jtype )
809 *
810 * Compute norm
811 *
812  GO TO ( 40, 50, 60 )kmagn( jtype )
813 *
814  40 CONTINUE
815  anorm = one
816  GO TO 70
817 *
818  50 CONTINUE
819  anorm = ( rtovfl*ulp )*aninv
820  GO TO 70
821 *
822  60 CONTINUE
823  anorm = rtunfl*n*ulpinv
824  GO TO 70
825 *
826  70 CONTINUE
827 *
828  CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
829  iinfo = 0
830  IF( jtype.LE.15 ) THEN
831  cond = ulpinv
832  ELSE
833  cond = ulpinv*aninv / ten
834  END IF
835 *
836 * Special Matrices -- Identity & Jordan block
837 *
838 * Zero
839 *
840  IF( itype.EQ.1 ) THEN
841  iinfo = 0
842 *
843  ELSE IF( itype.EQ.2 ) THEN
844 *
845 * Identity
846 *
847  DO 80 jc = 1, n
848  a( jc, jc ) = anorm
849  80 CONTINUE
850 *
851  ELSE IF( itype.EQ.4 ) THEN
852 *
853 * Diagonal Matrix, [Eigen]values Specified
854 *
855  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
856  $ anorm, 0, 0, 'N', a, lda, work( n+1 ),
857  $ iinfo )
858 *
859 *
860  ELSE IF( itype.EQ.5 ) THEN
861 *
862 * Symmetric, eigenvalues specified
863 *
864  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
865  $ anorm, n, n, 'N', a, lda, work( n+1 ),
866  $ iinfo )
867 *
868  ELSE IF( itype.EQ.7 ) THEN
869 *
870 * Diagonal, random eigenvalues
871 *
872  CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
873  $ 'T', 'N', work( n+1 ), 1, one,
874  $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
875  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
876 *
877  ELSE IF( itype.EQ.8 ) THEN
878 *
879 * Symmetric, random eigenvalues
880 *
881  CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
882  $ 'T', 'N', work( n+1 ), 1, one,
883  $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
884  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
885 *
886  ELSE IF( itype.EQ.9 ) THEN
887 *
888 * Positive definite, eigenvalues specified.
889 *
890  CALL dlatms( n, n, 'S', iseed, 'P', work, imode, cond,
891  $ anorm, n, n, 'N', a, lda, work( n+1 ),
892  $ iinfo )
893 *
894  ELSE IF( itype.EQ.10 ) THEN
895 *
896 * Positive definite tridiagonal, eigenvalues specified.
897 *
898  CALL dlatms( n, n, 'S', iseed, 'P', work, imode, cond,
899  $ anorm, 1, 1, 'N', a, lda, work( n+1 ),
900  $ iinfo )
901  DO 90 i = 2, n
902  temp1 = abs( a( i-1, i ) ) /
903  $ sqrt( abs( a( i-1, i-1 )*a( i, i ) ) )
904  IF( temp1.GT.half ) THEN
905  a( i-1, i ) = half*sqrt( abs( a( i-1, i-1 )*a( i,
906  $ i ) ) )
907  a( i, i-1 ) = a( i-1, i )
908  END IF
909  90 CONTINUE
910 *
911  ELSE
912 *
913  iinfo = 1
914  END IF
915 *
916  IF( iinfo.NE.0 ) THEN
917  WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
918  $ ioldsd
919  info = abs( iinfo )
920  RETURN
921  END IF
922 *
923  100 CONTINUE
924 *
925 * Call DSYTRD and DORGTR to compute S and U from
926 * upper triangle.
927 *
928  CALL dlacpy( 'U', n, n, a, lda, v, ldu )
929 *
930  ntest = 1
931  CALL dsytrd( 'U', n, v, ldu, sd, se, tau, work, lwork,
932  $ iinfo )
933 *
934  IF( iinfo.NE.0 ) THEN
935  WRITE( nounit, fmt = 9999 )'DSYTRD(U)', iinfo, n, jtype,
936  $ ioldsd
937  info = abs( iinfo )
938  IF( iinfo.LT.0 ) THEN
939  RETURN
940  ELSE
941  result( 1 ) = ulpinv
942  GO TO 280
943  END IF
944  END IF
945 *
946  CALL dlacpy( 'U', n, n, v, ldu, u, ldu )
947 *
948  ntest = 2
949  CALL dorgtr( 'U', n, u, ldu, tau, work, lwork, iinfo )
950  IF( iinfo.NE.0 ) THEN
951  WRITE( nounit, fmt = 9999 )'DORGTR(U)', iinfo, n, jtype,
952  $ ioldsd
953  info = abs( iinfo )
954  IF( iinfo.LT.0 ) THEN
955  RETURN
956  ELSE
957  result( 2 ) = ulpinv
958  GO TO 280
959  END IF
960  END IF
961 *
962 * Do tests 1 and 2
963 *
964  CALL dsyt21( 2, 'Upper', n, 1, a, lda, sd, se, u, ldu, v,
965  $ ldu, tau, work, result( 1 ) )
966  CALL dsyt21( 3, 'Upper', n, 1, a, lda, sd, se, u, ldu, v,
967  $ ldu, tau, work, result( 2 ) )
968 *
969 * Compute D1 the eigenvalues resulting from the tridiagonal
970 * form using the standard 1-stage algorithm and use it as a
971 * reference to compare with the 2-stage technique
972 *
973 * Compute D1 from the 1-stage and used as reference for the
974 * 2-stage
975 *
976  CALL dcopy( n, sd, 1, d1, 1 )
977  IF( n.GT.0 )
978  $ CALL dcopy( n-1, se, 1, work, 1 )
979 *
980  CALL dsteqr( 'N', n, d1, work, work( n+1 ), ldu,
981  $ work( n+1 ), iinfo )
982  IF( iinfo.NE.0 ) THEN
983  WRITE( nounit, fmt = 9999 )'DSTEQR(N)', iinfo, n, jtype,
984  $ ioldsd
985  info = abs( iinfo )
986  IF( iinfo.LT.0 ) THEN
987  RETURN
988  ELSE
989  result( 3 ) = ulpinv
990  GO TO 280
991  END IF
992  END IF
993 *
994 * 2-STAGE TRD Upper case is used to compute D2.
995 * Note to set SD and SE to zero to be sure not reusing
996 * the one from above. Compare it with D1 computed
997 * using the 1-stage.
998 *
999  CALL dlaset( 'Full', n, 1, zero, zero, sd, n )
1000  CALL dlaset( 'Full', n, 1, zero, zero, se, n )
1001  CALL dlacpy( "U", n, n, a, lda, v, ldu )
1002  lh = max(1, 4*n)
1003  lw = lwork - lh
1004  CALL dsytrd_2stage( 'N', "U", n, v, ldu, sd, se, tau,
1005  $ work, lh, work( lh+1 ), lw, iinfo )
1006 *
1007 * Compute D2 from the 2-stage Upper case
1008 *
1009  CALL dcopy( n, sd, 1, d2, 1 )
1010  IF( n.GT.0 )
1011  $ CALL dcopy( n-1, se, 1, work, 1 )
1012 *
1013  CALL dsteqr( 'N', n, d2, work, work( n+1 ), ldu,
1014  $ work( n+1 ), iinfo )
1015  IF( iinfo.NE.0 ) THEN
1016  WRITE( nounit, fmt = 9999 )'DSTEQR(N)', iinfo, n, jtype,
1017  $ ioldsd
1018  info = abs( iinfo )
1019  IF( iinfo.LT.0 ) THEN
1020  RETURN
1021  ELSE
1022  result( 3 ) = ulpinv
1023  GO TO 280
1024  END IF
1025  END IF
1026 *
1027 * 2-STAGE TRD Lower case is used to compute D3.
1028 * Note to set SD and SE to zero to be sure not reusing
1029 * the one from above. Compare it with D1 computed
1030 * using the 1-stage.
1031 *
1032  CALL dlaset( 'Full', n, 1, zero, zero, sd, n )
1033  CALL dlaset( 'Full', n, 1, zero, zero, se, n )
1034  CALL dlacpy( "L", n, n, a, lda, v, ldu )
1035  CALL dsytrd_2stage( 'N', "L", n, v, ldu, sd, se, tau,
1036  $ work, lh, work( lh+1 ), lw, iinfo )
1037 *
1038 * Compute D3 from the 2-stage Upper case
1039 *
1040  CALL dcopy( n, sd, 1, d3, 1 )
1041  IF( n.GT.0 )
1042  $ CALL dcopy( n-1, se, 1, work, 1 )
1043 *
1044  CALL dsteqr( 'N', n, d3, work, work( n+1 ), ldu,
1045  $ work( n+1 ), iinfo )
1046  IF( iinfo.NE.0 ) THEN
1047  WRITE( nounit, fmt = 9999 )'DSTEQR(N)', iinfo, n, jtype,
1048  $ ioldsd
1049  info = abs( iinfo )
1050  IF( iinfo.LT.0 ) THEN
1051  RETURN
1052  ELSE
1053  result( 4 ) = ulpinv
1054  GO TO 280
1055  END IF
1056  END IF
1057 *
1058 * Do Tests 3 and 4 which are similar to 11 and 12 but with the
1059 * D1 computed using the standard 1-stage reduction as reference
1060 *
1061  ntest = 4
1062  temp1 = zero
1063  temp2 = zero
1064  temp3 = zero
1065  temp4 = zero
1066 *
1067  DO 151 j = 1, n
1068  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
1069  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1070  temp3 = max( temp3, abs( d1( j ) ), abs( d3( j ) ) )
1071  temp4 = max( temp4, abs( d1( j )-d3( j ) ) )
1072  151 CONTINUE
1073 *
1074  result( 3 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
1075  result( 4 ) = temp4 / max( unfl, ulp*max( temp3, temp4 ) )
1076 *
1077 * Store the upper triangle of A in AP
1078 *
1079  i = 0
1080  DO 120 jc = 1, n
1081  DO 110 jr = 1, jc
1082  i = i + 1
1083  ap( i ) = a( jr, jc )
1084  110 CONTINUE
1085  120 CONTINUE
1086 *
1087 * Call DSPTRD and DOPGTR to compute S and U from AP
1088 *
1089  CALL dcopy( nap, ap, 1, vp, 1 )
1090 *
1091  ntest = 5
1092  CALL dsptrd( 'U', n, vp, sd, se, tau, iinfo )
1093 *
1094  IF( iinfo.NE.0 ) THEN
1095  WRITE( nounit, fmt = 9999 )'DSPTRD(U)', iinfo, n, jtype,
1096  $ ioldsd
1097  info = abs( iinfo )
1098  IF( iinfo.LT.0 ) THEN
1099  RETURN
1100  ELSE
1101  result( 5 ) = ulpinv
1102  GO TO 280
1103  END IF
1104  END IF
1105 *
1106  ntest = 6
1107  CALL dopgtr( 'U', n, vp, tau, u, ldu, work, iinfo )
1108  IF( iinfo.NE.0 ) THEN
1109  WRITE( nounit, fmt = 9999 )'DOPGTR(U)', iinfo, n, jtype,
1110  $ ioldsd
1111  info = abs( iinfo )
1112  IF( iinfo.LT.0 ) THEN
1113  RETURN
1114  ELSE
1115  result( 6 ) = ulpinv
1116  GO TO 280
1117  END IF
1118  END IF
1119 *
1120 * Do tests 5 and 6
1121 *
1122  CALL dspt21( 2, 'Upper', n, 1, ap, sd, se, u, ldu, vp, tau,
1123  $ work, result( 5 ) )
1124  CALL dspt21( 3, 'Upper', n, 1, ap, sd, se, u, ldu, vp, tau,
1125  $ work, result( 6 ) )
1126 *
1127 * Store the lower triangle of A in AP
1128 *
1129  i = 0
1130  DO 140 jc = 1, n
1131  DO 130 jr = jc, n
1132  i = i + 1
1133  ap( i ) = a( jr, jc )
1134  130 CONTINUE
1135  140 CONTINUE
1136 *
1137 * Call DSPTRD and DOPGTR to compute S and U from AP
1138 *
1139  CALL dcopy( nap, ap, 1, vp, 1 )
1140 *
1141  ntest = 7
1142  CALL dsptrd( 'L', n, vp, sd, se, tau, iinfo )
1143 *
1144  IF( iinfo.NE.0 ) THEN
1145  WRITE( nounit, fmt = 9999 )'DSPTRD(L)', iinfo, n, jtype,
1146  $ ioldsd
1147  info = abs( iinfo )
1148  IF( iinfo.LT.0 ) THEN
1149  RETURN
1150  ELSE
1151  result( 7 ) = ulpinv
1152  GO TO 280
1153  END IF
1154  END IF
1155 *
1156  ntest = 8
1157  CALL dopgtr( 'L', n, vp, tau, u, ldu, work, iinfo )
1158  IF( iinfo.NE.0 ) THEN
1159  WRITE( nounit, fmt = 9999 )'DOPGTR(L)', iinfo, n, jtype,
1160  $ ioldsd
1161  info = abs( iinfo )
1162  IF( iinfo.LT.0 ) THEN
1163  RETURN
1164  ELSE
1165  result( 8 ) = ulpinv
1166  GO TO 280
1167  END IF
1168  END IF
1169 *
1170  CALL dspt21( 2, 'Lower', n, 1, ap, sd, se, u, ldu, vp, tau,
1171  $ work, result( 7 ) )
1172  CALL dspt21( 3, 'Lower', n, 1, ap, sd, se, u, ldu, vp, tau,
1173  $ work, result( 8 ) )
1174 *
1175 * Call DSTEQR to compute D1, D2, and Z, do tests.
1176 *
1177 * Compute D1 and Z
1178 *
1179  CALL dcopy( n, sd, 1, d1, 1 )
1180  IF( n.GT.0 )
1181  $ CALL dcopy( n-1, se, 1, work, 1 )
1182  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1183 *
1184  ntest = 9
1185  CALL dsteqr( 'V', n, d1, work, z, ldu, work( n+1 ), iinfo )
1186  IF( iinfo.NE.0 ) THEN
1187  WRITE( nounit, fmt = 9999 )'DSTEQR(V)', iinfo, n, jtype,
1188  $ ioldsd
1189  info = abs( iinfo )
1190  IF( iinfo.LT.0 ) THEN
1191  RETURN
1192  ELSE
1193  result( 9 ) = ulpinv
1194  GO TO 280
1195  END IF
1196  END IF
1197 *
1198 * Compute D2
1199 *
1200  CALL dcopy( n, sd, 1, d2, 1 )
1201  IF( n.GT.0 )
1202  $ CALL dcopy( n-1, se, 1, work, 1 )
1203 *
1204  ntest = 11
1205  CALL dsteqr( 'N', n, d2, work, work( n+1 ), ldu,
1206  $ work( n+1 ), iinfo )
1207  IF( iinfo.NE.0 ) THEN
1208  WRITE( nounit, fmt = 9999 )'DSTEQR(N)', iinfo, n, jtype,
1209  $ ioldsd
1210  info = abs( iinfo )
1211  IF( iinfo.LT.0 ) THEN
1212  RETURN
1213  ELSE
1214  result( 11 ) = ulpinv
1215  GO TO 280
1216  END IF
1217  END IF
1218 *
1219 * Compute D3 (using PWK method)
1220 *
1221  CALL dcopy( n, sd, 1, d3, 1 )
1222  IF( n.GT.0 )
1223  $ CALL dcopy( n-1, se, 1, work, 1 )
1224 *
1225  ntest = 12
1226  CALL dsterf( n, d3, work, iinfo )
1227  IF( iinfo.NE.0 ) THEN
1228  WRITE( nounit, fmt = 9999 )'DSTERF', iinfo, n, jtype,
1229  $ ioldsd
1230  info = abs( iinfo )
1231  IF( iinfo.LT.0 ) THEN
1232  RETURN
1233  ELSE
1234  result( 12 ) = ulpinv
1235  GO TO 280
1236  END IF
1237  END IF
1238 *
1239 * Do Tests 9 and 10
1240 *
1241  CALL dstt21( n, 0, sd, se, d1, dumma, z, ldu, work,
1242  $ result( 9 ) )
1243 *
1244 * Do Tests 11 and 12
1245 *
1246  temp1 = zero
1247  temp2 = zero
1248  temp3 = zero
1249  temp4 = zero
1250 *
1251  DO 150 j = 1, n
1252  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
1253  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1254  temp3 = max( temp3, abs( d1( j ) ), abs( d3( j ) ) )
1255  temp4 = max( temp4, abs( d1( j )-d3( j ) ) )
1256  150 CONTINUE
1257 *
1258  result( 11 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
1259  result( 12 ) = temp4 / max( unfl, ulp*max( temp3, temp4 ) )
1260 *
1261 * Do Test 13 -- Sturm Sequence Test of Eigenvalues
1262 * Go up by factors of two until it succeeds
1263 *
1264  ntest = 13
1265  temp1 = thresh*( half-ulp )
1266 *
1267  DO 160 j = 0, log2ui
1268  CALL dstech( n, sd, se, d1, temp1, work, iinfo )
1269  IF( iinfo.EQ.0 )
1270  $ GO TO 170
1271  temp1 = temp1*two
1272  160 CONTINUE
1273 *
1274  170 CONTINUE
1275  result( 13 ) = temp1
1276 *
1277 * For positive definite matrices ( JTYPE.GT.15 ) call DPTEQR
1278 * and do tests 14, 15, and 16 .
1279 *
1280  IF( jtype.GT.15 ) THEN
1281 *
1282 * Compute D4 and Z4
1283 *
1284  CALL dcopy( n, sd, 1, d4, 1 )
1285  IF( n.GT.0 )
1286  $ CALL dcopy( n-1, se, 1, work, 1 )
1287  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1288 *
1289  ntest = 14
1290  CALL dpteqr( 'V', n, d4, work, z, ldu, work( n+1 ),
1291  $ iinfo )
1292  IF( iinfo.NE.0 ) THEN
1293  WRITE( nounit, fmt = 9999 )'DPTEQR(V)', iinfo, n,
1294  $ jtype, ioldsd
1295  info = abs( iinfo )
1296  IF( iinfo.LT.0 ) THEN
1297  RETURN
1298  ELSE
1299  result( 14 ) = ulpinv
1300  GO TO 280
1301  END IF
1302  END IF
1303 *
1304 * Do Tests 14 and 15
1305 *
1306  CALL dstt21( n, 0, sd, se, d4, dumma, z, ldu, work,
1307  $ result( 14 ) )
1308 *
1309 * Compute D5
1310 *
1311  CALL dcopy( n, sd, 1, d5, 1 )
1312  IF( n.GT.0 )
1313  $ CALL dcopy( n-1, se, 1, work, 1 )
1314 *
1315  ntest = 16
1316  CALL dpteqr( 'N', n, d5, work, z, ldu, work( n+1 ),
1317  $ iinfo )
1318  IF( iinfo.NE.0 ) THEN
1319  WRITE( nounit, fmt = 9999 )'DPTEQR(N)', iinfo, n,
1320  $ jtype, ioldsd
1321  info = abs( iinfo )
1322  IF( iinfo.LT.0 ) THEN
1323  RETURN
1324  ELSE
1325  result( 16 ) = ulpinv
1326  GO TO 280
1327  END IF
1328  END IF
1329 *
1330 * Do Test 16
1331 *
1332  temp1 = zero
1333  temp2 = zero
1334  DO 180 j = 1, n
1335  temp1 = max( temp1, abs( d4( j ) ), abs( d5( j ) ) )
1336  temp2 = max( temp2, abs( d4( j )-d5( j ) ) )
1337  180 CONTINUE
1338 *
1339  result( 16 ) = temp2 / max( unfl,
1340  $ hun*ulp*max( temp1, temp2 ) )
1341  ELSE
1342  result( 14 ) = zero
1343  result( 15 ) = zero
1344  result( 16 ) = zero
1345  END IF
1346 *
1347 * Call DSTEBZ with different options and do tests 17-18.
1348 *
1349 * If S is positive definite and diagonally dominant,
1350 * ask for all eigenvalues with high relative accuracy.
1351 *
1352  vl = zero
1353  vu = zero
1354  il = 0
1355  iu = 0
1356  IF( jtype.EQ.21 ) THEN
1357  ntest = 17
1358  abstol = unfl + unfl
1359  CALL dstebz( 'A', 'E', n, vl, vu, il, iu, abstol, sd, se,
1360  $ m, nsplit, wr, iwork( 1 ), iwork( n+1 ),
1361  $ work, iwork( 2*n+1 ), iinfo )
1362  IF( iinfo.NE.0 ) THEN
1363  WRITE( nounit, fmt = 9999 )'DSTEBZ(A,rel)', iinfo, n,
1364  $ jtype, ioldsd
1365  info = abs( iinfo )
1366  IF( iinfo.LT.0 ) THEN
1367  RETURN
1368  ELSE
1369  result( 17 ) = ulpinv
1370  GO TO 280
1371  END IF
1372  END IF
1373 *
1374 * Do test 17
1375 *
1376  temp2 = two*( two*n-one )*ulp*( one+eight*half**2 ) /
1377  $ ( one-half )**4
1378 *
1379  temp1 = zero
1380  DO 190 j = 1, n
1381  temp1 = max( temp1, abs( d4( j )-wr( n-j+1 ) ) /
1382  $ ( abstol+abs( d4( j ) ) ) )
1383  190 CONTINUE
1384 *
1385  result( 17 ) = temp1 / temp2
1386  ELSE
1387  result( 17 ) = zero
1388  END IF
1389 *
1390 * Now ask for all eigenvalues with high absolute accuracy.
1391 *
1392  ntest = 18
1393  abstol = unfl + unfl
1394  CALL dstebz( 'A', 'E', n, vl, vu, il, iu, abstol, sd, se, m,
1395  $ nsplit, wa1, iwork( 1 ), iwork( n+1 ), work,
1396  $ iwork( 2*n+1 ), iinfo )
1397  IF( iinfo.NE.0 ) THEN
1398  WRITE( nounit, fmt = 9999 )'DSTEBZ(A)', iinfo, n, jtype,
1399  $ ioldsd
1400  info = abs( iinfo )
1401  IF( iinfo.LT.0 ) THEN
1402  RETURN
1403  ELSE
1404  result( 18 ) = ulpinv
1405  GO TO 280
1406  END IF
1407  END IF
1408 *
1409 * Do test 18
1410 *
1411  temp1 = zero
1412  temp2 = zero
1413  DO 200 j = 1, n
1414  temp1 = max( temp1, abs( d3( j ) ), abs( wa1( j ) ) )
1415  temp2 = max( temp2, abs( d3( j )-wa1( j ) ) )
1416  200 CONTINUE
1417 *
1418  result( 18 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
1419 *
1420 * Choose random values for IL and IU, and ask for the
1421 * IL-th through IU-th eigenvalues.
1422 *
1423  ntest = 19
1424  IF( n.LE.1 ) THEN
1425  il = 1
1426  iu = n
1427  ELSE
1428  il = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1429  iu = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1430  IF( iu.LT.il ) THEN
1431  itemp = iu
1432  iu = il
1433  il = itemp
1434  END IF
1435  END IF
1436 *
1437  CALL dstebz( 'I', 'E', n, vl, vu, il, iu, abstol, sd, se,
1438  $ m2, nsplit, wa2, iwork( 1 ), iwork( n+1 ),
1439  $ work, iwork( 2*n+1 ), iinfo )
1440  IF( iinfo.NE.0 ) THEN
1441  WRITE( nounit, fmt = 9999 )'DSTEBZ(I)', iinfo, n, jtype,
1442  $ ioldsd
1443  info = abs( iinfo )
1444  IF( iinfo.LT.0 ) THEN
1445  RETURN
1446  ELSE
1447  result( 19 ) = ulpinv
1448  GO TO 280
1449  END IF
1450  END IF
1451 *
1452 * Determine the values VL and VU of the IL-th and IU-th
1453 * eigenvalues and ask for all eigenvalues in this range.
1454 *
1455  IF( n.GT.0 ) THEN
1456  IF( il.NE.1 ) THEN
1457  vl = wa1( il ) - max( half*( wa1( il )-wa1( il-1 ) ),
1458  $ ulp*anorm, two*rtunfl )
1459  ELSE
1460  vl = wa1( 1 ) - max( half*( wa1( n )-wa1( 1 ) ),
1461  $ ulp*anorm, two*rtunfl )
1462  END IF
1463  IF( iu.NE.n ) THEN
1464  vu = wa1( iu ) + max( half*( wa1( iu+1 )-wa1( iu ) ),
1465  $ ulp*anorm, two*rtunfl )
1466  ELSE
1467  vu = wa1( n ) + max( half*( wa1( n )-wa1( 1 ) ),
1468  $ ulp*anorm, two*rtunfl )
1469  END IF
1470  ELSE
1471  vl = zero
1472  vu = one
1473  END IF
1474 *
1475  CALL dstebz( 'V', 'E', n, vl, vu, il, iu, abstol, sd, se,
1476  $ m3, nsplit, wa3, iwork( 1 ), iwork( n+1 ),
1477  $ work, iwork( 2*n+1 ), iinfo )
1478  IF( iinfo.NE.0 ) THEN
1479  WRITE( nounit, fmt = 9999 )'DSTEBZ(V)', iinfo, n, jtype,
1480  $ ioldsd
1481  info = abs( iinfo )
1482  IF( iinfo.LT.0 ) THEN
1483  RETURN
1484  ELSE
1485  result( 19 ) = ulpinv
1486  GO TO 280
1487  END IF
1488  END IF
1489 *
1490  IF( m3.EQ.0 .AND. n.NE.0 ) THEN
1491  result( 19 ) = ulpinv
1492  GO TO 280
1493  END IF
1494 *
1495 * Do test 19
1496 *
1497  temp1 = dsxt1( 1, wa2, m2, wa3, m3, abstol, ulp, unfl )
1498  temp2 = dsxt1( 1, wa3, m3, wa2, m2, abstol, ulp, unfl )
1499  IF( n.GT.0 ) THEN
1500  temp3 = max( abs( wa1( n ) ), abs( wa1( 1 ) ) )
1501  ELSE
1502  temp3 = zero
1503  END IF
1504 *
1505  result( 19 ) = ( temp1+temp2 ) / max( unfl, temp3*ulp )
1506 *
1507 * Call DSTEIN to compute eigenvectors corresponding to
1508 * eigenvalues in WA1. (First call DSTEBZ again, to make sure
1509 * it returns these eigenvalues in the correct order.)
1510 *
1511  ntest = 21
1512  CALL dstebz( 'A', 'B', n, vl, vu, il, iu, abstol, sd, se, m,
1513  $ nsplit, wa1, iwork( 1 ), iwork( n+1 ), work,
1514  $ iwork( 2*n+1 ), iinfo )
1515  IF( iinfo.NE.0 ) THEN
1516  WRITE( nounit, fmt = 9999 )'DSTEBZ(A,B)', iinfo, n,
1517  $ jtype, ioldsd
1518  info = abs( iinfo )
1519  IF( iinfo.LT.0 ) THEN
1520  RETURN
1521  ELSE
1522  result( 20 ) = ulpinv
1523  result( 21 ) = ulpinv
1524  GO TO 280
1525  END IF
1526  END IF
1527 *
1528  CALL dstein( n, sd, se, m, wa1, iwork( 1 ), iwork( n+1 ), z,
1529  $ ldu, work, iwork( 2*n+1 ), iwork( 3*n+1 ),
1530  $ iinfo )
1531  IF( iinfo.NE.0 ) THEN
1532  WRITE( nounit, fmt = 9999 )'DSTEIN', iinfo, n, jtype,
1533  $ ioldsd
1534  info = abs( iinfo )
1535  IF( iinfo.LT.0 ) THEN
1536  RETURN
1537  ELSE
1538  result( 20 ) = ulpinv
1539  result( 21 ) = ulpinv
1540  GO TO 280
1541  END IF
1542  END IF
1543 *
1544 * Do tests 20 and 21
1545 *
1546  CALL dstt21( n, 0, sd, se, wa1, dumma, z, ldu, work,
1547  $ result( 20 ) )
1548 *
1549 * Call DSTEDC(I) to compute D1 and Z, do tests.
1550 *
1551 * Compute D1 and Z
1552 *
1553  CALL dcopy( n, sd, 1, d1, 1 )
1554  IF( n.GT.0 )
1555  $ CALL dcopy( n-1, se, 1, work, 1 )
1556  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1557 *
1558  ntest = 22
1559  CALL dstedc( 'I', n, d1, work, z, ldu, work( n+1 ), lwedc-n,
1560  $ iwork, liwedc, iinfo )
1561  IF( iinfo.NE.0 ) THEN
1562  WRITE( nounit, fmt = 9999 )'DSTEDC(I)', iinfo, n, jtype,
1563  $ ioldsd
1564  info = abs( iinfo )
1565  IF( iinfo.LT.0 ) THEN
1566  RETURN
1567  ELSE
1568  result( 22 ) = ulpinv
1569  GO TO 280
1570  END IF
1571  END IF
1572 *
1573 * Do Tests 22 and 23
1574 *
1575  CALL dstt21( n, 0, sd, se, d1, dumma, z, ldu, work,
1576  $ result( 22 ) )
1577 *
1578 * Call DSTEDC(V) to compute D1 and Z, do tests.
1579 *
1580 * Compute D1 and Z
1581 *
1582  CALL dcopy( n, sd, 1, d1, 1 )
1583  IF( n.GT.0 )
1584  $ CALL dcopy( n-1, se, 1, work, 1 )
1585  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1586 *
1587  ntest = 24
1588  CALL dstedc( 'V', n, d1, work, z, ldu, work( n+1 ), lwedc-n,
1589  $ iwork, liwedc, iinfo )
1590  IF( iinfo.NE.0 ) THEN
1591  WRITE( nounit, fmt = 9999 )'DSTEDC(V)', iinfo, n, jtype,
1592  $ ioldsd
1593  info = abs( iinfo )
1594  IF( iinfo.LT.0 ) THEN
1595  RETURN
1596  ELSE
1597  result( 24 ) = ulpinv
1598  GO TO 280
1599  END IF
1600  END IF
1601 *
1602 * Do Tests 24 and 25
1603 *
1604  CALL dstt21( n, 0, sd, se, d1, dumma, z, ldu, work,
1605  $ result( 24 ) )
1606 *
1607 * Call DSTEDC(N) to compute D2, do tests.
1608 *
1609 * Compute D2
1610 *
1611  CALL dcopy( n, sd, 1, d2, 1 )
1612  IF( n.GT.0 )
1613  $ CALL dcopy( n-1, se, 1, work, 1 )
1614  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1615 *
1616  ntest = 26
1617  CALL dstedc( 'N', n, d2, work, z, ldu, work( n+1 ), lwedc-n,
1618  $ iwork, liwedc, iinfo )
1619  IF( iinfo.NE.0 ) THEN
1620  WRITE( nounit, fmt = 9999 )'DSTEDC(N)', iinfo, n, jtype,
1621  $ ioldsd
1622  info = abs( iinfo )
1623  IF( iinfo.LT.0 ) THEN
1624  RETURN
1625  ELSE
1626  result( 26 ) = ulpinv
1627  GO TO 280
1628  END IF
1629  END IF
1630 *
1631 * Do Test 26
1632 *
1633  temp1 = zero
1634  temp2 = zero
1635 *
1636  DO 210 j = 1, n
1637  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
1638  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1639  210 CONTINUE
1640 *
1641  result( 26 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
1642 *
1643 * Only test DSTEMR if IEEE compliant
1644 *
1645  IF( ilaenv( 10, 'DSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND.
1646  $ ilaenv( 11, 'DSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN
1647 *
1648 * Call DSTEMR, do test 27 (relative eigenvalue accuracy)
1649 *
1650 * If S is positive definite and diagonally dominant,
1651 * ask for all eigenvalues with high relative accuracy.
1652 *
1653  vl = zero
1654  vu = zero
1655  il = 0
1656  iu = 0
1657  IF( jtype.EQ.21 .AND. srel ) THEN
1658  ntest = 27
1659  abstol = unfl + unfl
1660  CALL dstemr( 'V', 'A', n, sd, se, vl, vu, il, iu,
1661  $ m, wr, z, ldu, n, iwork( 1 ), tryrac,
1662  $ work, lwork, iwork( 2*n+1 ), lwork-2*n,
1663  $ iinfo )
1664  IF( iinfo.NE.0 ) THEN
1665  WRITE( nounit, fmt = 9999 )'DSTEMR(V,A,rel)',
1666  $ iinfo, n, jtype, ioldsd
1667  info = abs( iinfo )
1668  IF( iinfo.LT.0 ) THEN
1669  RETURN
1670  ELSE
1671  result( 27 ) = ulpinv
1672  GO TO 270
1673  END IF
1674  END IF
1675 *
1676 * Do test 27
1677 *
1678  temp2 = two*( two*n-one )*ulp*( one+eight*half**2 ) /
1679  $ ( one-half )**4
1680 *
1681  temp1 = zero
1682  DO 220 j = 1, n
1683  temp1 = max( temp1, abs( d4( j )-wr( n-j+1 ) ) /
1684  $ ( abstol+abs( d4( j ) ) ) )
1685  220 CONTINUE
1686 *
1687  result( 27 ) = temp1 / temp2
1688 *
1689  il = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1690  iu = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1691  IF( iu.LT.il ) THEN
1692  itemp = iu
1693  iu = il
1694  il = itemp
1695  END IF
1696 *
1697  IF( srange ) THEN
1698  ntest = 28
1699  abstol = unfl + unfl
1700  CALL dstemr( 'V', 'I', n, sd, se, vl, vu, il, iu,
1701  $ m, wr, z, ldu, n, iwork( 1 ), tryrac,
1702  $ work, lwork, iwork( 2*n+1 ),
1703  $ lwork-2*n, iinfo )
1704 *
1705  IF( iinfo.NE.0 ) THEN
1706  WRITE( nounit, fmt = 9999 )'DSTEMR(V,I,rel)',
1707  $ iinfo, n, jtype, ioldsd
1708  info = abs( iinfo )
1709  IF( iinfo.LT.0 ) THEN
1710  RETURN
1711  ELSE
1712  result( 28 ) = ulpinv
1713  GO TO 270
1714  END IF
1715  END IF
1716 *
1717 * Do test 28
1718 *
1719  temp2 = two*( two*n-one )*ulp*
1720  $ ( one+eight*half**2 ) / ( one-half )**4
1721 *
1722  temp1 = zero
1723  DO 230 j = il, iu
1724  temp1 = max( temp1, abs( wr( j-il+1 )-d4( n-j+
1725  $ 1 ) ) / ( abstol+abs( wr( j-il+1 ) ) ) )
1726  230 CONTINUE
1727 *
1728  result( 28 ) = temp1 / temp2
1729  ELSE
1730  result( 28 ) = zero
1731  END IF
1732  ELSE
1733  result( 27 ) = zero
1734  result( 28 ) = zero
1735  END IF
1736 *
1737 * Call DSTEMR(V,I) to compute D1 and Z, do tests.
1738 *
1739 * Compute D1 and Z
1740 *
1741  CALL dcopy( n, sd, 1, d5, 1 )
1742  IF( n.GT.0 )
1743  $ CALL dcopy( n-1, se, 1, work, 1 )
1744  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1745 *
1746  IF( srange ) THEN
1747  ntest = 29
1748  il = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1749  iu = 1 + ( n-1 )*int( dlarnd( 1, iseed2 ) )
1750  IF( iu.LT.il ) THEN
1751  itemp = iu
1752  iu = il
1753  il = itemp
1754  END IF
1755  CALL dstemr( 'V', 'I', n, d5, work, vl, vu, il, iu,
1756  $ m, d1, z, ldu, n, iwork( 1 ), tryrac,
1757  $ work( n+1 ), lwork-n, iwork( 2*n+1 ),
1758  $ liwork-2*n, iinfo )
1759  IF( iinfo.NE.0 ) THEN
1760  WRITE( nounit, fmt = 9999 )'DSTEMR(V,I)', iinfo,
1761  $ n, jtype, ioldsd
1762  info = abs( iinfo )
1763  IF( iinfo.LT.0 ) THEN
1764  RETURN
1765  ELSE
1766  result( 29 ) = ulpinv
1767  GO TO 280
1768  END IF
1769  END IF
1770 *
1771 * Do Tests 29 and 30
1772 *
1773  CALL dstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work,
1774  $ m, result( 29 ) )
1775 *
1776 * Call DSTEMR to compute D2, do tests.
1777 *
1778 * Compute D2
1779 *
1780  CALL dcopy( n, sd, 1, d5, 1 )
1781  IF( n.GT.0 )
1782  $ CALL dcopy( n-1, se, 1, work, 1 )
1783 *
1784  ntest = 31
1785  CALL dstemr( 'N', 'I', n, d5, work, vl, vu, il, iu,
1786  $ m, d2, z, ldu, n, iwork( 1 ), tryrac,
1787  $ work( n+1 ), lwork-n, iwork( 2*n+1 ),
1788  $ liwork-2*n, iinfo )
1789  IF( iinfo.NE.0 ) THEN
1790  WRITE( nounit, fmt = 9999 )'DSTEMR(N,I)', iinfo,
1791  $ n, jtype, ioldsd
1792  info = abs( iinfo )
1793  IF( iinfo.LT.0 ) THEN
1794  RETURN
1795  ELSE
1796  result( 31 ) = ulpinv
1797  GO TO 280
1798  END IF
1799  END IF
1800 *
1801 * Do Test 31
1802 *
1803  temp1 = zero
1804  temp2 = zero
1805 *
1806  DO 240 j = 1, iu - il + 1
1807  temp1 = max( temp1, abs( d1( j ) ),
1808  $ abs( d2( j ) ) )
1809  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1810  240 CONTINUE
1811 *
1812  result( 31 ) = temp2 / max( unfl,
1813  $ ulp*max( temp1, temp2 ) )
1814 *
1815 * Call DSTEMR(V,V) to compute D1 and Z, do tests.
1816 *
1817 * Compute D1 and Z
1818 *
1819  CALL dcopy( n, sd, 1, d5, 1 )
1820  IF( n.GT.0 )
1821  $ CALL dcopy( n-1, se, 1, work, 1 )
1822  CALL dlaset( 'Full', n, n, zero, one, z, ldu )
1823 *
1824  ntest = 32
1825 *
1826  IF( n.GT.0 ) THEN
1827  IF( il.NE.1 ) THEN
1828  vl = d2( il ) - max( half*
1829  $ ( d2( il )-d2( il-1 ) ), ulp*anorm,
1830  $ two*rtunfl )
1831  ELSE
1832  vl = d2( 1 ) - max( half*( d2( n )-d2( 1 ) ),
1833  $ ulp*anorm, two*rtunfl )
1834  END IF
1835  IF( iu.NE.n ) THEN
1836  vu = d2( iu ) + max( half*
1837  $ ( d2( iu+1 )-d2( iu ) ), ulp*anorm,
1838  $ two*rtunfl )
1839  ELSE
1840  vu = d2( n ) + max( half*( d2( n )-d2( 1 ) ),
1841  $ ulp*anorm, two*rtunfl )
1842  END IF
1843  ELSE
1844  vl = zero
1845  vu = one
1846  END IF
1847 *
1848  CALL dstemr( 'V', 'V', n, d5, work, vl, vu, il, iu,
1849  $ m, d1, z, ldu, n, iwork( 1 ), tryrac,
1850  $ work( n+1 ), lwork-n, iwork( 2*n+1 ),
1851  $ liwork-2*n, iinfo )
1852  IF( iinfo.NE.0 ) THEN
1853  WRITE( nounit, fmt = 9999 )'DSTEMR(V,V)', iinfo,
1854  $ n, jtype, ioldsd
1855  info = abs( iinfo )
1856  IF( iinfo.LT.0 ) THEN
1857  RETURN
1858  ELSE
1859  result( 32 ) = ulpinv
1860  GO TO 280
1861  END IF
1862  END IF
1863 *
1864 * Do Tests 32 and 33
1865 *
1866  CALL dstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work,
1867  $ m, result( 32 ) )
1868 *
1869 * Call DSTEMR to compute D2, do tests.
1870 *
1871 * Compute D2
1872 *
1873  CALL dcopy( n, sd, 1, d5, 1 )
1874  IF( n.GT.0 )
1875  $ CALL dcopy( n-1, se, 1, work, 1 )
1876 *
1877  ntest = 34
1878  CALL dstemr( 'N', 'V', n, d5, work, vl, vu, il, iu,
1879  $ m, d2, z, ldu, n, iwork( 1 ), tryrac,
1880  $ work( n+1 ), lwork-n, iwork( 2*n+1 ),
1881  $ liwork-2*n, iinfo )
1882  IF( iinfo.NE.0 ) THEN
1883  WRITE( nounit, fmt = 9999 )'DSTEMR(N,V)', iinfo,
1884  $ n, jtype, ioldsd
1885  info = abs( iinfo )
1886  IF( iinfo.LT.0 ) THEN
1887  RETURN
1888  ELSE
1889  result( 34 ) = ulpinv
1890  GO TO 280
1891  END IF
1892  END IF
1893 *
1894 * Do Test 34
1895 *
1896  temp1 = zero
1897  temp2 = zero
1898 *
1899  DO 250 j = 1, iu - il + 1
1900  temp1 = max( temp1, abs( d1( j ) ),
1901  $ abs( d2( j ) ) )
1902  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1903  250 CONTINUE
1904 *
1905  result( 34 ) = temp2 / max( unfl,
1906  $ ulp*max( temp1, temp2 ) )
1907  ELSE
1908  result( 29 ) = zero
1909  result( 30 ) = zero
1910  result( 31 ) = zero
1911  result( 32 ) = zero
1912  result( 33 ) = zero
1913  result( 34 ) = zero
1914  END IF
1915 *
1916 * Call DSTEMR(V,A) to compute D1 and Z, do tests.
1917 *
1918 * Compute D1 and Z
1919 *
1920  CALL dcopy( n, sd, 1, d5, 1 )
1921  IF( n.GT.0 )
1922  $ CALL dcopy( n-1, se, 1, work, 1 )
1923 *
1924  ntest = 35
1925 *
1926  CALL dstemr( 'V', 'A', n, d5, work, vl, vu, il, iu,
1927  $ m, d1, z, ldu, n, iwork( 1 ), tryrac,
1928  $ work( n+1 ), lwork-n, iwork( 2*n+1 ),
1929  $ liwork-2*n, iinfo )
1930  IF( iinfo.NE.0 ) THEN
1931  WRITE( nounit, fmt = 9999 )'DSTEMR(V,A)', iinfo, n,
1932  $ jtype, ioldsd
1933  info = abs( iinfo )
1934  IF( iinfo.LT.0 ) THEN
1935  RETURN
1936  ELSE
1937  result( 35 ) = ulpinv
1938  GO TO 280
1939  END IF
1940  END IF
1941 *
1942 * Do Tests 35 and 36
1943 *
1944  CALL dstt22( n, m, 0, sd, se, d1, dumma, z, ldu, work, m,
1945  $ result( 35 ) )
1946 *
1947 * Call DSTEMR to compute D2, do tests.
1948 *
1949 * Compute D2
1950 *
1951  CALL dcopy( n, sd, 1, d5, 1 )
1952  IF( n.GT.0 )
1953  $ CALL dcopy( n-1, se, 1, work, 1 )
1954 *
1955  ntest = 37
1956  CALL dstemr( 'N', 'A', n, d5, work, vl, vu, il, iu,
1957  $ m, d2, z, ldu, n, iwork( 1 ), tryrac,
1958  $ work( n+1 ), lwork-n, iwork( 2*n+1 ),
1959  $ liwork-2*n, iinfo )
1960  IF( iinfo.NE.0 ) THEN
1961  WRITE( nounit, fmt = 9999 )'DSTEMR(N,A)', iinfo, n,
1962  $ jtype, ioldsd
1963  info = abs( iinfo )
1964  IF( iinfo.LT.0 ) THEN
1965  RETURN
1966  ELSE
1967  result( 37 ) = ulpinv
1968  GO TO 280
1969  END IF
1970  END IF
1971 *
1972 * Do Test 37
1973 *
1974  temp1 = zero
1975  temp2 = zero
1976 *
1977  DO 260 j = 1, n
1978  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
1979  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
1980  260 CONTINUE
1981 *
1982  result( 37 ) = temp2 / max( unfl,
1983  $ ulp*max( temp1, temp2 ) )
1984  END IF
1985  270 CONTINUE
1986  280 CONTINUE
1987  ntestt = ntestt + ntest
1988 *
1989 * End of Loop -- Check for RESULT(j) > THRESH
1990 *
1991 * Print out tests which fail.
1992 *
1993  DO 290 jr = 1, ntest
1994  IF( result( jr ).GE.thresh ) THEN
1995 *
1996 * If this is the first test to fail,
1997 * print a header to the data file.
1998 *
1999  IF( nerrs.EQ.0 ) THEN
2000  WRITE( nounit, fmt = 9998 )'DST'
2001  WRITE( nounit, fmt = 9997 )
2002  WRITE( nounit, fmt = 9996 )
2003  WRITE( nounit, fmt = 9995 )'Symmetric'
2004  WRITE( nounit, fmt = 9994 )
2005 *
2006 * Tests performed
2007 *
2008  WRITE( nounit, fmt = 9988 )
2009  END IF
2010  nerrs = nerrs + 1
2011  WRITE( nounit, fmt = 9990 )n, ioldsd, jtype, jr,
2012  $ result( jr )
2013  END IF
2014  290 CONTINUE
2015  300 CONTINUE
2016  310 CONTINUE
2017 *
2018 * Summary
2019 *
2020  CALL dlasum( 'DST', nounit, nerrs, ntestt )
2021  RETURN
2022 *
2023  9999 FORMAT( ' DCHKST2STG: ', a, ' returned INFO=', i6, '.', / 9x,
2024  $ 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
2025 *
2026  9998 FORMAT( / 1x, a3, ' -- Real Symmetric eigenvalue problem' )
2027  9997 FORMAT( ' Matrix types (see DCHKST2STG for details): ' )
2028 *
2029  9996 FORMAT( / ' Special Matrices:',
2030  $ / ' 1=Zero matrix. ',
2031  $ ' 5=Diagonal: clustered entries.',
2032  $ / ' 2=Identity matrix. ',
2033  $ ' 6=Diagonal: large, evenly spaced.',
2034  $ / ' 3=Diagonal: evenly spaced entries. ',
2035  $ ' 7=Diagonal: small, evenly spaced.',
2036  $ / ' 4=Diagonal: geometr. spaced entries.' )
2037  9995 FORMAT( ' Dense ', a, ' Matrices:',
2038  $ / ' 8=Evenly spaced eigenvals. ',
2039  $ ' 12=Small, evenly spaced eigenvals.',
2040  $ / ' 9=Geometrically spaced eigenvals. ',
2041  $ ' 13=Matrix with random O(1) entries.',
2042  $ / ' 10=Clustered eigenvalues. ',
2043  $ ' 14=Matrix with large random entries.',
2044  $ / ' 11=Large, evenly spaced eigenvals. ',
2045  $ ' 15=Matrix with small random entries.' )
2046  9994 FORMAT( ' 16=Positive definite, evenly spaced eigenvalues',
2047  $ / ' 17=Positive definite, geometrically spaced eigenvlaues',
2048  $ / ' 18=Positive definite, clustered eigenvalues',
2049  $ / ' 19=Positive definite, small evenly spaced eigenvalues',
2050  $ / ' 20=Positive definite, large evenly spaced eigenvalues',
2051  $ / ' 21=Diagonally dominant tridiagonal, geometrically',
2052  $ ' spaced eigenvalues' )
2053 *
2054  9990 FORMAT( ' N=', i5, ', seed=', 4( i4, ',' ), ' type ', i2,
2055  $ ', test(', i2, ')=', g10.3 )
2056 *
2057  9988 FORMAT( / 'Test performed: see DCHKST2STG for details.', / )
2058 *
2059 * End of DCHKST2STG
2060 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:273
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DSTEQR
Definition: dsteqr.f:131
subroutine dstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEDC
Definition: dstedc.f:188
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
double precision function dsxt1(IJOB, D1, N1, D2, N2, ABSTOL, ULP, UNFL)
DSXT1
Definition: dsxt1.f:106
subroutine dlasum(TYPE, IOUNIT, IE, NRUN)
DLASUM
Definition: dlasum.f:43
subroutine dstt22(N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK, LDWORK, RESULT)
DSTT22
Definition: dstt22.f:139
subroutine dspt21(ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, TAU, WORK, RESULT)
DSPT21
Definition: dspt21.f:221
subroutine dstech(N, A, B, EIG, TOL, WORK, INFO)
DSTECH
Definition: dstech.f:101
subroutine dsyt21(ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, LDV, TAU, WORK, RESULT)
DSYT21
Definition: dsyt21.f:207
subroutine dstt21(N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RESULT)
DSTT21
Definition: dstt21.f:127
subroutine dlatmr(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, RSIGN, GRADE, DL, MODEL, CONDL, DR, MODER, CONDR, PIVTNG, IPIVOT, KL, KU, SPARSE, ANORM, PACK, A, LDA, IWORK, INFO)
DLATMR
Definition: dlatmr.f:471
double precision function dlarnd(IDIST, ISEED)
DLARND
Definition: dlarnd.f:73
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
Definition: dlatms.f:321
subroutine dstemr(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEMR
Definition: dstemr.f:321
subroutine dsptrd(UPLO, N, AP, D, E, TAU, INFO)
DSPTRD
Definition: dsptrd.f:150
subroutine dstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
DSTEIN
Definition: dstein.f:174
subroutine dopgtr(UPLO, N, AP, TAU, Q, LDQ, WORK, INFO)
DOPGTR
Definition: dopgtr.f:114
subroutine dorgtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
DORGTR
Definition: dorgtr.f:123
subroutine dpteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
DPTEQR
Definition: dpteqr.f:145
subroutine dsytrd_2stage(VECT, UPLO, N, A, LDA, D, E, TAU, HOUS2, LHOUS2, WORK, LWORK, INFO)
DSYTRD_2STAGE
subroutine dsytrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
DSYTRD
Definition: dsytrd.f:192
Here is the call graph for this function:
Here is the caller graph for this function: