SUBROUTINE SSPGVX_F95( AP, BP, W, ITYPE, UPLO, Z, VL, VU, IL, IU, &
     &  M, IFAIL, ABSTOL, INFO )
!  .. USE STATEMENTS ..
      USE LA_PRECISION, ONLY: WP => SP
      USE LA_AUXMOD, ONLY: ERINFO, LSAME
      USE F77_LAPACK, ONLY: LAMCH_F77 => SLAMCH
      USE F77_LAPACK, ONLY: SPGVX_F77 => LA_SPGVX
!  .. IMPLICIT STATEMENT ..
      IMPLICIT NONE
!  .. CHARACTER ARGUMENTS ..
      CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: UPLO
!  .. SCALAR ARGUMENTS ..
      INTEGER, INTENT(IN), OPTIONAL :: IL, IU, ITYPE
      INTEGER, INTENT(OUT), OPTIONAL :: INFO, M
      REAL(WP), INTENT(IN), OPTIONAL :: ABSTOL, VL, VU
!  .. ARRAY ARGUMENTS ..
      INTEGER, INTENT(OUT), OPTIONAL, TARGET :: IFAIL(:)
      REAL(WP), INTENT(OUT), OPTIONAL, TARGET :: Z(:,:)
      REAL(WP), INTENT(INOUT) :: AP(:), BP(:)
      REAL(WP), INTENT(OUT) :: W(:)
!----------------------------------------------------------------------
! 
! Purpose
! =======
! 
!    LA_SPGVX and LA_HPGVX compute selected eigenvalues and, optionally, 
! the corresponding eigenvectors of generalized eigenvalue problems of 
! the form
!       A*z = lambda*B*z, A*B*z = lambda*z, and B*A*z = lambda*z,
! where A and B are real symmetric in the case of LA_SPGVX and complex
! Hermitian in the case of LA_HPGVX. In both cases B is positive
! definite. Eigenvalues and eigenvectors can be selected by specifying
! either a range of values or a range of indices for the desired 
! eigenvalues. Matrices A and B are stored in a packed format.
! 
! =========
! 
!       SUBROUTINE LA_SPGVX / LA_HPGVX( AP, BP, W, ITYPE= itype, &
!             UPLO= uplo, Z= z, VL= vl, VU= vu, IL= il, IU= iu, M= m, &
!             IFAIL= ifail, ABSTOL= abstol, INFO= info )
!         <type>(<wp>), INTENT(INOUT) :: AP(:), BP(:)
!         REAL(<wp>), INTENT(OUT) :: W(:)
!         INTEGER, INTENT(IN), OPTIONAL :: ITYPE
!         CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: UPLO
!         <type>(<wp>), INTENT(OUT), OPTIONAL :: Z(:,:)
!         REAL(<wp>), INTENT(IN), OPTIONAL :: VL, VU
!         INTEGER, INTENT(IN), OPTIONAL :: IL, IU
!         INTEGER, INTENT(OUT), OPTIONAL :: M
!         INTEGER, INTENT(OUT), OPTIONAL :: IFAIL(:)
!         REAL(<wp>), INTENT(IN), OPTIONAL :: ABSTOL
!         INTEGER, INTENT(OUT), OPTIONAL :: INFO
!       where
!         <type> ::= REAL | COMPLEX
!         <wp>   ::= KIND(1.0) | KIND(1.0D0)
! 
! Arguments
! =========
! 
! AP       (input/output) REAL or COMPLEX array, shape (:) with 
!          size(AP) = n*(n + 1)/2, where n is the order of A and B.
!          On entry, the upper or lower triangle of matrix A in packed
! 	   storage. The elements are stored columnwise as follows:
! 	   if UPLO = 'U', AP(i +(j-1)*j/2) = A(i,j) for 1<=i<=j<=n;
! 	   if UPLO = 'L', AP(i +(j-1)*(2*n-j)/2) = A(i,j) for 1<=j<=i<=n.
!          On exit, the contents of AP are destroyed.
! BP       (input/output) REAL or COMPLEX array, shape (:) with 
!          size(BP) = size(AP).
!          On entry, the upper or lower triangle of matrix B in packed 
! 	   storage. The elements are stored columnwise as follows:
! 	   if UPLO = 'U', BP(i +(j-1)*j/2) = B(i,j) for 1<=i<=j<=n;
! 	   if UPLO = 'L', BP(i +(j-1)*(2*n-j)/2) = B(i,j) for 1<=j<=i<=n.
!          On exit, the triangular factor U or L of the Cholesky 
! 	   factorization B = U^T*U or B = L*L^T, in the same storage 
! 	   format as B.
! W        (output) REAL array, shape (:) with size(W) = n.
!          The eigenvalues in ascending order.
! ITYPE    Optional (input) INTEGER.
!          Specifies the problem type to be solved:
!             = 1: A*z = lambda*B*z
!             = 2: A*B*z = lambda*z
!             = 3: B*A*z = lambda*z
!          Default value: 1.
! UPLO     Optional (input) CHARACTER(LEN=1).
!             = 'U': Upper triangles of A and B are stored;
!             = 'L': Lower triangles of A and B are stored.
!          Default value: 'U'.
! Z        Optional (output) REAL or COMPLEX rectangular array, shape 
!          (:,:) with size(Z,1) = n and size(Z,2) = M.
!          The first M columns of Z contain the orthonormal eigenvectors
! 	   corresponding to the selected eigenvalues, with the i-th 
! 	   column of Z holding the eigenvector associated with the 
! 	   eigenvalue in W(i). The eigenvectors are normalized as 
! 	   follows:
!            if ITYPE = 1 or 2: Z^H * B * Z = I ,
!            if ITYPE = 3: Z^H * B^-1 * Z = I .
!          If an eigenvector fails to converge, then that column of Z
! 	   contains the latest approximation to the eigenvector and the
! 	   index of the eigenvector is returned in IFAIL.
! 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.
! IFAIL    Optional (output) INTEGER array, shape (:) with size(IFAIL) =
!          size(A,1).
!          If INFO = 0, the first M elements of IFAIL are zero.
!          If INFO > 0, then IFAIL contains the indices of the 
! 	   eigenvectors that failed to converge.
!          Note: If Z is present then IFAIL should also be present.
! 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 the generalized eigenvalue problem to tridiagonal
! 	   form. Eigenvalues will be computed most accurately when
! 	   ABSTOL is set to twice the underflow threshold 
! 	   2 * LA_LAMCH(1.0_<wp>, 'S'), not zero.
!          Default value: 0.0_<wp>.
!          Note: If this routine returns with 0 < INFO <= n, then some 
! 	   eigenvectors did not converge.
!          Try setting ABSTOL to 2 * LA_LAMCH(1.0_<wp>, 'S').
! INFO     Optional (output) INTEGER.
!          = 0: successful exit.
!          < 0: if INFO = -i, the i-th argument had an illegal value
!          > 0: the algorithm failed to converge or matrix B is not 
! 	     positive definite:
!             <= n: the algorithm failed to converge; if INFO = i, then 
! 	       i eigenvectors failed to converge. Their indices are 
! 	       stored in array IFAIL.
!             > n: if INFO = n + i, for 1 <= i <= n, then the leading
! 	       minor of order i of B is not positive definite. The 
! 	       factorization of B could not be completed and no
!                eigenvalues or eigenvectors were computed.
!          If INFO is not present and an error occurs, then the program
! 	   is terminated with an error message.
!----------------------------------------------------------------------
!  .. LOCAL PARAMETERS ..
      CHARACTER(LEN=8), PARAMETER :: SRNAME = 'LA_SPGVX'
      CHARACTER(LEN=1) :: LJOBZ, LUPLO, LRANGE
!  .. LOCAL SCALARS ..
      INTEGER :: N, LINFO, LIL, LIU, LM, ISTAT, &
     &  SIFAIL, S1Z, S2Z, NN, LITYPE
      INTEGER, TARGET :: ISTAT1(1)
      REAL(WP), TARGET :: LLZ(1,1)
      REAL(WP) :: LABSTOL, LVL, LVU
      COMPLEX(WP) :: WW
!  .. LOCAL ARRAYS ..
      INTEGER, POINTER :: IWORK(:), LIFAIL(:)
      REAL(WP), POINTER :: WORK(:)
!  .. INTRINSIC FUNCTIONS ..
      INTRINSIC HUGE, PRESENT, SIZE
!  .. EXECUTABLE STATEMENTS ..
   LINFO = 0; ISTAT = 0; NN = SIZE(AP)
   WW = (-1+SQRT(1+8*REAL(NN,WP)))*0.5; N = INT(WW)
   IF( PRESENT(ITYPE) )THEN; LITYPE = ITYPE; ELSE; LITYPE = 1; END IF 
   IF( PRESENT(IFAIL) )THEN; SIFAIL = SIZE(IFAIL); ELSE; SIFAIL = N; END IF
   IF( PRESENT(UPLO) ) THEN; LUPLO = UPLO; ELSE; LUPLO = 'U'; END IF
   IF( PRESENT(VL) )THEN; LVL = VL; ELSE; LVL = -HUGE(LVL); ENDIF
   IF( PRESENT(VU) )THEN; LVU = VU; ELSE; LVU = HUGE(LVU); ENDIF
   IF( PRESENT(IL) )THEN; LIL = IL; ELSE; LIL = 1; ENDIF
   IF( PRESENT(IU) )THEN; LIU = IU; ELSE; LIU = N; ENDIF
   IF( PRESENT(Z) )THEN; S1Z = SIZE(Z,1); S2Z = SIZE(Z,2)
    ELSE; S1Z = 1; S2Z = 1; ENDIF
!  .. TEST THE ARGUMENTS
   IF( PRESENT(VL) .OR. PRESENT(VU) )THEN; LRANGE = 'V'
   ELSE IF( PRESENT(IL) .OR. PRESENT(IU) )THEN; LRANGE = 'I'
   ELSE ; LRANGE = 'A'; END IF                                                 

   IF( NN < 0 .OR. AIMAG(WW) /= 0 .OR. REAL(N,WP) /= REAL(WW) ) THEN; LINFO = -1
   ELSE IF( SIZE(BP) /= SIZE(AP)  )THEN; LINFO = -2 
   ELSE IF( SIZE( W ) /= N )THEN; LINFO = -3
   ELSE IF (LITYPE <1 .OR. LITYPE >3) THEN; LINFO = -4
   ELSE IF( .NOT.LSAME(LUPLO,'U') .AND. .NOT.LSAME(LUPLO,'L') )THEN; LINFO = -5
   ELSE IF( PRESENT(Z) .AND. ( S1Z /= N .OR. S2Z /= N ) )THEN; LINFO = -6
   ELSE IF( LVU < LVL )THEN; LINFO = -7
   ELSE IF( (PRESENT(VL) .OR. PRESENT(VU)) .AND. &
     &  (PRESENT(IL) .OR. PRESENT(IU)) )THEN; LINFO = -8
   ELSE IF( LSAME(LRANGE,'I') .AND. ( LIU < MIN(N, LIL) .OR. LIU > N ))THEN; LINFO = -9
   ELSE IF( N < LIU )THEN; LINFO = -10
   ELSE IF( SIFAIL /= N .OR. PRESENT(IFAIL).AND..NOT.PRESENT(Z) )THEN; LINFO = -12
   ELSE IF( N > 0 )THEN
      IF( PRESENT(Z) ) THEN; LJOBZ = 'V'
        IF( PRESENT(IFAIL) )THEN; LIFAIL => IFAIL
        ELSE; ALLOCATE( LIFAIL(N), STAT=ISTAT ); END IF
        ELSE; LJOBZ = 'N'; LIFAIL => ISTAT1; ENDIF
!     .. DETERMINE THE WORKSPACE
       IF( ISTAT == 0 ) THEN
         ALLOCATE(IWORK(MAX(1,5*N)), WORK(MAX(1,8*N)), STAT=ISTAT)
         IF( ISTAT == 0 )THEN
           IF( PRESENT(ABSTOL) )THEN; LABSTOL = ABSTOL
           ELSE; LABSTOL = 2*LAMCH_F77('Safe minimum'); ENDIF
             IF (PRESENT (Z)) THEN
                  CALL SPGVX_F77( LITYPE, LJOBZ, LRANGE, LUPLO, N, AP, BP, LVL, &
     &              LVU, LIL, LIU, LABSTOL, LM, W, Z, S1Z, WORK, IWORK, LIFAIL, &
     &              LINFO )
             ELSE
	          CALL SPGVX_F77( LITYPE, LJOBZ, LRANGE, LUPLO, N, AP, BP, LVL, &
     &              LVU, LIL, LIU, LABSTOL, LM, W, LLZ, S1Z, WORK, IWORK, LIFAIL, &
     &              LINFO )
             ENDIF
             IF( PRESENT(M) ) M = LM
           ELSE; LINFO = -100; END IF
           END IF
           IF( PRESENT(Z) .AND. .NOT.PRESENT(IFAIL) ) DEALLOCATE( LIFAIL, STAT=ISTAT1(1) )
           DEALLOCATE(IWORK, WORK, STAT=ISTAT1(1))
         END IF
         CALL ERINFO(LINFO,SRNAME,INFO,ISTAT)
END SUBROUTINE SSPGVX_F95