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