geneigbrobl.F90 6.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
MODULE m_geneigprobl
  use m_juDFT
  !**********************************************************
  !     Solve the generalized eigenvalue problem
  !     Frank Freimuth, November 2006
  !**********************************************************
CONTAINS
  SUBROUTINE geneigprobl(nbasfcn, nsize,neigd,l_J,a,b,z,eig,ne)
#include"cpp_double.h"
    IMPLICIT NONE

    ! ... Arguments ...

    INTEGER, INTENT (IN) :: nbasfcn
    INTEGER, INTENT (IN) :: neigd
    INTEGER, INTENT (IN) :: nsize
    LOGICAL, INTENT (IN) :: l_J

    REAL,    INTENT(OUT) :: eig(:)
    INTEGER, INTENT(OUT) :: ne
#ifdef CPP_F90

#ifdef CPP_INVERSION
    REAL,  INTENT (INOUT) :: a(:),b(:)
    REAL,  INTENT (INOUT) :: z(:,:)
#else
    COMPLEX, INTENT (INOUT)::a(:),b(:)
    COMPLEX, INTENT (INOUT) :: z(:,:)
#endif

#else

#ifdef CPP_INVERSION
    REAL, ALLOCATABLE, INTENT (INOUT) :: a(:),b(:)
    REAL, ALLOCATABLE, INTENT (INOUT) :: z(:,:)
#else
    COMPLEX, ALLOCATABLE, INTENT (INOUT) :: a(:),b(:)
    COMPLEX, ALLOCATABLE, INTENT (INOUT) :: z(:,:)
#endif

#endif

    ! ... Local Variables ..

45
    INTEGER iind,ind1,ind2,info,lwork,liwork,lrwork,err,i
46 47 48 49 50
    INTEGER sizez,iu 
    REAL :: lb,ub
    ! 'sizez' is needed, as some compilers sometimes produce errors,
    ! if the size command is used directly as a lapack argument.  

51
    REAL toler, eigTemp(nsize)
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130

    REAL,    ALLOCATABLE :: work(:)
    INTEGER, ALLOCATABLE :: iwork(:),isuppz(:)
#ifdef CPP_INVERSION
    REAL, ALLOCATABLE :: largea(:,:),largeb(:,:)
#else
    COMPLEX, ALLOCATABLE :: largea(:,:),largeb(:,:)
    COMPLEX,ALLOCATABLE :: cwork(:)
#endif

    !**********************************
    !expand from packed to full storage: full storage lapack-routines
    !are faster than the packed lapack-routines.
    !**********************************
    !hamiltonian
    ALLOCATE ( largea(nsize,nsize), stat=err )
    IF (err/=0)  CALL juDFT_error("error allocating largea",calledby&
         &     ="geneigprobl")
    iind = 0
    DO ind1 = 1, nsize
       DO ind2 = 1, ind1
          iind = iind+1
          largea(ind2,ind1) = a(iind)
       ENDDO
    ENDDO
    !save some storage by deallocation of unused array
#ifndef CPP_F90
    DEALLOCATE (a)
#endif
    !metric
    ALLOCATE ( largeb(nsize,nsize), stat=err )
    IF (err/=0)  CALL juDFT_error("error allocating largeb",calledby ="geneigprobl")
    iind=0
    DO ind1 = 1, nsize
       DO ind2 = 1, ind1
          iind = iind+1
          largeb(ind2,ind1) = b(iind)
       ENDDO
    ENDDO
    !save some storage by deallocation of unused array
#ifndef CPP_F90
    DEALLOCATE (b)
#endif



#ifdef CPP_INVERSION
    CALL CPP_LAPACK_spotrf('U',nsize,largeb,nsize,info)
    IF (info /= 0)  CALL juDFT_error("error in spotrf",calledby ="geneigprobl")

    CALL CPP_LAPACK_ssygst(1,'U',nsize,largea,nsize,largeb,nsize,info)
    IF (info /= 0)  CALL juDFT_error("error in ssygst",calledby ="geneigprobl")

    toler = 2.0*tiny(toler)
    liwork = 10*nsize
    ALLOCATE ( iwork(liwork), stat=err )
    IF (err/=0)  CALL juDFT_error("error allocating iwork",calledby ="geneigprobl")

    lwork = 26*nsize
    ALLOCATE ( work(lwork), stat=err )
    IF (err/=0)  CALL juDFT_error(" error allocating work",calledby ="geneigprobl")
    ALLOCATE ( isuppz(2*nsize), stat=err )
    IF (err /= 0)  CALL juDFT_error("error allocating isuppz",calledby ="geneigprobl")
#ifndef CPP_F90
    IF (allocated(z)) THEN
       IF (.not.(size(z,1)==nbasfcn.and.size(z,2)==neigd)) deallocate(z)
    ENDIF
    IF (.not.allocated(z)) THEN
       ALLOCATE ( z(nbasfcn,neigd), stat=err )
       IF (err/=0) THEN
          write(*,*) nbasfcn,neigd,err
          CALL juDFT_error("error allocating z",calledby ="geneigprobl")
       ENDIF
    ENDIF
#endif
    sizez= size(z,1) 
    iu   = min(nsize,neigd)
#ifndef CPP_F90
    IF (l_J) THEN
131
       CALL CPP_LAPACK_ssyevr('N','I','U',nsize,largea, nsize,lb,ub,1,iu,toler,ne,eigTemp,z,&
132 133
            sizez,isuppz,work,lwork,iwork,liwork,info)
    ELSE
134
       CALL CPP_LAPACK_ssyevr('V','I','U',nsize,largea,nsize,lb,ub,1,iu,toler,ne,eigTemp,z,&
135 136 137 138
            sizez,isuppz,work,lwork,iwork,liwork,info)
    ENDIF
#else
    eig = 0.0
139
    eigTemp = 0.0
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
#endif
    IF (info /= 0)  CALL juDFT_error("error in ssyevr",calledby ="geneigprobl")
    DEALLOCATE (isuppz,work,iwork)

    CALL CPP_LAPACK_strtrs('U','N','N',nsize,ne,largeb, nsize,z,sizez,info)
    IF (info /= 0)  CALL juDFT_error("error in strtrs",calledby&
         &     ="geneigprobl")
#else

    CALL CPP_LAPACK_cpotrf('U',nsize,largeb,nsize,info)
    IF (info /= 0)  CALL juDFT_error("error in cpotrf",calledby ="geneigprobl")

    CALL CPP_LAPACK_chegst(1,'U',nsize,largea,nsize,largeb,nsize,info)
    IF (info /= 0)  CALL juDFT_error(" error in chegst",calledby ="geneigprobl")

    toler = 2.0*tiny(toler)
    liwork = 50*nsize
    ALLOCATE ( iwork(liwork), stat=err )
    IF (err/=0)  CALL juDFT_error("error allocating iwork",calledby ="geneigprobl")

    lwork = 20*nsize
    ALLOCATE( cwork(lwork), stat=err )
    IF (err/=0)  CALL juDFT_error("error allocating cwork",calledby ="geneigprobl")
    ALLOCATE( isuppz(10*nsize), stat=err )
    IF (err/=0)  CALL juDFT_error("error allocating isuppz",calledby ="geneigprobl")

    lrwork = 84*nsize
    ALLOCATE (work(lrwork), stat=err )
    IF (err/=0)  CALL juDFT_error(" error allocating work",calledby ="geneigprobl")
#ifndef CPP_F90
    IF (allocated(z)) THEN
       IF (.not.(size(z,1)==nbasfcn.and.size(z,2)==neigd)) deallocate(z)
    ENDIF
    IF (.not.allocated(z)) THEN
       ALLOCATE ( z(nbasfcn,neigd), stat=err )
       IF (err/=0) THEN
          write(*,*) nbasfcn,neigd,err
          CALL juDFT_error("error allocating z",calledby ="geneigprobl")
       ENDIF
    ENDIF
#endif
    sizez= size(z,1) 
    iu   = min(nsize,neigd)
#ifndef CPP_F90
    IF (l_J) THEN
185
       CALL CPP_LAPACK_cheevr('N','I','U',nsize,largea, nsize,lb,ub,1,iu,toler,ne,eigTemp,z,&
186 187 188
            sizez,isuppz,cwork,lwork,work,lrwork,iwork,liwork,info)
    ELSE
#if (1==1)
189
       CALL CPP_LAPACK_cheevr('V','I','U',nsize,largea, nsize,lb,ub,1,iu,toler,ne,eigTemp,z,&
190 191 192
            sizez,isuppz,cwork,lwork,work,lrwork,iwork,liwork,info)
#else

193
       CALL CPP_LAPACK_cheevx('V','I','U',nsize,largea, nsize,lb,ub,1,iu,toler,ne,eigTemp,z,&
194 195 196 197 198
            sizez,cwork,lwork,work,iwork,isuppz,info)
#endif
    ENDIF
#else
    eig = 0.0
199
    eigTemp = 0.0
200 201 202 203 204 205 206 207 208 209 210 211
#endif
    IF (info /= 0)  CALL juDFT_error("error in cheevr",calledby ="geneigprobl")
    DEALLOCATE ( isuppz )
    deallocate ( work   )
    deallocate ( iwork  )
    deallocate ( cwork  )

    CALL CPP_LAPACK_ctrtrs('U','N','N',nsize,ne,largeb, nsize,z,sizez,info)
    IF (info /= 0)  CALL juDFT_error("error in ctrtrs",calledby ="geneigprobl")
#endif
    DEALLOCATE ( largea,largeb )

212 213 214 215
    DO i = 1, neigd
       eig(i) = eigTemp(i)
    END DO

216 217 218
  END SUBROUTINE geneigprobl
END MODULE m_geneigprobl