aline.F90 7.79 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, bkpt,eig,ne,zMat,a_r,b_r,a_c,b_c)
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 48
    TYPE(t_zMat),INTENT(INOUT)     :: zMat

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 ..
Daniel Wortmann's avatar
Daniel Wortmann committed
56
    REAL,    INTENT (IN)  :: el(0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd)
57
    REAL,    INTENT (OUT) :: eig(DIMENSION%neigd),bkpt(3)
58 59 60
    REAL,OPTIONAL,    INTENT (IN)  :: a_r(:),b_r(:)!(matsize)
    COMPLEX,OPTIONAL, INTENT (IN)  :: a_c(:),b_c(:)!(matsize)

61 62 63 64 65 66 67 68
    !     ..
    !     .. Local Scalars ..
    INTEGER lhelp
    INTEGER i,info,j 
    !     ..
    !     .. Local Arrays ..
    INTEGER kveclo(atoms%nlotot)
    COMPLEX, ALLOCATABLE :: acof(:,:,:),bcof(:,:,:),ccof(:,:,:,:)
69 70

    REAL,    ALLOCATABLE :: help_r(:),h_r(:,:),s_r(:,:) 
71 72
    REAL     CPP_BLAS_sdot
    EXTERNAL CPP_BLAS_sdot,CPP_BLAS_sspmv
73 74 75

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

80
    LOGICAL:: l_real
81
    l_real=zMat%l_real
82 83


84
    lhelp= MAX(lapw%nmat,(DIMENSION%neigd+2)*DIMENSION%neigd)
85 86 87
#ifdef __PGI
        STOP "CODE NOT WORKING WITH PGI"
#else
88
    IF (l_real) THEN
89
       CALL read_eig(eig_id,nk,jsp,bk=bkpt,neig=ne,nv=lapw%nv(jsp),nmat=lapw%nmat, eig=eig,kveclo=kveclo,z=zMat%z_r)
90 91 92
       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) )
93
    ELSE
94
       CALL read_eig(eig_id,nk,jsp,bk=bkpt,neig=ne,nv=lapw%nv(jsp),nmat=lapw%nmat, eig=eig,kveclo=kveclo,z=zMat%z_c)
95 96
       !     in outeig z is complex conjugated to make it usable for abcof. Here we 
       !                       first have to undo this  complex conjugation for the 
97 98 99
       ! multiplication with a and b matrices.

       zmat%z_c=conjg(zmat%z_c)
100 101 102 103
       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
104
#endif
105 106
    !
    DO i = 1,ne
107
       IF (l_real) THEN
108
          CALL CPP_BLAS_sspmv('U',lapw%nmat,1.0,a_r,zMat%z_r(1,i),1,0.0,help_r,1)
109
       ELSE
110
          CALL CPP_BLAS_chpmv('U',lapw%nmat,one_c,a_c,zMat%z_c(1,i),1,zro_c,help_c,1)
111
       ENDIF
112
       DO j = i,ne
113
          IF (l_real) THEN
114
             h_r(j,i) = CPP_BLAS_sdot(lapw%nmat,zMat%z_r(1,j),1,help_r,1)
115
          ELSE
116
             h_c(j,i) = CPP_BLAS_cdotc(lapw%nmat,zMat%z_c(1,j),1,help_c,1)
117
          ENDIF
118 119 120 121
       END DO
    END DO

    DO i = 1,ne
122
       IF (l_real) THEN
123
          CALL CPP_BLAS_sspmv('U',lapw%nmat,1.0,b_r,zMat%z_r(1,i),1,0.0,help_r,1)
124
       ELSE
125
          CALL CPP_BLAS_chpmv('U',lapw%nmat,one_c,b_c,zMat%z_c(1,i),1,zro_c,help_c,1)
126
       ENDIF
127
       DO j = i,ne
128
          IF (l_real) THEN
129
             s_r(j,i) = CPP_BLAS_sdot(lapw%nmat,zMat%z_r(1,j),1,help_r,1)
130
          ELSE
131
             s_c(j,i) = CPP_BLAS_cdotc(lapw%nmat,zMat%z_c(1,j),1,help_c,1)
132
          ENDIF
133 134 135
       END DO
    END DO

Daniel Wortmann's avatar
Daniel Wortmann committed
136 137
    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) ) 
138

139
    !     conjugate again for use with abcof; finally use cdotc to revert again
140
    IF (.NOT.l_real) zMat%z_c = CONJG(zMat%z_c)
141 142 143
    if (noco%l_soc)  CALL juDFT_error("no SOC & reduced diagonalization",calledby="aline")

    CALL abcof(input,atoms,DIMENSION%neigd,sym,cell, bkpt,lapw,ne,&
144
         usdus,noco,1,kveclo,oneD,acof,bcof,ccof,zMat)  ! ispin = 1&
145 146


147 148
    !
    CALL timestart("aline: hssr_wu")
149 150 151 152 153 154 155
    IF (l_real) THEN
       CALL hssr_wu(atoms,DIMENSION,sym, jsp,el,ne,usdus,lapw,&
            tlmplm, acof,bcof,ccof, h_r,s_r)
    ELSE
       CALL hssr_wu(atoms,DIMENSION,sym, jsp,el,ne,usdus,lapw,&
            tlmplm, acof,bcof,ccof, h_c=h_c,s_c=s_c)
    ENDIF
156 157 158 159 160 161

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

    !
162 163 164 165 166 167 168 169
    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
170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
    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
185
       IF (l_real) THEN
186
          help_r(:ne)=zMat%z_r(i,:ne)
187
       ELSE
188
          help_c(:ne)=zMat%z_c(i,:ne)
189
       END IF
190
       DO j = 1,ne
191 192
          IF (l_real) THEN
             !--->       for LAPACK call
193
             zMat%z_r(i,j) = CPP_BLAS_sdot(ne,help_r,1,h_r(1,j),1)
194
          ELSE
195
             zMat%z_c(i,j) = CPP_BLAS_cdotc(ne,help_c,1,h_c(1,j),1)
196
          ENDIF
197 198 199 200 201
       END DO
    END DO

  END SUBROUTINE aline
END MODULE m_aline