Provided by: liblapack-doc_3.9.0-1build1_all bug

NAME

       complex16HEeigen

SYNOPSIS

   Functions
       subroutine zheev (JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO)
            ZHEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors
           for HE matrices
       subroutine zheev_2stage (JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO)
            ZHEEV_2STAGE computes the eigenvalues and, optionally, the left and/or right
           eigenvectors for HE matrices
       subroutine zheevd (JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK,
           INFO)
            ZHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors
           for HE matrices
       subroutine zheevd_2stage (JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, LRWORK, IWORK,
           LIWORK, INFO)
            ZHEEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right
           eigenvectors for HE matrices
       subroutine zheevr (JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ,
           ISUPPZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
            ZHEEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors
           for HE matrices
       subroutine zheevr_2stage (JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z,
           LDZ, ISUPPZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
            ZHEEVR_2STAGE computes the eigenvalues and, optionally, the left and/or right
           eigenvectors for HE matrices
       subroutine zheevx (JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ,
           WORK, LWORK, RWORK, IWORK, IFAIL, INFO)
            ZHEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors
           for HE matrices
       subroutine zheevx_2stage (JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z,
           LDZ, WORK, LWORK, RWORK, IWORK, IFAIL, INFO)
            ZHEEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right
           eigenvectors for HE matrices
       subroutine zhegv (ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, RWORK, INFO)
           ZHEGV
       subroutine zhegv_2stage (ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, RWORK,
           INFO)
           ZHEGV_2STAGE
       subroutine zhegvd (ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, RWORK, LRWORK,
           IWORK, LIWORK, INFO)
           ZHEGVD
       subroutine zhegvx (ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, VL, VU, IL, IU, ABSTOL, M,
           W, Z, LDZ, WORK, LWORK, RWORK, IWORK, IFAIL, INFO)
           ZHEGVX

Detailed Description

       This is the group of complex16 eigenvalue driver functions for HE matrices

Function Documentation

   subroutine zheev (character JOBZ, character UPLO, integer N, complex*16, dimension( lda, * )
       A, integer LDA, double precision, dimension( * ) W, complex*16, dimension( * ) WORK,
       integer LWORK, double precision, dimension( * ) RWORK, integer INFO)
        ZHEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE
       matrices

       Purpose:

            ZHEEV computes all eigenvalues and, optionally, eigenvectors of a
            complex Hermitian matrix A.

       Parameters
           JOBZ

                     JOBZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only;
                     = 'V':  Compute eigenvalues and eigenvectors.

           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  Upper triangle of A is stored;
                     = 'L':  Lower triangle of A is stored.

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, the Hermitian 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.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,N).

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     If INFO = 0, the eigenvalues in ascending order.

           WORK

                     WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The length of the array WORK.  LWORK >= max(1,2*N-1).
                     For optimal efficiency, LWORK >= (NB+1)*N,
                     where NB is the blocksize for ZHETRD 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.

           RWORK

                     RWORK is DOUBLE PRECISION array, dimension (max(1, 3*N-2))

           INFO

                     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.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date
           December 2016

   subroutine zheev_2stage (character JOBZ, character UPLO, integer N, complex*16, dimension(
       lda, * ) A, integer LDA, double precision, dimension( * ) W, complex*16, dimension( * )
       WORK, integer LWORK, double precision, dimension( * ) RWORK, integer INFO)
        ZHEEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors
       for HE matrices

       Purpose:

            ZHEEV_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
            complex Hermitian matrix A using the 2stage technique for
            the reduction to tridiagonal.

       Parameters
           JOBZ

                     JOBZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only;
                     = 'V':  Compute eigenvalues and eigenvectors.
                             Not available in this release.

           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  Upper triangle of A is stored;
                     = 'L':  Lower triangle of A is stored.

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, the Hermitian 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.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,N).

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     If INFO = 0, the eigenvalues in ascending order.

           WORK

                     WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The length of the array WORK. LWORK >= 1, when N <= 1;
                     otherwise
                     If JOBZ = 'N' and N > 1, LWORK must be queried.
                                              LWORK = MAX(1, dimension) where
                                              dimension = max(stage1,stage2) + (KD+1)*N + N
                                                        = N*KD + N*max(KD+1,FACTOPTNB)
                                                          + max(2*KD*KD, KD*NTHREADS)
                                                          + (KD+1)*N + N
                                              where KD is the blocking size of the reduction,
                                              FACTOPTNB is the blocking used by the QR or LQ
                                              algorithm, usually FACTOPTNB=128 is a good choice
                                              NTHREADS is the number of threads used when
                                              openMP compilation is enabled, otherwise =1.
                     If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available

                     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.

           RWORK

                     RWORK is DOUBLE PRECISION array, dimension (max(1, 3*N-2))

           INFO

                     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.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date
           November 2017

       Further Details:

             All details about the 2stage techniques are available in:

             Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
             Parallel reduction to condensed forms for symmetric eigenvalue problems
             using aggregated fine-grained and memory-aware kernels. In Proceedings
             of 2011 International Conference for High Performance Computing,
             Networking, Storage and Analysis (SC '11), New York, NY, USA,
             Article 8 , 11 pages.
             http://doi.acm.org/10.1145/2063384.2063394

             A. Haidar, J. Kurzak, P. Luszczek, 2013.
             An improved parallel singular value algorithm and its implementation
             for multicore hardware, In Proceedings of 2013 International Conference
             for High Performance Computing, Networking, Storage and Analysis (SC '13).
             Denver, Colorado, USA, 2013.
             Article 90, 12 pages.
             http://doi.acm.org/10.1145/2503210.2503292

             A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
             A novel hybrid CPU-GPU generalized eigensolver for electronic structure
             calculations based on fine-grained memory aware tasks.
             International Journal of High Performance Computing Applications.
             Volume 28 Issue 2, Pages 196-209, May 2014.
             http://hpc.sagepub.com/content/28/2/196

   subroutine zheevd (character JOBZ, character UPLO, integer N, complex*16, dimension( lda, * )
       A, integer LDA, double precision, dimension( * ) W, complex*16, dimension( * ) WORK,
       integer LWORK, double precision, dimension( * ) RWORK, integer LRWORK, integer, dimension(
       * ) IWORK, integer LIWORK, integer INFO)
        ZHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for
       HE matrices

       Purpose:

            ZHEEVD computes all eigenvalues and, optionally, eigenvectors of a
            complex Hermitian matrix A.  If eigenvectors are desired, it uses a
            divide and conquer algorithm.

            The divide and conquer algorithm makes very mild assumptions about
            floating point arithmetic. It will work on machines with a guard
            digit in add/subtract, or on those binary machines without guard
            digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
            Cray-2. It could conceivably fail on hexadecimal or decimal machines
            without guard digits, but we know of none.

       Parameters
           JOBZ

                     JOBZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only;
                     = 'V':  Compute eigenvalues and eigenvectors.

           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  Upper triangle of A is stored;
                     = 'L':  Lower triangle of A is stored.

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, the Hermitian 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.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,N).

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     If INFO = 0, the eigenvalues in ascending order.

           WORK

                     WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The length of the array WORK.
                     If N <= 1,                LWORK must be at least 1.
                     If JOBZ  = 'N' and N > 1, LWORK must be at least N + 1.
                     If JOBZ  = 'V' and N > 1, LWORK must be at least 2*N + N**2.

                     If LWORK = -1, then a workspace query is assumed; the routine
                     only calculates the optimal sizes of the WORK, RWORK and
                     IWORK arrays, returns these values as the first entries of
                     the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           RWORK

                     RWORK is DOUBLE PRECISION array,
                                                    dimension (LRWORK)
                     On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.

           LRWORK

                     LRWORK is INTEGER
                     The dimension of the array RWORK.
                     If N <= 1,                LRWORK must be at least 1.
                     If JOBZ  = 'N' and N > 1, LRWORK must be at least N.
                     If JOBZ  = 'V' and N > 1, LRWORK must be at least
                                    1 + 5*N + 2*N**2.

                     If LRWORK = -1, then a workspace query is assumed; the
                     routine only calculates the optimal sizes of the WORK, RWORK
                     and IWORK arrays, returns these values as the first entries
                     of the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           IWORK

                     IWORK is INTEGER array, dimension (MAX(1,LIWORK))
                     On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.

           LIWORK

                     LIWORK is INTEGER
                     The dimension of the array IWORK.
                     If N <= 1,                LIWORK must be at least 1.
                     If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
                     If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.

                     If LIWORK = -1, then a workspace query is assumed; the
                     routine only calculates the optimal sizes of the WORK, RWORK
                     and IWORK arrays, returns these values as the first entries
                     of the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  if INFO = i and JOBZ = 'N', then the algorithm failed
                           to converge; i off-diagonal elements of an intermediate
                           tridiagonal form did not converge to zero;
                           if INFO = i and JOBZ = 'V', then the algorithm failed
                           to compute an eigenvalue while working on the submatrix
                           lying in rows and columns INFO/(N+1) through
                           mod(INFO,N+1).

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date
           December 2016

       Further Details:
           Modified description of INFO. Sven, 16 Feb 05.

       Contributors:
           Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

   subroutine zheevd_2stage (character JOBZ, character UPLO, integer N, complex*16, dimension(
       lda, * ) A, integer LDA, double precision, dimension( * ) W, complex*16, dimension( * )
       WORK, integer LWORK, double precision, dimension( * ) RWORK, integer LRWORK, integer,
       dimension( * ) IWORK, integer LIWORK, integer INFO)
        ZHEEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right
       eigenvectors for HE matrices

       Purpose:

            ZHEEVD_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
            complex Hermitian matrix A using the 2stage technique for
            the reduction to tridiagonal.  If eigenvectors are desired, it uses a
            divide and conquer algorithm.

            The divide and conquer algorithm makes very mild assumptions about
            floating point arithmetic. It will work on machines with a guard
            digit in add/subtract, or on those binary machines without guard
            digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
            Cray-2. It could conceivably fail on hexadecimal or decimal machines
            without guard digits, but we know of none.

       Parameters
           JOBZ

                     JOBZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only;
                     = 'V':  Compute eigenvalues and eigenvectors.
                             Not available in this release.

           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  Upper triangle of A is stored;
                     = 'L':  Lower triangle of A is stored.

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, the Hermitian 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.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,N).

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     If INFO = 0, the eigenvalues in ascending order.

           WORK

                     WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The dimension of the array WORK.
                     If N <= 1,               LWORK must be at least 1.
                     If JOBZ = 'N' and N > 1, LWORK must be queried.
                                              LWORK = MAX(1, dimension) where
                                              dimension = max(stage1,stage2) + (KD+1)*N + N+1
                                                        = N*KD + N*max(KD+1,FACTOPTNB)
                                                          + max(2*KD*KD, KD*NTHREADS)
                                                          + (KD+1)*N + N+1
                                              where KD is the blocking size of the reduction,
                                              FACTOPTNB is the blocking used by the QR or LQ
                                              algorithm, usually FACTOPTNB=128 is a good choice
                                              NTHREADS is the number of threads used when
                                              openMP compilation is enabled, otherwise =1.
                     If JOBZ = 'V' and N > 1, LWORK must be at least 2*N + N**2

                     If LWORK = -1, then a workspace query is assumed; the routine
                     only calculates the optimal sizes of the WORK, RWORK and
                     IWORK arrays, returns these values as the first entries of
                     the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           RWORK

                     RWORK is DOUBLE PRECISION array,
                                                    dimension (LRWORK)
                     On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.

           LRWORK

                     LRWORK is INTEGER
                     The dimension of the array RWORK.
                     If N <= 1,                LRWORK must be at least 1.
                     If JOBZ  = 'N' and N > 1, LRWORK must be at least N.
                     If JOBZ  = 'V' and N > 1, LRWORK must be at least
                                    1 + 5*N + 2*N**2.

                     If LRWORK = -1, then a workspace query is assumed; the
                     routine only calculates the optimal sizes of the WORK, RWORK
                     and IWORK arrays, returns these values as the first entries
                     of the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           IWORK

                     IWORK is INTEGER array, dimension (MAX(1,LIWORK))
                     On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.

           LIWORK

                     LIWORK is INTEGER
                     The dimension of the array IWORK.
                     If N <= 1,                LIWORK must be at least 1.
                     If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
                     If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.

                     If LIWORK = -1, then a workspace query is assumed; the
                     routine only calculates the optimal sizes of the WORK, RWORK
                     and IWORK arrays, returns these values as the first entries
                     of the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  if INFO = i and JOBZ = 'N', then the algorithm failed
                           to converge; i off-diagonal elements of an intermediate
                           tridiagonal form did not converge to zero;
                           if INFO = i and JOBZ = 'V', then the algorithm failed
                           to compute an eigenvalue while working on the submatrix
                           lying in rows and columns INFO/(N+1) through
                           mod(INFO,N+1).

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date
           November 2017

       Further Details:
           Modified description of INFO. Sven, 16 Feb 05.

       Contributors:
           Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

       Further Details:

             All details about the 2stage techniques are available in:

             Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
             Parallel reduction to condensed forms for symmetric eigenvalue problems
             using aggregated fine-grained and memory-aware kernels. In Proceedings
             of 2011 International Conference for High Performance Computing,
             Networking, Storage and Analysis (SC '11), New York, NY, USA,
             Article 8 , 11 pages.
             http://doi.acm.org/10.1145/2063384.2063394

             A. Haidar, J. Kurzak, P. Luszczek, 2013.
             An improved parallel singular value algorithm and its implementation
             for multicore hardware, In Proceedings of 2013 International Conference
             for High Performance Computing, Networking, Storage and Analysis (SC '13).
             Denver, Colorado, USA, 2013.
             Article 90, 12 pages.
             http://doi.acm.org/10.1145/2503210.2503292

             A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
             A novel hybrid CPU-GPU generalized eigensolver for electronic structure
             calculations based on fine-grained memory aware tasks.
             International Journal of High Performance Computing Applications.
             Volume 28 Issue 2, Pages 196-209, May 2014.
             http://hpc.sagepub.com/content/28/2/196

   subroutine zheevr (character JOBZ, character RANGE, character UPLO, integer N, complex*16,
       dimension( lda, * ) A, integer LDA, double precision VL, double precision VU, integer IL,
       integer IU, double precision ABSTOL, integer M, double precision, dimension( * ) W,
       complex*16, dimension( ldz, * ) Z, integer LDZ, integer, dimension( * ) ISUPPZ,
       complex*16, dimension( * ) WORK, integer LWORK, double precision, dimension( * ) RWORK,
       integer LRWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO)
        ZHEEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for
       HE matrices

       Purpose:

            ZHEEVR computes selected eigenvalues and, optionally, eigenvectors
            of a 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.

            ZHEEVR first reduces the matrix A to tridiagonal form T with a call
            to ZHETRD.  Then, whenever possible, ZHEEVR calls ZSTEMR to compute
            eigenspectrum using Relatively Robust Representations.  ZSTEMR
            computes eigenvalues by the dqds algorithm, while orthogonal
            eigenvectors are computed from various "good" L D L^T representations
            (also known as Relatively Robust Representations). Gram-Schmidt
            orthogonalization is avoided as far as possible. More specifically,
            the various steps of the algorithm are as follows.

            For each unreduced block (submatrix) of T,
               (a) Compute T - sigma I  = L D L^T, so that L and D
                   define all the wanted eigenvalues to high relative accuracy.
                   This means that small relative changes in the entries of D and L
                   cause only small relative changes in the eigenvalues and
                   eigenvectors. The standard (unfactored) representation of the
                   tridiagonal matrix T does not have this property in general.
               (b) Compute the eigenvalues to suitable accuracy.
                   If the eigenvectors are desired, the algorithm attains full
                   accuracy of the computed eigenvalues only right before
                   the corresponding vectors have to be computed, see steps c) and d).
               (c) For each cluster of close eigenvalues, select a new
                   shift close to the cluster, find a new factorization, and refine
                   the shifted eigenvalues to suitable accuracy.
               (d) For each eigenvalue with a large enough relative separation compute
                   the corresponding eigenvector by forming a rank revealing twisted
                   factorization. Go back to (c) for any clusters that remain.

            The desired accuracy of the output can be specified by the input
            parameter ABSTOL.

            For more details, see DSTEMR's documentation and:
            - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
              to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
              Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
            - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
              Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
              2004.  Also LAPACK Working Note 154.
            - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
              tridiagonal eigenvalue/eigenvector problem",
              Computer Science Division Technical Report No. UCB/CSD-97-971,
              UC Berkeley, May 1997.

            Note 1 : ZHEEVR calls ZSTEMR when the full spectrum is requested
            on machines which conform to the ieee-754 floating point standard.
            ZHEEVR calls DSTEBZ and ZSTEIN on non-ieee machines and
            when partial spectrum requests are made.

            Normal execution of ZSTEMR may create NaNs and infinities and
            hence may abort due to a floating point exception in environments
            which do not handle NaNs and infinities in the ieee standard default
            manner.

       Parameters
           JOBZ

                     JOBZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only;
                     = 'V':  Compute eigenvalues and eigenvectors.

           RANGE

                     RANGE is CHARACTER*1
                     = 'A': all eigenvalues will be found.
                     = 'V': all eigenvalues in the half-open interval (VL,VU]
                            will be found.
                     = 'I': the IL-th through IU-th eigenvalues will be found.
                     For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
                     ZSTEIN are called

           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  Upper triangle of A is stored;
                     = 'L':  Lower triangle of A is stored.

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, the Hermitian 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, the lower triangle (if UPLO='L') or the upper
                     triangle (if UPLO='U') of A, including the diagonal, is
                     destroyed.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,N).

           VL

                     VL is DOUBLE PRECISION
                     If RANGE='V', the lower bound of the interval to
                     be searched for eigenvalues. VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           VU

                     VU is DOUBLE PRECISION
                     If RANGE='V', the upper bound of the interval to
                     be searched for eigenvalues. VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           IL

                     IL is INTEGER
                     If RANGE='I', the index of the
                     smallest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           IU

                     IU is INTEGER
                     If RANGE='I', the index of the
                     largest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           ABSTOL

                     ABSTOL is DOUBLE PRECISION
                     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 + EPS *   max( |a|,|b| ) ,

                     where EPS is the machine precision.  If ABSTOL is less than
                     or equal to zero, then  EPS*|T|  will be used in its place,
                     where |T| is the 1-norm of the tridiagonal matrix obtained
                     by reducing A to tridiagonal form.

                     See "Computing Small Singular Values of Bidiagonal Matrices
                     with Guaranteed High Relative Accuracy," by Demmel and
                     Kahan, LAPACK Working Note #3.

                     If high relative accuracy is important, set ABSTOL to
                     DLAMCH( 'Safe minimum' ).  Doing so will guarantee that
                     eigenvalues are computed to high relative accuracy when
                     possible in future releases.  The current code does not
                     make any guarantees about high relative accuracy, but
                     future releases will. See J. Barlow and J. Demmel,
                     "Computing Accurate Eigensystems of Scaled Diagonally
                     Dominant Matrices", LAPACK Working Note #7, for a discussion
                     of which matrices define their eigenvalues to high relative
                     accuracy.

           M

                     M is INTEGER
                     The total number of eigenvalues found.  0 <= M <= N.
                     If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     The first M elements contain the selected eigenvalues in
                     ascending order.

           Z

                     Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
                     If JOBZ = 'V', then if INFO = 0, the first M columns of Z
                     contain the orthonormal eigenvectors of the matrix A
                     corresponding to the selected eigenvalues, with the i-th
                     column of Z holding the eigenvector associated with W(i).
                     If JOBZ = 'N', then Z is not referenced.
                     Note: the user must ensure that at least max(1,M) columns are
                     supplied in the array Z; if RANGE = 'V', the exact value of M
                     is not known in advance and an upper bound must be used.

           LDZ

                     LDZ is INTEGER
                     The leading dimension of the array Z.  LDZ >= 1, and if
                     JOBZ = 'V', LDZ >= max(1,N).

           ISUPPZ

                     ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
                     The support of the eigenvectors in Z, i.e., the indices
                     indicating the nonzero elements in Z. The i-th eigenvector
                     is nonzero only in elements ISUPPZ( 2*i-1 ) through
                     ISUPPZ( 2*i ). This is an output of ZSTEMR (tridiagonal
                     matrix). The support of the eigenvectors of A is typically
                     1:N because of the unitary transformations applied by ZUNMTR.
                     Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1

           WORK

                     WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The length of the array WORK.  LWORK >= max(1,2*N).
                     For optimal efficiency, LWORK >= (NB+1)*N,
                     where NB is the max of the blocksize for ZHETRD and for
                     ZUNMTR as returned by ILAENV.

                     If LWORK = -1, then a workspace query is assumed; the routine
                     only calculates the optimal sizes of the WORK, RWORK and
                     IWORK arrays, returns these values as the first entries of
                     the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           RWORK

                     RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
                     On exit, if INFO = 0, RWORK(1) returns the optimal
                     (and minimal) LRWORK.

           LRWORK

                     LRWORK is INTEGER
                     The length of the array RWORK.  LRWORK >= max(1,24*N).

                     If LRWORK = -1, then a workspace query is assumed; the
                     routine only calculates the optimal sizes of the WORK, RWORK
                     and IWORK arrays, returns these values as the first entries
                     of the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           IWORK

                     IWORK is INTEGER array, dimension (MAX(1,LIWORK))
                     On exit, if INFO = 0, IWORK(1) returns the optimal
                     (and minimal) LIWORK.

           LIWORK

                     LIWORK is INTEGER
                     The dimension of the array IWORK.  LIWORK >= max(1,10*N).

                     If LIWORK = -1, then a workspace query is assumed; the
                     routine only calculates the optimal sizes of the WORK, RWORK
                     and IWORK arrays, returns these values as the first entries
                     of the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  Internal error

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date
           June 2016

       Contributors:
           Inderjit Dhillon, IBM Almaden, USA
            Osni Marques, LBNL/NERSC, USA
            Ken Stanley, Computer Science Division, University of California at Berkeley, USA
            Jason Riedy, Computer Science Division, University of California at Berkeley, USA

   subroutine zheevr_2stage (character JOBZ, character RANGE, character UPLO, integer N,
       complex*16, dimension( lda, * ) A, integer LDA, double precision VL, double precision VU,
       integer IL, integer IU, double precision ABSTOL, integer M, double precision, dimension( *
       ) W, complex*16, dimension( ldz, * ) Z, integer LDZ, integer, dimension( * ) ISUPPZ,
       complex*16, dimension( * ) WORK, integer LWORK, double precision, dimension( * ) RWORK,
       integer LRWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO)
        ZHEEVR_2STAGE computes the eigenvalues and, optionally, the left and/or right
       eigenvectors for HE matrices

       Purpose:

            ZHEEVR_2STAGE computes selected eigenvalues and, optionally, eigenvectors
            of a complex Hermitian matrix A using the 2stage technique for
            the reduction to tridiagonal.  Eigenvalues and eigenvectors can
            be selected by specifying either a range of values or a range of
            indices for the desired eigenvalues.

            ZHEEVR_2STAGE first reduces the matrix A to tridiagonal form T with a call
            to ZHETRD.  Then, whenever possible, ZHEEVR_2STAGE calls ZSTEMR to compute
            eigenspectrum using Relatively Robust Representations.  ZSTEMR
            computes eigenvalues by the dqds algorithm, while orthogonal
            eigenvectors are computed from various "good" L D L^T representations
            (also known as Relatively Robust Representations). Gram-Schmidt
            orthogonalization is avoided as far as possible. More specifically,
            the various steps of the algorithm are as follows.

            For each unreduced block (submatrix) of T,
               (a) Compute T - sigma I  = L D L^T, so that L and D
                   define all the wanted eigenvalues to high relative accuracy.
                   This means that small relative changes in the entries of D and L
                   cause only small relative changes in the eigenvalues and
                   eigenvectors. The standard (unfactored) representation of the
                   tridiagonal matrix T does not have this property in general.
               (b) Compute the eigenvalues to suitable accuracy.
                   If the eigenvectors are desired, the algorithm attains full
                   accuracy of the computed eigenvalues only right before
                   the corresponding vectors have to be computed, see steps c) and d).
               (c) For each cluster of close eigenvalues, select a new
                   shift close to the cluster, find a new factorization, and refine
                   the shifted eigenvalues to suitable accuracy.
               (d) For each eigenvalue with a large enough relative separation compute
                   the corresponding eigenvector by forming a rank revealing twisted
                   factorization. Go back to (c) for any clusters that remain.

            The desired accuracy of the output can be specified by the input
            parameter ABSTOL.

            For more details, see DSTEMR's documentation and:
            - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
              to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
              Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
            - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
              Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
              2004.  Also LAPACK Working Note 154.
            - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
              tridiagonal eigenvalue/eigenvector problem",
              Computer Science Division Technical Report No. UCB/CSD-97-971,
              UC Berkeley, May 1997.

            Note 1 : ZHEEVR_2STAGE calls ZSTEMR when the full spectrum is requested
            on machines which conform to the ieee-754 floating point standard.
            ZHEEVR_2STAGE calls DSTEBZ and ZSTEIN on non-ieee machines and
            when partial spectrum requests are made.

            Normal execution of ZSTEMR may create NaNs and infinities and
            hence may abort due to a floating point exception in environments
            which do not handle NaNs and infinities in the ieee standard default
            manner.

       Parameters
           JOBZ

                     JOBZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only;
                     = 'V':  Compute eigenvalues and eigenvectors.
                             Not available in this release.

           RANGE

                     RANGE is CHARACTER*1
                     = 'A': all eigenvalues will be found.
                     = 'V': all eigenvalues in the half-open interval (VL,VU]
                            will be found.
                     = 'I': the IL-th through IU-th eigenvalues will be found.
                     For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
                     ZSTEIN are called

           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  Upper triangle of A is stored;
                     = 'L':  Lower triangle of A is stored.

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, the Hermitian 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, the lower triangle (if UPLO='L') or the upper
                     triangle (if UPLO='U') of A, including the diagonal, is
                     destroyed.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,N).

           VL

                     VL is DOUBLE PRECISION
                     If RANGE='V', the lower bound of the interval to
                     be searched for eigenvalues. VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           VU

                     VU is DOUBLE PRECISION
                     If RANGE='V', the upper bound of the interval to
                     be searched for eigenvalues. VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           IL

                     IL is INTEGER
                     If RANGE='I', the index of the
                     smallest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           IU

                     IU is INTEGER
                     If RANGE='I', the index of the
                     largest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           ABSTOL

                     ABSTOL is DOUBLE PRECISION
                     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 + EPS *   max( |a|,|b| ) ,

                     where EPS is the machine precision.  If ABSTOL is less than
                     or equal to zero, then  EPS*|T|  will be used in its place,
                     where |T| is the 1-norm of the tridiagonal matrix obtained
                     by reducing A to tridiagonal form.

                     See "Computing Small Singular Values of Bidiagonal Matrices
                     with Guaranteed High Relative Accuracy," by Demmel and
                     Kahan, LAPACK Working Note #3.

                     If high relative accuracy is important, set ABSTOL to
                     DLAMCH( 'Safe minimum' ).  Doing so will guarantee that
                     eigenvalues are computed to high relative accuracy when
                     possible in future releases.  The current code does not
                     make any guarantees about high relative accuracy, but
                     future releases will. See J. Barlow and J. Demmel,
                     "Computing Accurate Eigensystems of Scaled Diagonally
                     Dominant Matrices", LAPACK Working Note #7, for a discussion
                     of which matrices define their eigenvalues to high relative
                     accuracy.

           M

                     M is INTEGER
                     The total number of eigenvalues found.  0 <= M <= N.
                     If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     The first M elements contain the selected eigenvalues in
                     ascending order.

           Z

                     Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
                     If JOBZ = 'V', then if INFO = 0, the first M columns of Z
                     contain the orthonormal eigenvectors of the matrix A
                     corresponding to the selected eigenvalues, with the i-th
                     column of Z holding the eigenvector associated with W(i).
                     If JOBZ = 'N', then Z is not referenced.
                     Note: the user must ensure that at least max(1,M) columns are
                     supplied in the array Z; if RANGE = 'V', the exact value of M
                     is not known in advance and an upper bound must be used.

           LDZ

                     LDZ is INTEGER
                     The leading dimension of the array Z.  LDZ >= 1, and if
                     JOBZ = 'V', LDZ >= max(1,N).

           ISUPPZ

                     ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
                     The support of the eigenvectors in Z, i.e., the indices
                     indicating the nonzero elements in Z. The i-th eigenvector
                     is nonzero only in elements ISUPPZ( 2*i-1 ) through
                     ISUPPZ( 2*i ). This is an output of ZSTEMR (tridiagonal
                     matrix). The support of the eigenvectors of A is typically
                     1:N because of the unitary transformations applied by ZUNMTR.
                     Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1

           WORK

                     WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The dimension of the array WORK.
                     If JOBZ = 'N' and N > 1, LWORK must be queried.
                                              LWORK = MAX(1, 26*N, dimension) where
                                              dimension = max(stage1,stage2) + (KD+1)*N + N
                                                        = N*KD + N*max(KD+1,FACTOPTNB)
                                                          + max(2*KD*KD, KD*NTHREADS)
                                                          + (KD+1)*N + N
                                              where KD is the blocking size of the reduction,
                                              FACTOPTNB is the blocking used by the QR or LQ
                                              algorithm, usually FACTOPTNB=128 is a good choice
                                              NTHREADS is the number of threads used when
                                              openMP compilation is enabled, otherwise =1.
                     If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available

                     If LWORK = -1, then a workspace query is assumed; the routine
                     only calculates the optimal sizes of the WORK, RWORK and
                     IWORK arrays, returns these values as the first entries of
                     the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           RWORK

                     RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
                     On exit, if INFO = 0, RWORK(1) returns the optimal
                     (and minimal) LRWORK.

           LRWORK

                     LRWORK is INTEGER
                     The length of the array RWORK.  LRWORK >= max(1,24*N).

                     If LRWORK = -1, then a workspace query is assumed; the
                     routine only calculates the optimal sizes of the WORK, RWORK
                     and IWORK arrays, returns these values as the first entries
                     of the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           IWORK

                     IWORK is INTEGER array, dimension (MAX(1,LIWORK))
                     On exit, if INFO = 0, IWORK(1) returns the optimal
                     (and minimal) LIWORK.

           LIWORK

                     LIWORK is INTEGER
                     The dimension of the array IWORK.  LIWORK >= max(1,10*N).

                     If LIWORK = -1, then a workspace query is assumed; the
                     routine only calculates the optimal sizes of the WORK, RWORK
                     and IWORK arrays, returns these values as the first entries
                     of the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  Internal error

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date
           June 2016

       Contributors:

           Inderjit Dhillon, IBM Almaden, USA
           Osni Marques, LBNL/NERSC, USA
           Ken Stanley, Computer Science Division, University of
             California at Berkeley, USA
           Jason Riedy, Computer Science Division, University of
             California at Berkeley, USA

       Further Details:

             All details about the 2stage techniques are available in:

             Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
             Parallel reduction to condensed forms for symmetric eigenvalue problems
             using aggregated fine-grained and memory-aware kernels. In Proceedings
             of 2011 International Conference for High Performance Computing,
             Networking, Storage and Analysis (SC '11), New York, NY, USA,
             Article 8 , 11 pages.
             http://doi.acm.org/10.1145/2063384.2063394

             A. Haidar, J. Kurzak, P. Luszczek, 2013.
             An improved parallel singular value algorithm and its implementation
             for multicore hardware, In Proceedings of 2013 International Conference
             for High Performance Computing, Networking, Storage and Analysis (SC '13).
             Denver, Colorado, USA, 2013.
             Article 90, 12 pages.
             http://doi.acm.org/10.1145/2503210.2503292

             A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
             A novel hybrid CPU-GPU generalized eigensolver for electronic structure
             calculations based on fine-grained memory aware tasks.
             International Journal of High Performance Computing Applications.
             Volume 28 Issue 2, Pages 196-209, May 2014.
             http://hpc.sagepub.com/content/28/2/196

   subroutine zheevx (character JOBZ, character RANGE, character UPLO, integer N, complex*16,
       dimension( lda, * ) A, integer LDA, double precision VL, double precision VU, integer IL,
       integer IU, double precision ABSTOL, integer M, double precision, dimension( * ) W,
       complex*16, dimension( ldz, * ) Z, integer LDZ, complex*16, dimension( * ) WORK, integer
       LWORK, double precision, dimension( * ) RWORK, integer, dimension( * ) IWORK, integer,
       dimension( * ) IFAIL, integer INFO)
        ZHEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for
       HE matrices

       Purpose:

            ZHEEVX computes selected eigenvalues and, optionally, eigenvectors
            of a 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.

       Parameters
           JOBZ

                     JOBZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only;
                     = 'V':  Compute eigenvalues and eigenvectors.

           RANGE

                     RANGE is CHARACTER*1
                     = 'A': all eigenvalues will be found.
                     = 'V': all eigenvalues in the half-open interval (VL,VU]
                            will be found.
                     = 'I': the IL-th through IU-th eigenvalues will be found.

           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  Upper triangle of A is stored;
                     = 'L':  Lower triangle of A is stored.

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, the Hermitian 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, the lower triangle (if UPLO='L') or the upper
                     triangle (if UPLO='U') of A, including the diagonal, is
                     destroyed.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,N).

           VL

                     VL is DOUBLE PRECISION
                     If RANGE='V', the lower bound of the interval to
                     be searched for eigenvalues. VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           VU

                     VU is DOUBLE PRECISION
                     If RANGE='V', the upper bound of the interval to
                     be searched for eigenvalues. VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           IL

                     IL is INTEGER
                     If RANGE='I', the index of the
                     smallest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           IU

                     IU is INTEGER
                     If RANGE='I', the index of the
                     largest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           ABSTOL

                     ABSTOL is DOUBLE PRECISION
                     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 + EPS *   max( |a|,|b| ) ,

                     where EPS is the machine precision.  If ABSTOL is less than
                     or equal to zero, then  EPS*|T|  will be used in its place,
                     where |T| is the 1-norm of the tridiagonal matrix obtained
                     by reducing A to tridiagonal form.

                     Eigenvalues will be computed most accurately when ABSTOL is
                     set to twice the underflow threshold 2*DLAMCH('S'), not zero.
                     If this routine returns with INFO>0, indicating that some
                     eigenvectors did not converge, try setting ABSTOL to
                     2*DLAMCH('S').

                     See "Computing Small Singular Values of Bidiagonal Matrices
                     with Guaranteed High Relative Accuracy," by Demmel and
                     Kahan, LAPACK Working Note #3.

           M

                     M is INTEGER
                     The total number of eigenvalues found.  0 <= M <= N.
                     If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     On normal exit, the first M elements contain the selected
                     eigenvalues in ascending order.

           Z

                     Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
                     If JOBZ = 'V', then if INFO = 0, the first M columns of Z
                     contain the orthonormal eigenvectors of the matrix A
                     corresponding to the selected eigenvalues, with the i-th
                     column of Z holding the eigenvector associated with W(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.
                     If JOBZ = 'N', then Z is not referenced.
                     Note: the user must ensure that at least max(1,M) columns are
                     supplied in the array Z; if RANGE = 'V', the exact value of M
                     is not known in advance and an upper bound must be used.

           LDZ

                     LDZ is INTEGER
                     The leading dimension of the array Z.  LDZ >= 1, and if
                     JOBZ = 'V', LDZ >= max(1,N).

           WORK

                     WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The length of the array WORK.  LWORK >= 1, when N <= 1;
                     otherwise 2*N.
                     For optimal efficiency, LWORK >= (NB+1)*N,
                     where NB is the max of the blocksize for ZHETRD and for
                     ZUNMTR as 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.

           RWORK

                     RWORK is DOUBLE PRECISION array, dimension (7*N)

           IWORK

                     IWORK is INTEGER array, dimension (5*N)

           IFAIL

                     IFAIL is INTEGER array, dimension (N)
                     If JOBZ = 'V', then 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.
                     If JOBZ = 'N', then IFAIL is not referenced.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  if INFO = i, then i eigenvectors failed to converge.
                           Their indices are stored in array IFAIL.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date
           June 2016

   subroutine zheevx_2stage (character JOBZ, character RANGE, character UPLO, integer N,
       complex*16, dimension( lda, * ) A, integer LDA, double precision VL, double precision VU,
       integer IL, integer IU, double precision ABSTOL, integer M, double precision, dimension( *
       ) W, complex*16, dimension( ldz, * ) Z, integer LDZ, complex*16, dimension( * ) WORK,
       integer LWORK, double precision, dimension( * ) RWORK, integer, dimension( * ) IWORK,
       integer, dimension( * ) IFAIL, integer INFO)
        ZHEEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right
       eigenvectors for HE matrices

       Purpose:

            ZHEEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
            of a complex Hermitian matrix A using the 2stage technique for
            the reduction to tridiagonal.  Eigenvalues and eigenvectors can
            be selected by specifying either a range of values or a range of
            indices for the desired eigenvalues.

       Parameters
           JOBZ

                     JOBZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only;
                     = 'V':  Compute eigenvalues and eigenvectors.
                             Not available in this release.

           RANGE

                     RANGE is CHARACTER*1
                     = 'A': all eigenvalues will be found.
                     = 'V': all eigenvalues in the half-open interval (VL,VU]
                            will be found.
                     = 'I': the IL-th through IU-th eigenvalues will be found.

           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  Upper triangle of A is stored;
                     = 'L':  Lower triangle of A is stored.

           N

                     N is INTEGER
                     The order of the matrix A.  N >= 0.

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, the Hermitian 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, the lower triangle (if UPLO='L') or the upper
                     triangle (if UPLO='U') of A, including the diagonal, is
                     destroyed.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,N).

           VL

                     VL is DOUBLE PRECISION
                     If RANGE='V', the lower bound of the interval to
                     be searched for eigenvalues. VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           VU

                     VU is DOUBLE PRECISION
                     If RANGE='V', the upper bound of the interval to
                     be searched for eigenvalues. VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           IL

                     IL is INTEGER
                     If RANGE='I', the index of the
                     smallest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           IU

                     IU is INTEGER
                     If RANGE='I', the index of the
                     largest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           ABSTOL

                     ABSTOL is DOUBLE PRECISION
                     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 + EPS *   max( |a|,|b| ) ,

                     where EPS is the machine precision.  If ABSTOL is less than
                     or equal to zero, then  EPS*|T|  will be used in its place,
                     where |T| is the 1-norm of the tridiagonal matrix obtained
                     by reducing A to tridiagonal form.

                     Eigenvalues will be computed most accurately when ABSTOL is
                     set to twice the underflow threshold 2*DLAMCH('S'), not zero.
                     If this routine returns with INFO>0, indicating that some
                     eigenvectors did not converge, try setting ABSTOL to
                     2*DLAMCH('S').

                     See "Computing Small Singular Values of Bidiagonal Matrices
                     with Guaranteed High Relative Accuracy," by Demmel and
                     Kahan, LAPACK Working Note #3.

           M

                     M is INTEGER
                     The total number of eigenvalues found.  0 <= M <= N.
                     If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     On normal exit, the first M elements contain the selected
                     eigenvalues in ascending order.

           Z

                     Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
                     If JOBZ = 'V', then if INFO = 0, the first M columns of Z
                     contain the orthonormal eigenvectors of the matrix A
                     corresponding to the selected eigenvalues, with the i-th
                     column of Z holding the eigenvector associated with W(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.
                     If JOBZ = 'N', then Z is not referenced.
                     Note: the user must ensure that at least max(1,M) columns are
                     supplied in the array Z; if RANGE = 'V', the exact value of M
                     is not known in advance and an upper bound must be used.

           LDZ

                     LDZ is INTEGER
                     The leading dimension of the array Z.  LDZ >= 1, and if
                     JOBZ = 'V', LDZ >= max(1,N).

           WORK

                     WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The length of the array WORK. LWORK >= 1, when N <= 1;
                     otherwise
                     If JOBZ = 'N' and N > 1, LWORK must be queried.
                                              LWORK = MAX(1, 8*N, dimension) where
                                              dimension = max(stage1,stage2) + (KD+1)*N + N
                                                        = N*KD + N*max(KD+1,FACTOPTNB)
                                                          + max(2*KD*KD, KD*NTHREADS)
                                                          + (KD+1)*N + N
                                              where KD is the blocking size of the reduction,
                                              FACTOPTNB is the blocking used by the QR or LQ
                                              algorithm, usually FACTOPTNB=128 is a good choice
                                              NTHREADS is the number of threads used when
                                              openMP compilation is enabled, otherwise =1.
                     If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available

                     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.

           RWORK

                     RWORK is DOUBLE PRECISION array, dimension (7*N)

           IWORK

                     IWORK is INTEGER array, dimension (5*N)

           IFAIL

                     IFAIL is INTEGER array, dimension (N)
                     If JOBZ = 'V', then 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.
                     If JOBZ = 'N', then IFAIL is not referenced.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  if INFO = i, then i eigenvectors failed to converge.
                           Their indices are stored in array IFAIL.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date
           June 2016

       Further Details:

             All details about the 2stage techniques are available in:

             Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
             Parallel reduction to condensed forms for symmetric eigenvalue problems
             using aggregated fine-grained and memory-aware kernels. In Proceedings
             of 2011 International Conference for High Performance Computing,
             Networking, Storage and Analysis (SC '11), New York, NY, USA,
             Article 8 , 11 pages.
             http://doi.acm.org/10.1145/2063384.2063394

             A. Haidar, J. Kurzak, P. Luszczek, 2013.
             An improved parallel singular value algorithm and its implementation
             for multicore hardware, In Proceedings of 2013 International Conference
             for High Performance Computing, Networking, Storage and Analysis (SC '13).
             Denver, Colorado, USA, 2013.
             Article 90, 12 pages.
             http://doi.acm.org/10.1145/2503210.2503292

             A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
             A novel hybrid CPU-GPU generalized eigensolver for electronic structure
             calculations based on fine-grained memory aware tasks.
             International Journal of High Performance Computing Applications.
             Volume 28 Issue 2, Pages 196-209, May 2014.
             http://hpc.sagepub.com/content/28/2/196

   subroutine zhegv (integer ITYPE, character JOBZ, character UPLO, integer N, complex*16,
       dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB, double
       precision, dimension( * ) W, complex*16, dimension( * ) WORK, integer LWORK, double
       precision, dimension( * ) RWORK, integer INFO)
       ZHEGV

       Purpose:

            ZHEGV computes all the eigenvalues, and optionally, the eigenvectors
            of a complex generalized Hermitian-definite eigenproblem, of the form
            A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
            Here A and B are assumed to be Hermitian and B is also
            positive definite.

       Parameters
           ITYPE

                     ITYPE is INTEGER
                     Specifies the problem type to be solved:
                     = 1:  A*x = (lambda)*B*x
                     = 2:  A*B*x = (lambda)*x
                     = 3:  B*A*x = (lambda)*x

           JOBZ

                     JOBZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only;
                     = 'V':  Compute eigenvalues and eigenvectors.

           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  Upper triangles of A and B are stored;
                     = 'L':  Lower triangles of A and B are stored.

           N

                     N is INTEGER
                     The order of the matrices A and B.  N >= 0.

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, the Hermitian 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
                     matrix Z of eigenvectors.  The eigenvectors are normalized
                     as follows:
                     if ITYPE = 1 or 2, Z**H*B*Z = I;
                     if ITYPE = 3, Z**H*inv(B)*Z = I.
                     If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
                     or the lower triangle (if UPLO='L') of A, including the
                     diagonal, is destroyed.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,N).

           B

                     B is COMPLEX*16 array, dimension (LDB, N)
                     On entry, the Hermitian positive definite matrix B.
                     If UPLO = 'U', the leading N-by-N upper triangular part of B
                     contains the upper triangular part of the matrix B.
                     If UPLO = 'L', the leading N-by-N lower triangular part of B
                     contains the lower triangular part of the matrix B.

                     On exit, if INFO <= N, the part of B containing the matrix is
                     overwritten by the triangular factor U or L from the Cholesky
                     factorization B = U**H*U or B = L*L**H.

           LDB

                     LDB is INTEGER
                     The leading dimension of the array B.  LDB >= max(1,N).

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     If INFO = 0, the eigenvalues in ascending order.

           WORK

                     WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The length of the array WORK.  LWORK >= max(1,2*N-1).
                     For optimal efficiency, LWORK >= (NB+1)*N,
                     where NB is the blocksize for ZHETRD 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.

           RWORK

                     RWORK is DOUBLE PRECISION array, dimension (max(1, 3*N-2))

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  ZPOTRF or ZHEEV returned an error code:
                        <= N:  if INFO = i, ZHEEV failed to converge;
                               i off-diagonal elements of an intermediate
                               tridiagonal form did not converge to zero;
                        > 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.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date
           December 2016

   subroutine zhegv_2stage (integer ITYPE, character JOBZ, character UPLO, integer N, complex*16,
       dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB, double
       precision, dimension( * ) W, complex*16, dimension( * ) WORK, integer LWORK, double
       precision, dimension( * ) RWORK, integer INFO)
       ZHEGV_2STAGE

       Purpose:

            ZHEGV_2STAGE computes all the eigenvalues, and optionally, the eigenvectors
            of a complex generalized Hermitian-definite eigenproblem, of the form
            A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
            Here A and B are assumed to be Hermitian and B is also
            positive definite.
            This routine use the 2stage technique for the reduction to tridiagonal
            which showed higher performance on recent architecture and for large
            sizes N>2000.

       Parameters
           ITYPE

                     ITYPE is INTEGER
                     Specifies the problem type to be solved:
                     = 1:  A*x = (lambda)*B*x
                     = 2:  A*B*x = (lambda)*x
                     = 3:  B*A*x = (lambda)*x

           JOBZ

                     JOBZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only;
                     = 'V':  Compute eigenvalues and eigenvectors.
                             Not available in this release.

           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  Upper triangles of A and B are stored;
                     = 'L':  Lower triangles of A and B are stored.

           N

                     N is INTEGER
                     The order of the matrices A and B.  N >= 0.

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, the Hermitian 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
                     matrix Z of eigenvectors.  The eigenvectors are normalized
                     as follows:
                     if ITYPE = 1 or 2, Z**H*B*Z = I;
                     if ITYPE = 3, Z**H*inv(B)*Z = I.
                     If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
                     or the lower triangle (if UPLO='L') of A, including the
                     diagonal, is destroyed.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,N).

           B

                     B is COMPLEX*16 array, dimension (LDB, N)
                     On entry, the Hermitian positive definite matrix B.
                     If UPLO = 'U', the leading N-by-N upper triangular part of B
                     contains the upper triangular part of the matrix B.
                     If UPLO = 'L', the leading N-by-N lower triangular part of B
                     contains the lower triangular part of the matrix B.

                     On exit, if INFO <= N, the part of B containing the matrix is
                     overwritten by the triangular factor U or L from the Cholesky
                     factorization B = U**H*U or B = L*L**H.

           LDB

                     LDB is INTEGER
                     The leading dimension of the array B.  LDB >= max(1,N).

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     If INFO = 0, the eigenvalues in ascending order.

           WORK

                     WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The length of the array WORK. LWORK >= 1, when N <= 1;
                     otherwise
                     If JOBZ = 'N' and N > 1, LWORK must be queried.
                                              LWORK = MAX(1, dimension) where
                                              dimension = max(stage1,stage2) + (KD+1)*N + N
                                                        = N*KD + N*max(KD+1,FACTOPTNB)
                                                          + max(2*KD*KD, KD*NTHREADS)
                                                          + (KD+1)*N + N
                                              where KD is the blocking size of the reduction,
                                              FACTOPTNB is the blocking used by the QR or LQ
                                              algorithm, usually FACTOPTNB=128 is a good choice
                                              NTHREADS is the number of threads used when
                                              openMP compilation is enabled, otherwise =1.
                     If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available

                     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.

           RWORK

                     RWORK is DOUBLE PRECISION array, dimension (max(1, 3*N-2))

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  ZPOTRF or ZHEEV returned an error code:
                        <= N:  if INFO = i, ZHEEV failed to converge;
                               i off-diagonal elements of an intermediate
                               tridiagonal form did not converge to zero;
                        > 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.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date
           November 2017

       Further Details:

             All details about the 2stage techniques are available in:

             Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
             Parallel reduction to condensed forms for symmetric eigenvalue problems
             using aggregated fine-grained and memory-aware kernels. In Proceedings
             of 2011 International Conference for High Performance Computing,
             Networking, Storage and Analysis (SC '11), New York, NY, USA,
             Article 8 , 11 pages.
             http://doi.acm.org/10.1145/2063384.2063394

             A. Haidar, J. Kurzak, P. Luszczek, 2013.
             An improved parallel singular value algorithm and its implementation
             for multicore hardware, In Proceedings of 2013 International Conference
             for High Performance Computing, Networking, Storage and Analysis (SC '13).
             Denver, Colorado, USA, 2013.
             Article 90, 12 pages.
             http://doi.acm.org/10.1145/2503210.2503292

             A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
             A novel hybrid CPU-GPU generalized eigensolver for electronic structure
             calculations based on fine-grained memory aware tasks.
             International Journal of High Performance Computing Applications.
             Volume 28 Issue 2, Pages 196-209, May 2014.
             http://hpc.sagepub.com/content/28/2/196

   subroutine zhegvd (integer ITYPE, character JOBZ, character UPLO, integer N, complex*16,
       dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB, double
       precision, dimension( * ) W, complex*16, dimension( * ) WORK, integer LWORK, double
       precision, dimension( * ) RWORK, integer LRWORK, integer, dimension( * ) IWORK, integer
       LIWORK, integer INFO)
       ZHEGVD

       Purpose:

            ZHEGVD computes all the eigenvalues, and optionally, the eigenvectors
            of a complex generalized Hermitian-definite eigenproblem, of the form
            A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
            B are assumed to be Hermitian and B is also positive definite.
            If eigenvectors are desired, it uses a divide and conquer algorithm.

            The divide and conquer algorithm makes very mild assumptions about
            floating point arithmetic. It will work on machines with a guard
            digit in add/subtract, or on those binary machines without guard
            digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
            Cray-2. It could conceivably fail on hexadecimal or decimal machines
            without guard digits, but we know of none.

       Parameters
           ITYPE

                     ITYPE is INTEGER
                     Specifies the problem type to be solved:
                     = 1:  A*x = (lambda)*B*x
                     = 2:  A*B*x = (lambda)*x
                     = 3:  B*A*x = (lambda)*x

           JOBZ

                     JOBZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only;
                     = 'V':  Compute eigenvalues and eigenvectors.

           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  Upper triangles of A and B are stored;
                     = 'L':  Lower triangles of A and B are stored.

           N

                     N is INTEGER
                     The order of the matrices A and B.  N >= 0.

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, the Hermitian 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
                     matrix Z of eigenvectors.  The eigenvectors are normalized
                     as follows:
                     if ITYPE = 1 or 2, Z**H*B*Z = I;
                     if ITYPE = 3, Z**H*inv(B)*Z = I.
                     If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
                     or the lower triangle (if UPLO='L') of A, including the
                     diagonal, is destroyed.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,N).

           B

                     B is COMPLEX*16 array, dimension (LDB, N)
                     On entry, the Hermitian matrix B.  If UPLO = 'U', the
                     leading N-by-N upper triangular part of B contains the
                     upper triangular part of the matrix B.  If UPLO = 'L',
                     the leading N-by-N lower triangular part of B contains
                     the lower triangular part of the matrix B.

                     On exit, if INFO <= N, the part of B containing the matrix is
                     overwritten by the triangular factor U or L from the Cholesky
                     factorization B = U**H*U or B = L*L**H.

           LDB

                     LDB is INTEGER
                     The leading dimension of the array B.  LDB >= max(1,N).

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     If INFO = 0, the eigenvalues in ascending order.

           WORK

                     WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The length of the array WORK.
                     If N <= 1,                LWORK >= 1.
                     If JOBZ  = 'N' and N > 1, LWORK >= N + 1.
                     If JOBZ  = 'V' and N > 1, LWORK >= 2*N + N**2.

                     If LWORK = -1, then a workspace query is assumed; the routine
                     only calculates the optimal sizes of the WORK, RWORK and
                     IWORK arrays, returns these values as the first entries of
                     the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           RWORK

                     RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
                     On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.

           LRWORK

                     LRWORK is INTEGER
                     The dimension of the array RWORK.
                     If N <= 1,                LRWORK >= 1.
                     If JOBZ  = 'N' and N > 1, LRWORK >= N.
                     If JOBZ  = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2.

                     If LRWORK = -1, then a workspace query is assumed; the
                     routine only calculates the optimal sizes of the WORK, RWORK
                     and IWORK arrays, returns these values as the first entries
                     of the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           IWORK

                     IWORK is INTEGER array, dimension (MAX(1,LIWORK))
                     On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.

           LIWORK

                     LIWORK is INTEGER
                     The dimension of the array IWORK.
                     If N <= 1,                LIWORK >= 1.
                     If JOBZ  = 'N' and N > 1, LIWORK >= 1.
                     If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.

                     If LIWORK = -1, then a workspace query is assumed; the
                     routine only calculates the optimal sizes of the WORK, RWORK
                     and IWORK arrays, returns these values as the first entries
                     of the WORK, RWORK and IWORK arrays, and no error message
                     related to LWORK or LRWORK or LIWORK is issued by XERBLA.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  ZPOTRF or ZHEEVD returned an error code:
                        <= N:  if INFO = i and JOBZ = 'N', then the algorithm
                               failed to converge; i off-diagonal elements of an
                               intermediate tridiagonal form did not converge to
                               zero;
                               if INFO = i and JOBZ = 'V', then the algorithm
                               failed to compute an eigenvalue while working on
                               the submatrix lying in rows and columns INFO/(N+1)
                               through mod(INFO,N+1);
                        > 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.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date
           December 2016

       Further Details:

             Modified so that no backsubstitution is performed if ZHEEVD fails to
             converge (NEIG in old code could be greater than N causing out of
             bounds reference to A - reported by Ralf Meyer).  Also corrected the
             description of INFO and the test on ITYPE. Sven, 16 Feb 05.

       Contributors:
           Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

   subroutine zhegvx (integer ITYPE, character JOBZ, character RANGE, character UPLO, integer N,
       complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer
       LDB, double precision VL, double precision VU, integer IL, integer IU, double precision
       ABSTOL, integer M, double precision, dimension( * ) W, complex*16, dimension( ldz, * ) Z,
       integer LDZ, complex*16, dimension( * ) WORK, integer LWORK, double precision, dimension(
       * ) RWORK, integer, dimension( * ) IWORK, integer, dimension( * ) IFAIL, integer INFO)
       ZHEGVX

       Purpose:

            ZHEGVX computes selected eigenvalues, and optionally, eigenvectors
            of a complex generalized Hermitian-definite eigenproblem, of the form
            A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
            B are assumed to be Hermitian and B is also positive definite.
            Eigenvalues and eigenvectors can be selected by specifying either a
            range of values or a range of indices for the desired eigenvalues.

       Parameters
           ITYPE

                     ITYPE is INTEGER
                     Specifies the problem type to be solved:
                     = 1:  A*x = (lambda)*B*x
                     = 2:  A*B*x = (lambda)*x
                     = 3:  B*A*x = (lambda)*x

           JOBZ

                     JOBZ is CHARACTER*1
                     = 'N':  Compute eigenvalues only;
                     = 'V':  Compute eigenvalues and eigenvectors.

           RANGE

                     RANGE is CHARACTER*1
                     = 'A': all eigenvalues will be found.
                     = 'V': all eigenvalues in the half-open interval (VL,VU]
                            will be found.
                     = 'I': the IL-th through IU-th eigenvalues will be found.

           UPLO

                     UPLO is CHARACTER*1
                     = 'U':  Upper triangles of A and B are stored;
                     = 'L':  Lower triangles of A and B are stored.

           N

                     N is INTEGER
                     The order of the matrices A and B.  N >= 0.

           A

                     A is COMPLEX*16 array, dimension (LDA, N)
                     On entry, the Hermitian 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,  the lower triangle (if UPLO='L') or the upper
                     triangle (if UPLO='U') of A, including the diagonal, is
                     destroyed.

           LDA

                     LDA is INTEGER
                     The leading dimension of the array A.  LDA >= max(1,N).

           B

                     B is COMPLEX*16 array, dimension (LDB, N)
                     On entry, the Hermitian matrix B.  If UPLO = 'U', the
                     leading N-by-N upper triangular part of B contains the
                     upper triangular part of the matrix B.  If UPLO = 'L',
                     the leading N-by-N lower triangular part of B contains
                     the lower triangular part of the matrix B.

                     On exit, if INFO <= N, the part of B containing the matrix is
                     overwritten by the triangular factor U or L from the Cholesky
                     factorization B = U**H*U or B = L*L**H.

           LDB

                     LDB is INTEGER
                     The leading dimension of the array B.  LDB >= max(1,N).

           VL

                     VL is DOUBLE PRECISION

                     If RANGE='V', the lower bound of the interval to
                     be searched for eigenvalues. VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           VU

                     VU is DOUBLE PRECISION

                     If RANGE='V', the upper bound of the interval to
                     be searched for eigenvalues. VL < VU.
                     Not referenced if RANGE = 'A' or 'I'.

           IL

                     IL is INTEGER

                     If RANGE='I', the index of the
                     smallest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           IU

                     IU is INTEGER

                     If RANGE='I', the index of the
                     largest eigenvalue to be returned.
                     1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                     Not referenced if RANGE = 'A' or 'V'.

           ABSTOL

                     ABSTOL is DOUBLE PRECISION
                     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 + EPS *   max( |a|,|b| ) ,

                     where EPS is the machine precision.  If ABSTOL is less than
                     or equal to zero, then  EPS*|T|  will be used in its place,
                     where |T| is the 1-norm of the tridiagonal matrix obtained
                     by reducing C to tridiagonal form, where C is the symmetric
                     matrix of the standard symmetric problem to which the
                     generalized problem is transformed.

                     Eigenvalues will be computed most accurately when ABSTOL is
                     set to twice the underflow threshold 2*DLAMCH('S'), not zero.
                     If this routine returns with INFO>0, indicating that some
                     eigenvectors did not converge, try setting ABSTOL to
                     2*DLAMCH('S').

           M

                     M is INTEGER
                     The total number of eigenvalues found.  0 <= M <= N.
                     If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.

           W

                     W is DOUBLE PRECISION array, dimension (N)
                     The first M elements contain the selected
                     eigenvalues in ascending order.

           Z

                     Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
                     If JOBZ = 'N', then Z is not referenced.
                     If JOBZ = 'V', then if INFO = 0, the first M columns of Z
                     contain the orthonormal eigenvectors of the matrix A
                     corresponding to the selected eigenvalues, with the i-th
                     column of Z holding the eigenvector associated with W(i).
                     The eigenvectors are normalized as follows:
                     if ITYPE = 1 or 2, Z**T*B*Z = I;
                     if ITYPE = 3, Z**T*inv(B)*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.
                     Note: the user must ensure that at least max(1,M) columns are
                     supplied in the array Z; if RANGE = 'V', the exact value of M
                     is not known in advance and an upper bound must be used.

           LDZ

                     LDZ is INTEGER
                     The leading dimension of the array Z.  LDZ >= 1, and if
                     JOBZ = 'V', LDZ >= max(1,N).

           WORK

                     WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

           LWORK

                     LWORK is INTEGER
                     The length of the array WORK.  LWORK >= max(1,2*N).
                     For optimal efficiency, LWORK >= (NB+1)*N,
                     where NB is the blocksize for ZHETRD 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.

           RWORK

                     RWORK is DOUBLE PRECISION array, dimension (7*N)

           IWORK

                     IWORK is INTEGER array, dimension (5*N)

           IFAIL

                     IFAIL is INTEGER array, dimension (N)
                     If JOBZ = 'V', then 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.
                     If JOBZ = 'N', then IFAIL is not referenced.

           INFO

                     INFO is INTEGER
                     = 0:  successful exit
                     < 0:  if INFO = -i, the i-th argument had an illegal value
                     > 0:  ZPOTRF or ZHEEVX returned an error code:
                        <= N:  if INFO = i, ZHEEVX failed to converge;
                               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.

       Author
           Univ. of Tennessee

           Univ. of California Berkeley

           Univ. of Colorado Denver

           NAG Ltd.

       Date
           June 2016

       Contributors:
           Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

Author

       Generated automatically by Doxygen for LAPACK from the source code.