aline.F90 7.25 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 8 9 10
MODULE m_aline
  USE m_juDFT
CONTAINS
  SUBROUTINE aline(eig_id, nk,atoms,DIMENSION,sym,&
11
       cell,input, jsp,el,usdus,lapw,tlmplm, noco, oneD,eig,ne,zMat,hmat,smat)
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 45 46
    !************************************************************************
    !*                                                                      *
    !*     eigensystem-solver for moderatly-well converged potentials       *
    !*     a*z=e*b*z is transformed to h*z'=e*s*z' , whereby                *
    !*     h=C^T*a*C, s=C^T*b*C and z'=C^(-1)*z, when C is z of the last    *
    !*     iteration (lapw%nv*ne-array)                                          *
    !*     For ne<<lapw%nv the matrixsize is significantly reduced               *
    !*     aline uses ESSL-calls (use LAPACK's reduc3, tred3, bisect,       *
    !*     tinvit, trback and rebk3  if no ESSL available):                 *
    !*     SSPMV:  matrix-vector multiplication for symmetric matrices      *
    !*             in packed storage.                                       *
    !*     SSYGV:  eigensystem-solver for symmetric, real h and positive    *
    !*             definite, real, symmetric s using Cholesky-factorisation *
    !*             tridiagonalisation and a QL-algorithm.                   *
    !*     For all eigenvalues are needed, DSYGV should perform better      *
    !*     then seclr4 (hope so)                                            *
    !*                                                     Gustav           * *                                                                      *
    !************************************************************************

#include"cpp_double.h"
    USE m_abcof
    USE m_hssrwu
    USE m_eig66_io
    USE m_types
    IMPLICIT NONE
    TYPE(t_dimension),INTENT(IN)   :: DIMENSION
    TYPE(t_oneD),INTENT(IN)        :: oneD
    TYPE(t_input),INTENT(IN)       :: input
    TYPE(t_noco),INTENT(IN)        :: noco
    TYPE(t_sym),INTENT(IN)         :: sym
    TYPE(t_cell),INTENT(IN)        :: cell
    TYPE(t_atoms),INTENT(IN)       :: atoms
    TYPE(t_usdus),INTENT(IN)       :: usdus
    TYPE(t_tlmplm),INTENT(IN)      :: tlmplm
    TYPE(t_lapw),INTENT(INOUT)     :: lapw
47
    TYPE(t_mat),INTENT(INOUT)      :: zMat
48

49 50 51 52 53 54 55
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: eig_id
    INTEGER, INTENT (IN) :: jsp,nk
    INTEGER, INTENT (OUT):: ne
    !     ..
    !     .. Array Arguments ..
56
    REAL,    INTENT (IN)  :: el(0:atoms%lmaxd,atoms%ntype,input%jspins)
57
    REAL,    INTENT (OUT) :: eig(DIMENSION%neigd)
58
    TYPE(t_mat),INTENT(IN):: hmat,smat
59

60 61 62 63 64 65 66
    !     ..
    !     .. Local Scalars ..
    INTEGER lhelp
    INTEGER i,info,j 
    !     ..
    !     .. Local Arrays ..
    COMPLEX, ALLOCATABLE :: acof(:,:,:),bcof(:,:,:),ccof(:,:,:,:)
67 68

    REAL,    ALLOCATABLE :: help_r(:),h_r(:,:),s_r(:,:) 
69 70
    REAL     CPP_BLAS_sdot
    EXTERNAL CPP_BLAS_sdot,CPP_BLAS_sspmv
71 72 73

    COMPLEX,   PARAMETER :: one_c=(1.0,0.0), zro_c=(0.0,0.0)
    COMPLEX, ALLOCATABLE :: help_c(:),h_c(:,:),s_c(:,:) 
74 75 76 77
    COMPLEX  CPP_BLAS_cdotc
    EXTERNAL CPP_BLAS_cdotc,CPP_BLAS_chpmv
    REAL,    ALLOCATABLE :: rwork(:)

78
    LOGICAL:: l_real
79
    l_real=zMat%l_real
80 81


82
    lhelp= MAX(lapw%nmat,(DIMENSION%neigd+2)*DIMENSION%neigd)
83
    CALL read_eig(eig_id,nk,jsp,neig=ne, eig=eig,zmat=zmat)
84 85 86 87
    IF (l_real) THEN
       ALLOCATE ( h_r(DIMENSION%neigd,DIMENSION%neigd),s_r(DIMENSION%neigd,DIMENSION%neigd) )
       h_r = 0.0 ; s_r=0.0
       ALLOCATE ( help_r(lhelp) )
88
    ELSE
89 90
       !     in outeig z is complex conjugated to make it usable for abcof. Here we 
       !                       first have to undo this  complex conjugation for the 
91 92
       ! multiplication with a and b matrices.

93
       zmat%data_c=conjg(zmat%data_c)
94 95 96 97
       ALLOCATE ( h_c(DIMENSION%neigd,DIMENSION%neigd),s_c(DIMENSION%neigd,DIMENSION%neigd) )
       h_c = 0.0 ; s_c=0.0
       ALLOCATE ( help_r(lhelp) )
    ENDIF
98 99
    !
    DO i = 1,ne
100
       IF (l_real) THEN
101
          help_r=MATMUL(hmat%data_r,zmat%data_r(:,i))
102
       ELSE
103
          help_c=MATMUL(hmat%data_c,zmat%data_c(:,i))
104
       ENDIF
105
       DO j = i,ne
106
          IF (l_real) THEN
107
             h_r(j,i)=dot_PRODUCT(zmat%data_r(:,j),help_r)
108
          ELSE
109
             h_c(j,i)=dot_PRODUCT(zmat%data_c(:,j),help_c)
110
          ENDIF
111 112 113 114
       END DO
    END DO

    DO i = 1,ne
115
       IF (l_real) THEN
116
          help_r=MATMUL(smat%data_r,zmat%data_r(:,i))
117
       ELSE
118
          help_c=MATMUL(smat%data_c,zmat%data_c(:,i))
119
       ENDIF
120
       DO j = i,ne
121
          IF (l_real) THEN
122
             s_r(j,i) = dot_product(zmat%data_r(:,j),help_r)
123
          ELSE
124
             s_c(j,i) =dot_PRODUCT(zmat%data_c(:,j),help_c)
125
          ENDIF
126 127 128
       END DO
    END DO

129 130
    ALLOCATE ( acof(DIMENSION%neigd,0:DIMENSION%lmd,atoms%nat),bcof(DIMENSION%neigd,0:DIMENSION%lmd,atoms%nat) )
    ALLOCATE ( ccof(-atoms%llod:atoms%llod,DIMENSION%neigd,atoms%nlod,atoms%nat) ) 
131

132
    !     conjugate again for use with abcof; finally use cdotc to revert again
133
    IF (.NOT.l_real) zMat%data_c = CONJG(zMat%data_c)
134 135
    if (noco%l_soc)  CALL juDFT_error("no SOC & reduced diagonalization",calledby="aline")

136 137
    CALL abcof(input,atoms,sym,cell,lapw,ne,&
         usdus,noco,1,oneD,acof,bcof,ccof,zMat)  ! ispin = 1&
138 139


140 141
    !
    CALL timestart("aline: hssr_wu")
142
    IF (l_real) THEN
143
       CALL hssr_wu(atoms,DIMENSION,sym, jsp,el,ne,usdus,lapw,input,&
144 145
            tlmplm, acof,bcof,ccof, h_r,s_r)
    ELSE
146
       CALL hssr_wu(atoms,DIMENSION,sym, jsp,el,ne,usdus,lapw,input,&
147 148
            tlmplm, acof,bcof,ccof, h_c=h_c,s_c=s_c)
    ENDIF
149 150 151 152 153 154

    DEALLOCATE ( ccof, bcof, acof )
    CALL timestop("aline: hssr_wu")
    CALL timestart("aline: seclr4")

    !
155 156 157 158 159 160 161 162
    IF (l_real) THEN
       !---> LAPACK call
       CALL CPP_LAPACK_ssygv(1,'V','L',ne,h_r,DIMENSION%neigd,s_r,DIMENSION%neigd,eig,help_r,lhelp,info)
    ELSE
       ALLOCATE ( rwork(MAX(1,3*ne-2)) )
       CALL CPP_LAPACK_chegv(1,'V','L',ne,h_c,DIMENSION%neigd,s_c,DIMENSION%neigd,eig,help_c,lhelp,rwork,info)
       DEALLOCATE ( rwork )
    ENDIF
163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
    IF (info /= 0) THEN
       WRITE (6,FMT=8000) info
       IF (i < 0) THEN
          WRITE(6,'(a7,i3,a22)') 'element',info,' has an illegal value'
       ELSEIF (i > ne) THEN
          WRITE(6,'(a2,i3,a22)') 's:',info-ne,' not positive definite'
       ELSE
          WRITE(6,'(a8,i3,a15)') 'argument',info,' not  converged'
       ENDIF
       CALL juDFT_error("Diagonalisation failed",calledby ='aline')
    ENDIF
8000 FORMAT (' AFTER CPP_LAPACK_ssygv: info=',i4)
    CALL timestop("aline: seclr4")

    DO i = 1,lapw%nmat
178
       IF (l_real) THEN
179
          help_r(:ne)=zMat%data_r(i,:ne)
180
       ELSE
181
          help_c(:ne)=zMat%data_c(i,:ne)
182
       END IF
183
       DO j = 1,ne
184 185
          IF (l_real) THEN
             !--->       for LAPACK call
186
             zMat%data_r(i,j) = CPP_BLAS_sdot(ne,help_r,1,h_r(1,j),1)
187
          ELSE
188
             zMat%data_c(i,j) = CPP_BLAS_cdotc(ne,help_c,1,h_c(1,j),1)
189
          ENDIF
190 191 192 193 194
       END DO
    END DO

  END SUBROUTINE aline
END MODULE m_aline