c============================================================================
c     Eigenissue solves the semi-definite, Hermitian, generalized 
c     eigenvalue problem: 
c
c         M.x=lambda.E.x    (where M,E Hermitian, E positive semidefinite)
C
C     For now, E is also diagonal
c
c     eigenissue(N,M,LDM,E,VALS,JOB,LOW,INFO,INFO2)
c
c     arguments:
c         N(input) INTEGER
c             size of matrix M
c         M(input/output) COMPLEX*16
c             on input dimension (LDM,LDM), M is a complex Hermitian matrix
c             on output it contains eigenvectors calculated by eigenmess
c         LDM(input) INTEGER
c             the leading order of M
c         E(input) DOUBLE PRECISION REAL
c             a vector length N containing the diagonal elements of E
c         VALS(output) DOUBLE PRECISION REAL
c             vector of length N, on output first LOW indices contain
c             the eigenvalues calculated by eigenmess
c         JOB(input) INTEGER
c             = 0:  Compute eigenvalues only;
c             = 1:  Compute eigenvalues and eigenvectors without constraints
c             = 2:  Compute eigenvalues and eigenvectors with constraints
c         LOW(input/output) INTEGER
c             on input number of lowest eigenvalues/vectors to be calculated
c             on output number of eigenvalues actually calculated
c         INFO(output) INTEGER
c             =0: then block matrix C was negative definite, invertible
c             <0: then info=-i, the ith argument to zposv was illegal
c             >0: if INFO = i, D(i,i) is exactly zero.  The factorization
c                has been completed, but the block diagonal matrix D(used
c                in the factorization of C for zhesv) is exactly singular,
c                so the solution could not be computed.
c         INFO2(output) INTEGER
c             =0: all eigenvectors converged, successful exit
c             <0: then info=-i, the ith argument to zheevx was illegal
c             >0: then info=i, then i eigenvectors failed to converge
c============================================================================

      subroutine eigenissue(N,M,LDM,E,VALS,JOB,LOW,INFO,INFO2)

* declare variables

      implicit none
      integer LDM,NUM,INFO,INFO2,N,JOB
      double complex M(LDM,LDM)
      double precision E(N),VALS(N)
      integer i,LOW


* test job input

      if(JOB.lt.0.OR.JOB.gt.2)then
         print*,'ERROR!! Invalid job: must be 0, 1, or 2.'
         return
      end if

* count number of zeros in E

      NUM=0
      do i=1,N
         if(E(i).eq.0)NUM=NUM+1
      end do

* calculate eigenvalues and eigenvectors

      call eigenmess(N,M,LDM,E,NUM,VALS,JOB,LOW,INFO,INFO2)
      return
      end

c==========================================================================
c     Eigenmess solves the semi-definite, Hermitian, generalized eigenvalue
c     problem, but only for the case where E is diagonal.
c
c     eigenmess(N,M,LDM,E,NUM,VALS,JOB,LOW,INFO,INFO2)
c
c     arguments:
c         N(input) INTEGER
c             size of matrix M
c         M(input/output) COMPLEX*16
c             on input dimension (LDM,LDM), M is a complex Hermitian matrix
c             on output it contains eigenvectors calculated by LAPACK
c             function zheevx
c         LDM(input) INTEGER
c             leading order of M
c         E(input) DOUBLE PRECISION REAL
c             a vector length N containing zeros and reals
c         NUM(input) INTEGER
c             number of zeros in E
c         VALS(output) DOUBLE PRECISION REAL
c             vector of length N, on output first LOW indices contain
c             the eigenvalues calculated by LAPACK function zheevx
c         JOB(input) INTEGER
c             = 0:  Compute eigenvalues only;
c             = 1:  Compute eigenvalues and eigenvectors without constraints
c             = 2:  Compute eigenvalues and eigenvectors with constraints
c         LOW(input) INTEGER
c             number of lowest eigenvalues/vectors to be calculated
c         INFO(output) INTEGER
c             =0: then block matrix C was negative definite, invertible
c             <0: then info=-i, the ith argument to zposv was illegal
c             >0: if INFO = i, D(i,i) is exactly zero.  The factorization
c                has been completed, but the block diagonal matrix D(used
c                in the factorization of C for zhesv) is exactly singular,
c                so the solution could not be computed.
c         INFO2(output) INTEGER
c             =0: all eigenvectors converged, successful exit
c             <0: then info=-i, the ith argument to zheevx was illegal
c             >0: then info=i, then i eigenvectors failed to converge
c============================================================================

      subroutine eigenmess(N,M,LDM,E,NUM,VALS,JOB,LOW,INFO,INFO2)

* declare variables

      implicit none
      integer NUM,N,LDM,LOW,IWORK(5*(N-NUM)),IFAIL(N-NUM),FOUND,
     *     PIV(NUM),JOB
      double complex M(LDM,LDM),PARTA(N-NUM,N-NUM),
     *     PARTB(N-NUM,NUM),PARTBdag(NUM,N-NUM),tmp(N,N),
     *     PARTC(NUM,NUM),WORK(33*(N-NUM)),WORK2(32*NUM),ALPHA,BETA
      double precision E(N),VALS(N),RWORK(7*(N-NUM))
      integer ZEROS(NUM),NONZEROS(N-NUM),INFO,i,j,zcount,nzcount,INFO2
      character jobz

* find zeros and nonzeros of E, store as two arrays

      zcount=1
      nzcount=1
      do i=1,N
         if(E(i).eq.0)then
            ZEROS(zcount)=i
            zcount=zcount+1
         else
            NONZEROS(nzcount)=i
            nzcount=nzcount+1
         end if
      end do

* make block matrices from large matrix M

      do i=1,N-NUM
         do j=i,N-NUM
            PARTA(i,j)=M(NONZEROS(i),NONZEROS(j))
            PARTA(j,i)=conjg(PARTA(i,j))
         end do
      end do
      do i=1,N-NUM
         do j=1,NUM
            if(NONZEROS(i).lt.ZEROS(j))then
               PARTB(i,j)=M(NONZEROS(i),ZEROS(j))
            else
               PARTB(i,j)=conjg(M(ZEROS(j),NONZEROS(i)))
            end if
         end do
      end do
      do i=1,NUM
         do j=i,NUM
            PARTC(i,j)=M(ZEROS(i),ZEROS(j))
         end do
      end do
      do i=1,NUM
         do j=1,N-NUM
            PARTBdag(i,j)=conjg(PARTB(j,i))
         end do
      end do

* calculate (A - B*(Cinverse)*Bdagger), stored in A to save memory

      ALPHA=(1D0,0)
      BETA=(0,0)
      call zhesv('U',NUM,N-NUM,PARTC,NUM,PIV,PARTBdag,NUM,
     *     WORK2,32*NUM,INFO)
      call zgemm('N','N',N-NUM,N-NUM,NUM,ALPHA,PARTB,N-NUM,
     *     PARTBdag,NUM,BETA,tmp,N)
      do i=1,N-NUM
         do j=1,N-NUM
            PARTA(i,j)=PARTA(i,j)-tmp(i,j)
         end do
      end do

* correct new A for rescaling

      do i=1,N-NUM
         do j=1,N-NUM
            PARTA(i,j)=PARTA(i,j)/
     *           sqrt(E(NONZEROS(i))*E(NONZEROS(j)))
         end do
      end do

* set job flag for zheevx (calculate eigenvectors or not)

      jobz='N'
      if(JOB.gt.0)jobz='V'

* solve eigensystem of new normalized A

      if(LOW.gt.N-NUM)LOW=N-NUM   !prevent from asking for more than possible
      call zheevx(jobz,'I','U',N-NUM,PARTA,N-NUM,0,0,1,LOW,
     *     0,FOUND,VALS,M,LDM,WORK,33*(N-NUM),RWORK,IWORK,IFAIL,INFO2)
      LOW=FOUND

** calculate eigenvectors with constraints if desired

      if(JOB.eq.2)then
* store X in bottom of tmp, so tmp column contains Y then X
         do i=1,N-NUM
            do j=1,LOW
               tmp(NUM+i,j)=M(i,j)
            end do
         end do
* resize eigenvectors X
         do i=1,N-NUM
            do j=1,LOW
               M(i,j)=M(i,j)/sqrt(E(NONZEROS(i)))
            end do
         end do
* calculate eigenvectors Y in tmp
         call zgemm('N','N',NUM,LOW,N-NUM,ALPHA,PARTBdag,
     *        NUM,M,LDM,BETA,tmp,N)
* sort X and Y back into M in the correct order
         do j=1,LOW
            do i=1,NUM
               M(ZEROS(i),j)=-tmp(i,j)
            end do
            do i=1,N-NUM
               M(NONZEROS(i),j)=tmp(NUM+i,j)
            end do
         end do
      endif

      return
      end

!Written by: Justin Jay Lambright between 21.6.2004 and 28.6.2004
!JJL added:  constraint eigenvectors calculation on 26.7.2004