!*> \brief <b> DSYEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices</b>
!*
!*  =========== DOCUMENTATION ===========
!*
!* Online html documentation available at
!*            http://www.netlib.org/lapack/explore-html/
!*
!*> \htmlonly
!*> Download DSYEV + dependencies
!*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyev.f">
!*> [TGZ]</a>
!*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyev.f">
!*> [ZIP]</a>
!*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyev.f">
!*> [TXT]</a>
!*> \endhtmlonly
!*
!*  Definition:
!*  ===========
!*
!*       SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
!*
!*       .. Scalar Arguments ..
!*       CHARACTER          JOBZ, UPLO
!*       INTEGER            INFO, LDA, LWORK, N
!*       ..
!*       .. Array Arguments ..
!*       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * )
!*       ..
!*
!*
!*> \par Purpose:
!*  =============
!*>
!*> \verbatim
!*>
!*> DSYEV computes all eigenvalues and, optionally, eigenvectors of a
!*> real symmetric matrix A.
!*> \endverbatim
!*
!*  Arguments:
!*  ==========
!*
!*> \param[in] JOBZ
!*> \verbatim
!*>          JOBZ is CHARACTER*1
!*>          = 'N':  Compute eigenvalues only;
!*>          = 'V':  Compute eigenvalues and eigenvectors.
!*> \endverbatim
!*>
!*> \param[in] UPLO
!*> \verbatim
!*>          UPLO is CHARACTER*1
!*>          = 'U':  Upper triangle of A is stored;
!*>          = 'L':  Lower triangle of A is stored.
!*> \endverbatim
!*>
!*> \param[in] N
!*> \verbatim
!*>          N is INTEGER
!*>          The order of the matrix A.  N >= 0.
!*> \endverbatim
!*>
!*> \param[in,out] A
!*> \verbatim
!*>          A is DOUBLE PRECISION array, dimension (LDA, N)
!*>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!*>          leading N-by-N upper triangular part of A contains the
!*>          upper triangular part of the matrix A.  If UPLO = 'L',
!*>          the leading N-by-N lower triangular part of A contains
!*>          the lower triangular part of the matrix A.
!*>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
!*>          orthonormal eigenvectors of the matrix A.
!*>          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
!*>          or the upper triangle (if UPLO='U') of A, including the
!*>          diagonal, is destroyed.
!*> \endverbatim
!*>
!*> \param[in] LDA
!*> \verbatim
!*>          LDA is INTEGER
!*>          The leading dimension of the array A.  LDA >= max(1,N).
!*> \endverbatim
!*>
!*> \param[out] W
!*> \verbatim
!*>          W is DOUBLE PRECISION array, dimension (N)
!*>          If INFO = 0, the eigenvalues in ascending order.
!*> \endverbatim
!*>
!*> \param[out] WORK
!*> \verbatim
!*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!*> \endverbatim
!*>
!*> \param[in] LWORK
!*> \verbatim
!*>          LWORK is INTEGER
!*>          The length of the array WORK.  LWORK >= max(1,3*N-1).
!*>          For optimal efficiency, LWORK >= (NB+2)*N,
!*>          where NB is the blocksize for DSYTRD returned by ILAENV.
!*>
!*>          If LWORK = -1, then a workspace query is assumed; the routine
!*>          only calculates the optimal size of the WORK array, returns
!*>          this value as the first entry of the WORK array, and no error
!*>          message related to LWORK is issued by XERBLA.
!*> \endverbatim
!*>
!*> \param[out] INFO
!*> \verbatim
!*>          INFO is INTEGER
!*>          = 0:  successful exit
!*>          < 0:  if INFO = -i, the i-th argument had an illegal value
!*>          > 0:  if INFO = i, the algorithm failed to converge; i
!*>                off-diagonal elements of an intermediate tridiagonal
!*>                form did not converge to zero.
!*> \endverbatim
!*
!*  Authors:
!*  ========
!*
!*> \author Univ. of Tennessee
!*> \author Univ. of California Berkeley
!*> \author Univ. of Colorado Denver
!*> \author NAG Ltd.
!*
!*> \date December 2016
!*
!*> \ingroup doubleSYeigen
!*
module DSYEV_PF
contains
!*  =====================================================================
   SUBROUTINE DSYEV_PF(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO)
!*
!*  -- LAPACK driver routine (version 3.7.0) --
!*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!*     December 2016
!*
!*     .. Scalar Arguments ..
      CHARACTER JOBZ, UPLO
      INTEGER INFO, LDA, LWORK, N
!*     ..
!*     .. Array Arguments ..
      DOUBLE PRECISION A(LDA, *), W(*), WORK(*)
!*     ..
!*
!*  =====================================================================
!*
!*     .. Parameters ..
      DOUBLE PRECISION ZERO, ONE
      PARAMETER(ZERO=0.0D0, ONE=1.0D0)
!*     ..
!*     .. Local Scalars ..
      LOGICAL LOWER, LQUERY, WANTZ
      INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE, &
         LLWORK, LWKOPT, NB
      DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, &
         SMLNUM
!*     ..
!*     .. External Functions ..
      LOGICAL LSAME
      INTEGER ILAENV
      DOUBLE PRECISION DLAMCH, DLANSY
      EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY
!*     ..
!*     .. External Subroutines ..
      EXTERNAL DLASCL, DORGTR, DSCAL, DSTEQR, DSTERF, DSYTRD, &
         XERBLA
!*     ..
!*     .. Intrinsic Functions ..
      INTRINSIC MAX, SQRT
!*     ..
!*     .. Executable Statements ..
!*
!*     Test the input parameters.
!*
      WANTZ = LSAME(JOBZ, 'V')
      LOWER = LSAME(UPLO, 'L')
      LQUERY = (LWORK .EQ. -1)
!*
      INFO = 0
      IF (.NOT. (WANTZ .OR. LSAME(JOBZ, 'N'))) THEN
         INFO = -1
      ELSE IF (.NOT. (LOWER .OR. LSAME(UPLO, 'U'))) THEN
         INFO = -2
      ELSE IF (N .LT. 0) THEN
         INFO = -3
      ELSE IF (LDA .LT. MAX(1, N)) THEN
         INFO = -5
      END IF
!*
      IF (INFO .EQ. 0) THEN
         NB = ILAENV(1, 'DSYTRD', UPLO, N, -1, -1, -1)
         LWKOPT = MAX(1, (NB + 2)*N)
         WORK(1) = LWKOPT
!*
         IF (LWORK .LT. MAX(1, 3*N - 1) .AND. .NOT. LQUERY) &
            INFO = -8
      END IF
!*
      IF (INFO .NE. 0) THEN
         CALL XERBLA('DSYEV ', -INFO)
         RETURN
      ELSE IF (LQUERY) THEN
         RETURN
      END IF
!*
!*     Quick return if possible
!*
      IF (N .EQ. 0) THEN
         RETURN
      END IF
!*
      IF (N .EQ. 1) THEN
         W(1) = A(1, 1)
         WORK(1) = 2
         IF (WANTZ) &
            A(1, 1) = ONE
         RETURN
      END IF
!*
!*     Get machine constants.
!*
      SAFMIN = DLAMCH('Safe minimum')
      EPS = DLAMCH('Precision')
      SMLNUM = SAFMIN/EPS
      BIGNUM = ONE/SMLNUM
      RMIN = SQRT(SMLNUM)
      RMAX = SQRT(BIGNUM)
!*
!*     Scale matrix to allowable range, if necessary.
!*
      ANRM = DLANSY('M', UPLO, N, A, LDA, WORK)
      ISCALE = 0
      IF (ANRM .GT. ZERO .AND. ANRM .LT. RMIN) THEN
         ISCALE = 1
         SIGMA = RMIN/ANRM
      ELSE IF (ANRM .GT. RMAX) THEN
         ISCALE = 1
         SIGMA = RMAX/ANRM
      END IF
      IF (ISCALE .EQ. 1) &
         CALL DLASCL(UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO)
!*
!*     Call DSYTRD to reduce symmetric matrix to tridiagonal form.
!*
      INDE = 1
      INDTAU = INDE + N
      INDWRK = INDTAU + N
      LLWORK = LWORK - INDWRK + 1
      CALL DSYTRD(UPLO, N, A, LDA, W, WORK(INDE), WORK(INDTAU), &
                  WORK(INDWRK), LLWORK, IINFO)
!*
!*     For eigenvalues only, call DSTERF.  For eigenvectors, first call
!*     DORGTR to generate the orthogonal matrix, then call DSTEQR.
!*
      IF (.NOT. WANTZ) THEN
         CALL DSTERF(N, W, WORK(INDE), INFO)
      ELSE
         CALL DORGTR(UPLO, N, A, LDA, WORK(INDTAU), WORK(INDWRK), &
                     LLWORK, IINFO)
         CALL DSTEQR(JOBZ, N, W, WORK(INDE), A, LDA, WORK(INDTAU), &
                     INFO)
      END IF
!*
!*     If matrix was scaled, then rescale eigenvalues appropriately.
!*
      IF (ISCALE .EQ. 1) THEN
         IF (INFO .EQ. 0) THEN
            IMAX = N
         ELSE
            IMAX = INFO - 1
         END IF
         CALL DSCAL(IMAX, ONE/SIGMA, W, 1)
      END IF
!*
!*     Set WORK(1) to optimal workspace size.
!*
      WORK(1) = LWKOPT
!*
      RETURN
!*
!*     End of DSYEV
!*
   END SUBROUTINE
end module DSYEV_PF