Purpose
 =======
       LA_STEVR computes selected eigenvalues and, optionally, the 
 corresponding eigenvectors of a real symmetric tridiagonal matrix A.
 Eigenvalues and eigenvectors can be selected by specifying either a 
 range of values or a range of indices for the desired eigenvalues.
       LA_STEVR uses a relatively robust representation (RRR) algorithm.
 It is usually the fastest algorithm of all and uses the least
 workspace.
 
 =========
 
         SUBROUTINE LA_STEVR ( D, E, W, Z=z, VL=vl, VU=vu, &
                         IL=il, IU=iu, M=m, ISUPPZ=isuppz, &
                         ABSTOL=abstol, INFO=info )
                 REAL(<wp>), INTENT(INOUT) :: D(:), E(:)
                 REAL(<wp>), INTENT(OUT) :: W(:)
                 REAL(<wp>), INTENT(OUT), OPTIONAL :: Z(:,:)
                 INTEGER, INTENT(OUT), OPTIONAL :: ISUPPZ(:)
                 REAL(<wp>), INTENT(IN), OPTIONAL :: VL, VU
                 INTEGER, INTENT(IN), OPTIONAL :: IL, IU
                 INTEGER, INTENT(OUT), OPTIONAL :: M
                 REAL(<wp>), INTENT(IN), OPTIONAL :: ABSTOL
                 INTEGER, INTENT(OUT), OPTIONAL :: INFO
         where
                 <wp> ::= KIND(1.0) | KIND(1.0D0)
 
 Arguments
 =========
 
 D      (input/output) REAL array, shape (:) with size(D) = n, where n 
        is the order of A.
        On entry, the diagonal elements of the matrix A.
        On exit, the original contents of D possibly multiplied by a 
        constant factor to avoid over/underflow in computing the 
        eigenvalues.
 E      (input/output) REAL array, shape (:) with size(E) = n.
        On entry, the n-1 subdiagonal elements of A in E(1) to E(n-1) .
        E(n) need not be set.
        On exit, the original contents of E possibly multiplied by a 
        constant factor to avoid over/underflow in computing the 
        eigenvalues.
 W      (output) REAL array with size(W) = n.
        The first M elements contain the selected eigenvalues in 
        ascending order.
 Z      Optional (output) REAL or COMPLEX array, shape (:,:) with 
        size(Z,1) = n and size(Z,2) = M.
        The first M columns of Z contain the orthonormal eigenvectors of
        A corresponding to the selected eigenvalues, with the i-th column 
        of Z containing the eigenvector associated with the eigenvalue in
        W(i).
        Note: The user must ensure that at least M columns are supplied 
        in the array Z. When the exact value of M is not known in advance,
        an upper bound must be used. In all cases M <= n.
 VL,VU  Optional (input) REAL.
        The lower and upper bounds of the interval to be searched for 
        eigenvalues. VL < VU.
        Default values: VL = -HUGE(<wp>) and VU = HUGE(<wp>), where 
        <wp> ::= KIND(1.0) | KIND(1.0D0).
        Note: Neither VL nor VU may be present if IL and/or IU is 
        present.
 IL,IU  Optional (input) INTEGER.
        The indices of the smallest and largest eigenvalues to be 
        returned. The IL-th through IU-th eigenvalues will be found. 
        1 <= IL <= IU <= n.
        Default values: IL = 1 and IU = n.
        Note: Neither IL nor IU may be present if VL and/or VU is 
        present.
        Note: All eigenvalues are calculated if none of the arguments 
        VL, VU, IL and IU are present.
 M      Optional (output) INTEGER.
        The total number of eigenvalues found. 0 <= M <=  n.
        Note: If IL and IU are present then M = IU - IL + 1.
 ISUPPZ Optional (output) INTEGER array, shape (:) with 
        size(ISUPPZ) = 2*max(1,M).
        The support of the eigenvectors in A, i.e., the indices 
        indicating the nonzero elements. The i-th eigenvector is nonzero
        only in elements ISUPPZ(2*i-1) through ISUPPZ(2*i).
 ABSTOL Optional (input) REAL.
        The absolute error tolerance for the eigenvalues. An approximate
        eigenvalue is accepted as converged when it is determined to lie
        in an interval [a, b] of width less than or equal to 
        ABSTOL + EPSILON(1.0_<wp>) * max(| a |, | b |),
        where <wp> is the working precision. If ABSTOL <= 0, then 
        EPSILON(1.0_<wp>)*||A||1 will be used in its place. Eigenvalues
        will be computed most accurately if ABSTOL is set to 
        LA_LAMCH( 1.0_<wp>, 'Safe minimum'), not zero.
        Default value: 0.0_<wp>.
 INFO   Optional (output) INTEGER
        = 0: successful exit.
        < 0: if INFO = -i, the i-th argument had an illegal value.
        > 0: an internal error occurred.
        If INFO is not present and an error occurs, then the program is
        terminated with an error message.