LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ zchkst2stg()

subroutine zchkst2stg ( integer  NSIZES,
integer, dimension( * )  NN,
integer  NTYPES,
logical, dimension( * )  DOTYPE,
integer, dimension( 4 )  ISEED,
double precision  THRESH,
integer  NOUNIT,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, 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,
complex*16, dimension( ldu, * )  U,
integer  LDU,
complex*16, dimension( ldu, * )  V,
complex*16, dimension( * )  VP,
complex*16, dimension( * )  TAU,
complex*16, dimension( ldu, * )  Z,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  LRWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
double precision, dimension( * )  RESULT,
integer  INFO 
)

ZCHKST2STG

Purpose:
 ZCHKST2STG  checks the Hermitian 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
 ZHETRD. For that, we call the standard ZHETRD and compute D1 using 
 DSTEQR, then we call the 2-stage ZHETRD_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 ZCHKST in the next 
 release when vectors and generation of Q will be implemented.

    ZHETRD factors A as  U S U* , where * means conjugate transpose,
    S is real symmetric tridiagonal, and U is unitary.
    ZHETRD can use either just the lower or just the upper triangle
    of A; ZCHKST2STG 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.

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

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

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

    ZSTEQR factors S as  Z D1 Z* , where Z is the unitary
    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.

    ZPTEQR factors S as  Z4 D4 Z4* , for a
    Hermitian 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.

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

    ZSTEDC factors S as Z D1 Z* , where Z is the unitary
    matrix of eigenvectors and D1 is a diagonal matrix with
    the eigenvalues on the diagonal ('I' option). It may also
    update an input unitary matrix, usually the output
    from ZHETRD/ZUNGTR or ZHPTRD/ZUPGTR ('V' option). It may
    also just compute eigenvalues ('N' option).

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

 When ZCHKST2STG 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 Hermitian eigenroutines.  For each matrix, a number
 of tests will be performed:

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

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

 (3)     | A - V S V* | / ( |A| n ulp ) ZHETRD( 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
         ZHETRD_2STAGE("N", "U",....). D1 and D2 are computed 
         via DSTEQR('N',...) 

 (4)     | I - UV* | / ( n ulp )        ZUNGTR( 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
         ZHETRD_2STAGE("N", "L",....). D1 and D3 are computed 
         via DSTEQR('N',...)  

 (5-8)   Same as 1-4, but for ZHPTRD and ZUPGTR.

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

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

 (11)    | D1 - D2 | / ( |D1| ulp )        ZSTEQR('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 ) ZPTEQR('V',...)

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

 (16)    | D4 - D5 | / ( 100 |D4| ulp )       ZPTEQR('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, ZSTEIN

 (21)    | I - Y Y* | / ( n ulp )          DSTEBZ, ZSTEIN

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

 (23)    | I - ZZ* | / ( n ulp )           ZSTEDC('I')

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

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

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

 Test 27 is disabled at the moment because ZSTEMR 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
                                              ZSTEMR('V', 'A')

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

 Tests 29 through 34 are disable at present because ZSTEMR
 does not handle partial specturm requests.

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

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

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

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

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

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

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

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

 (37)    ( max { min | WA2(i)-WA3(j) | } +
            i     j
           max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp )
            i     j
         ZSTEMR('N', 'A') vs. CSTEMR('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 unitary 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 unitary 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 unitary 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) Hermitian 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,
          ZCHKST2STG 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, ZCHKST2STG
          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 ZCHKST2STG 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 COMPLEX*16 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 COMPLEX*16 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 ZHETRD.
          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
          ZHETRD.  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 ZSTEQR 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 ZSTEQR 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 ZPTEQR(V).
          ZPTEQR 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 ZPTEQR(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 COMPLEX*16 array of
                             dimension( LDU, max(NN) ).
          The unitary matrix computed by ZHETRD + ZUNGTR.
[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 COMPLEX*16 array of
                             dimension( LDU, max(NN) ).
          The Housholder vectors computed by ZHETRD 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 ZHETRD, 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 ZUNGTR, set those entries to
          1 before using them, and then restore them later.)
[out]VP
          VP is COMPLEX*16 array of
                      dimension( max(NN)*max(NN+1)/2 )
          The matrix V stored in packed format.
[out]TAU
          TAU is COMPLEX*16 array of
                             dimension( max(NN) )
          The Householder factors computed by ZHETRD in reducing A
          to tridiagonal form.
[out]Z
          Z is COMPLEX*16 array of
                             dimension( LDU, max(NN) ).
          The unitary matrix of eigenvectors computed by ZSTEQR,
          ZPTEQR, and ZSTEIN.
[out]WORK
          WORK is COMPLEX*16 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]RWORK
          RWORK is DOUBLE PRECISION array
[in]LRWORK
          LRWORK is INTEGER
          The number of entries in LRWORK (dimension( ??? )
[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  ZLATMR, CLATMS, ZHETRD, ZUNGTR, ZSTEQR, DSTERF,
              or ZUNMC2 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.
Date
December 2016

Definition at line 627 of file zchkst2stg.f.

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