Purpose
 =======
 
       LA_SYEVR / LA_HEEVR compute selected eigenvalues and, optionally,
 the corresponding eigenvectors of a real symmetric/complex Hermitian 
 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_SYEVR and LA_HEEVR use a relatively robust representation 
 (RRR) algorithm. It is usually the fastest algorithm of all and uses 
 the least workspace.
 
 ========
 
        SUBROUTINE LA SYEVR / LA HEEVR( A, W, JOBZ=jobz, &
                 UPLO=uplo, VL=vl, VU=vu, IL=il, IU=iu, M=m, &
                 ISUPPZ= isuppz, ABSTOL=abstol, INFO=info )
           <type>(<wp>), INTENT(INOUT) :: A(:,:)
           REAL(<wp>), INTENT(OUT) :: W(:)
           CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: JOBZ, UPLO
           REAL(<wp>), INTENT(IN), OPTIONAL :: VL, VU
           INTEGER, INTENT(IN), OPTIONAL :: IL, IU
           INTEGER, INTENT(OUT), OPTIONAL :: M
           INTEGER, INTENT(OUT), OPTIONAL :: ISUPPZ(:)
           REAL(<wp>), INTENT(IN), OPTIONAL :: ABSTOL
           INTEGER, INTENT(OUT), OPTIONAL :: INFO
        where
           <type> ::= REAL | COMPLEX
           <wp>   ::= KIND(1.0) | KIND(1.0D0)
 
 Arguments
 =========
 
 A       (input/output) REAL or COMPLEX square array, shape (:,:).
         On entry, the matrix A.
         If UPLO = 'U', the upper triangular part of A contains the 
         upper triangular part of the matrix A. If UPLO = 'L', the 
         lower triangular part of A contains the lower triangular part
         of the matrix A.
         On exit:
         If JOBZ = 'V', then the first M columns of A contain the 
         orthonormal eigenvectors of the matrix A corresponding to the
         selected eigenvalues, with the i-th column of A containing the
         eigenvector associated with the eigenvalue in W(i).
         If JOBZ = 'N', the upper triangle (if UPLO = 'U') or the lower 
         triangle (if UPLO = 'L') of A, including the diagonal, is 
 	  destroyed.
 W       (output) REAL array, shape (:) with size(W) = size(A,1).
         The first M elements contain the selected eigenvalues in 
         ascending order.
 JOBZ    Optional (input) CHARACTER(LEN=1).
         = 'N': Computes eigenvalues only;
         = 'V': Computes eigenvalues and eigenvectors.
         Default value: 'N'.
 UPLO    Optional (input) CHARACTER(LEN=1).
         = 'U': Upper triangle of A is stored;
         = 'L': Lower triangle of A is stored.
         Default value: 'U'.
 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 <= size(A,1).
         Default values: IL = 1 and IU = size(A,1).
         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 <= size(A,1).
         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(2i-1) through ISUPPZ(2i) .
         Note: ISUPPZ must be absent if JOBZ = 'N'.
 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>)* ||T||1 will be used in its place, where 
 	  ||T||1 is the l1 norm of the tridiagonal matrix obtained by 
         reducing A to tridiagonal form.
         Default value: 0.0_<wp>.
         Note: Eigenvalues are computed most accurately if ABSTOL is set
         to LA_LAMCH(1.0_<wp>,'Safe minimum'), not zero.
 INFO    Optional (output) INTEGER
         = 0: successful exit.
         < 0: if INFO = -i, the i-th argument had an illegal value.
         If INFO is not present and an error occurs, then the program
         is terminated with an error message.