geneigbrobl.F90 6.18 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7
MODULE m_geneigprobl
8
  USE m_juDFT
9 10 11 12 13
  !**********************************************************
  !     Solve the generalized eigenvalue problem
  !     Frank Freimuth, November 2006
  !**********************************************************
CONTAINS
14
  SUBROUTINE geneigprobl(nbasfcn, nsize,neigd,l_J,eig,ne,a_r,b_r,z_r,a_c,b_c,z_c)
15
#include"cpp_double.h"
16
    USE m_packed_to_full
17 18 19 20 21 22 23 24 25 26 27 28
    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

29 30 31 32
    REAL,OPTIONAL, ALLOCATABLE, INTENT (INOUT) :: a_r(:),b_r(:)
    REAL,OPTIONAL, ALLOCATABLE, INTENT (INOUT) :: z_r(:,:)
    COMPLEX,OPTIONAL, ALLOCATABLE, INTENT (INOUT) :: a_c(:),b_c(:)
    COMPLEX,OPTIONAL, ALLOCATABLE, INTENT (INOUT) :: z_c(:,:)
33 34 35 36


    ! ... Local Variables ..

37
    INTEGER iind,ind1,ind2,info,lwork,liwork,lrwork,err,i
38 39 40 41 42
    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.  

43
    REAL toler, eigTemp(nsize)
44 45 46

    REAL,    ALLOCATABLE :: work(:)
    INTEGER, ALLOCATABLE :: iwork(:),isuppz(:)
47 48
    REAL, ALLOCATABLE :: largea_r(:,:),largeb_r(:,:)
    COMPLEX, ALLOCATABLE :: largea_c(:,:),largeb_c(:,:)
49
    COMPLEX,ALLOCATABLE :: cwork(:)
50 51 52 53 54

    LOGICAL :: l_real

    l_real=PRESENT(a_r)

55 56 57 58 59 60

    !**********************************
    !expand from packed to full storage: full storage lapack-routines
    !are faster than the packed lapack-routines.
    !**********************************
    !hamiltonian
61 62 63 64 65 66
    IF (l_real) THEN
       call packed_to_full(nsize,a_r,largea_r)
       DEALLOCATE (a_r)
       call packed_to_full(nsize,b_r,largeb_r)
       DEALLOCATE (b_r)
    ELSE
67 68 69 70
       call packed_to_full(nsize,a_c,largea_c)
       DEALLOCATE (a_c)
       call packed_to_full(nsize,b_c,largeb_c)
       DEALLOCATE (b_c)
71 72 73 74 75 76
    ENDIF



    IF (l_real) then
    CALL CPP_LAPACK_spotrf('U',nsize,largeb_r,nsize,info)
77 78
    IF (info /= 0)  CALL juDFT_error("error in spotrf",calledby ="geneigprobl")

79
    CALL CPP_LAPACK_ssygst(1,'U',nsize,largea_r,nsize,largeb_r,nsize,info)
80 81
    IF (info /= 0)  CALL juDFT_error("error in ssygst",calledby ="geneigprobl")

82
    toler = 2.0*TINY(toler)
83 84 85 86 87 88 89 90 91
    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")
92 93
    IF (ALLOCATED(z_r)) THEN
       IF (.NOT.(SIZE(z_r,1)==nbasfcn.AND.SIZE(z_r,2)==neigd)) DEALLOCATE(z_r)
94
    ENDIF
95 96
    IF (.NOT.ALLOCATED(z_r)) THEN
       ALLOCATE ( z_r(nbasfcn,neigd), stat=err )
97
       IF (err/=0) THEN
98
          WRITE(*,*) nbasfcn,neigd,err
99 100 101
          CALL juDFT_error("error allocating z",calledby ="geneigprobl")
       ENDIF
    ENDIF
102 103
    sizez= SIZE(z_r,1) 
    iu   = MIN(nsize,neigd)
104
    IF (l_J) THEN
105
       CALL CPP_LAPACK_ssyevr('N','I','U',nsize,largea_r, nsize,lb,ub,1,iu,toler,ne,eigTemp,z_r,&
106 107
            sizez,isuppz,work,lwork,iwork,liwork,info)
    ELSE
108
       CALL CPP_LAPACK_ssyevr('V','I','U',nsize,largea_r,nsize,lb,ub,1,iu,toler,ne,eigTemp,z_r,&
109 110 111 112 113
            sizez,isuppz,work,lwork,iwork,liwork,info)
    ENDIF
    IF (info /= 0)  CALL juDFT_error("error in ssyevr",calledby ="geneigprobl")
    DEALLOCATE (isuppz,work,iwork)

114 115 116
    CALL CPP_LAPACK_strtrs('U','N','N',nsize,ne,largeb_r, nsize,z_r,sizez,info)
    IF (info /= 0)  CALL juDFT_error("error in strtrs",calledby ="geneigprobl")
 ELSE
117

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

121
    CALL CPP_LAPACK_chegst(1,'U',nsize,largea_c,nsize,largeb_c,nsize,info)
122 123
    IF (info /= 0)  CALL juDFT_error(" error in chegst",calledby ="geneigprobl")

124
    toler = 2.0*TINY(toler)
125 126 127 128 129 130 131 132 133 134 135 136 137
    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")
138 139
    IF (ALLOCATED(z_c)) THEN
       IF (.NOT.(SIZE(z_c,1)==nbasfcn.AND.SIZE(z_c,2)==neigd)) DEALLOCATE(z_c)
140
    ENDIF
141 142
    IF (.NOT.ALLOCATED(z_c)) THEN
       ALLOCATE ( z_c(nbasfcn,neigd), stat=err )
143
       IF (err/=0) THEN
144
          WRITE(*,*) nbasfcn,neigd,err
145 146 147
          CALL juDFT_error("error allocating z",calledby ="geneigprobl")
       ENDIF
    ENDIF
148 149
    sizez= SIZE(z_c,1) 
    iu   = MIN(nsize,neigd)
150
    IF (l_J) THEN
151
       CALL CPP_LAPACK_cheevr('N','I','U',nsize,largea_c, nsize,lb,ub,1,iu,toler,ne,eigTemp,z_c,&
152 153
            sizez,isuppz,cwork,lwork,work,lrwork,iwork,liwork,info)
    ELSE
154
       CALL CPP_LAPACK_cheevr('V','I','U',nsize,largea_c, nsize,lb,ub,1,iu,toler,ne,eigTemp,z_c,&
155 156 157 158
            sizez,isuppz,cwork,lwork,work,lrwork,iwork,liwork,info)
    ENDIF
    IF (info /= 0)  CALL juDFT_error("error in cheevr",calledby ="geneigprobl")
    DEALLOCATE ( isuppz )
159 160 161
    DEALLOCATE ( work   )
    DEALLOCATE ( iwork  )
    DEALLOCATE ( cwork  )
162

163
    CALL CPP_LAPACK_ctrtrs('U','N','N',nsize,ne,largeb_c, nsize,z_c,sizez,info)
164
    IF (info /= 0)  CALL juDFT_error("error in ctrtrs",calledby ="geneigprobl")
165
 ENDIF
166

167 168 169
 DO i = 1, neigd
    eig(i) = eigTemp(i)
 END DO
170

171
END SUBROUTINE geneigprobl
172 173
END MODULE m_geneigprobl