eigen_diag.F90 9.91 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
MODULE m_eigen_diag
  USE m_juDFT
! the parameters are set to negative values to indicate that a particular solver is not compiled
10 11 12 13 14 15
#ifdef CPP_ELEMENTAL
    USE m_elemental
#endif
#ifdef CPP_SCALAPACK
    USE m_chani
#endif
16 17
#ifdef CPP_ELPA
    USE m_elpa
18 19 20
#endif
    IMPLICIT NONE
#ifdef CPP_ELPA
21 22 23 24 25 26 27 28 29 30 31 32 33
  INTEGER,PARAMETER:: diag_elpa=1
#else
  INTEGER,PARAMETER:: diag_elpa=-1
#endif
#ifdef CPP_ELEMENTAL
  INTEGER,PARAMETER:: diag_elemental=2
#else
  INTEGER,PARAMETER:: diag_elemental=-2
#endif
#ifdef CPP_SCALAPACK
  INTEGER,PARAMETER:: diag_scalapack=3
#else
  INTEGER,PARAMETER:: diag_scalapack=-3
34 35 36 37 38
#endif
#ifdef CPP_MAGMA
  INTEGER,PARAMETER:: diag_magma=6
#else
  INTEGER,PARAMETER:: diag_magma=-6
39 40 41 42 43
#endif
  INTEGER,PARAMETER:: diag_lapack=4
  INTEGER,PARAMETER:: diag_lapack2=5
CONTAINS
  SUBROUTINE eigen_diag(jsp,eig_id,it,atoms,dimension,matsize,mpi, n_rank,n_size,ne,nk,lapw,input,nred,sub_comm,&
44
       sym,matind,kveclo, noco,cell,bkpt,el,jij,l_wu,oneD,td,ud, eig,ne_found,hamOvlp,zMat,realdata)
45 46 47 48 49
    USE m_zsymsecloc
    USE m_aline
    USE m_alinemuff
    USE m_types
    USE m_franza
50 51 52
#ifdef CPP_MAGMA
    USE m_magma
#endif
53 54 55 56 57 58
#ifdef CPP_ELPA
    USE m_elpa
#endif
#ifdef CPP_SCALAPACK
    USE m_chani
#endif
59
#ifdef CPP_ELEMENTAL
60 61
    USE m_elemental
#endif
62
    IMPLICIT NONE
63 64 65
#ifdef CPP_MPI    
    include 'mpif.h'
#endif
66 67 68 69 70 71 72 73 74 75 76 77
    TYPE(t_mpi),INTENT(IN)        :: mpi
    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_jij),INTENT(IN)        :: jij
    TYPE(t_sym),INTENT(IN)        :: sym
    TYPE(t_cell),INTENT(IN)       :: cell
    TYPE(t_atoms),INTENT(IN)      :: atoms
    TYPE(t_lapw),INTENT(INOUT)    :: lapw !might be modified in aline
    TYPE(t_hamOvlp),INTENT(INOUT) :: hamOvlp
    TYPE(t_zMat),INTENT(INOUT)    :: zMat
78 79 80
    INTEGER, INTENT(IN) :: jsp,eig_id,it,matsize 
    INTEGER, INTENT(IN) :: n_rank,n_size  ,nk   ,nred,sub_comm
    INTEGER, INTENT(IN) :: matind(dimension%nbasfcn,2),kveclo(atoms%nlotot)
81 82
    INTEGER,INTENT(IN)  :: ne
    INTEGER,INTENT(OUT) :: ne_found
83 84 85 86 87 88 89 90
    REAL,INTENT(IN)     :: el(:,:,:)
    LOGICAL, INTENT(IN) :: l_wu   
    REAL,INTENT(INOUT)  :: bkpt(3)
    TYPE(t_tlmplm),INTENT(IN) :: td
    TYPE(t_usdus),INTENT(IN)  :: ud

    REAL,INTENT(OUT) :: eig(:)

91
    LOGICAL,OPTIONAL,INTENT(IN) :: realdata
92 93 94 95

    !Locals
    REAL :: time1
    REAL :: time2
96
    INTEGER :: ndim,err,n,nn,i,ndim1
97
    LOGICAL :: parallel
98
    CHARACTER(len=20)::f
99
    LOGICAL :: l_real
100 101

    l_real=hamOvlp%l_real
102
    if (present(realdata)) l_real=realdata
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124

#if 1==2
    !This is only needed for debugging
    print *,n_rank,lapw%nmat
    write(f,'(a,i0)') "a.",n_rank
    open(99,file=f)
    write(f,'(a,i0)') "b.",n_rank
    open(98,file=f)
    i=0
    DO n=n_rank+1,lapw%nmat,n_size
       DO nn=1,n
          i=i+1
          write(99,'(2(i10,1x),4(f15.4,1x))') n,nn,a(i)
          write(98,'(2(i10,1x),4(f15.4,1x))') n,nn,b(i)
       ENDDO
    ENDDO
    CALL MPI_BARRIER(MPI_COMM_WORLD,err)
    close(99)
     close(98)
    STOP 'DEBUG'
#endif
    
125 126 127 128 129
    !
    !  The array z will contain the eigenvectors and is allocated now
    !
    IF (n_size.NE.1) THEN
       ndim = CEILING(real(dimension%neigd)/n_size)
130 131
       ndim1=lapw%nmat
     ELSE
132
       ndim = dimension%neigd
133
       ndim1=dimension%nbasfcn   
134
    ENDIF
135
    if (l_real) THEN
136
       ALLOCATE ( zMat%z_r(ndim1,ndim), STAT = err )
137
    ELSE
138
       ALLOCATE ( zMat%z_c(ndim1,ndim), STAT = err )
139
    END if
140 141 142 143

    zMat%nbasfcn = ndim1
    zMat%nbands = ndim
    zMat%l_real = l_real
144
    
145 146 147 148 149
    IF (err.NE.0) THEN
       WRITE (*,*) 'eigen: error during allocation of the'
       WRITE (*,*) 'eigenvecs',err,'  size: ',dimension%nbasfcn*ndim
       CALL juDFT_error("eigen: Error during allocation of the eigenvecs",calledby ="eigen")
    ENDIF
150
    if (l_real) THEN
151
       zMat%z_r=0.0
152
    else
153
       zMat%z_c=0.0
154
    endif
155
   
156 157 158 159 160
    !l_wu selects a full diagonalization step or a direct call of aline with a subspace diagonalization only
    IF (.NOT.l_wu) THEN
       CALL timestart("Diagonalization")
       !Select the solver
       parallel=(n_size>1)
161
       ne_found=ne
162 163 164
       SELECT CASE (priv_select_solver(parallel))
#ifdef CPP_ELPA
       CASE (diag_elpa)
165
          CALL elpa(lapw%nmat,n,SUB_COMM,eig,ne_found,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c)
166 167 168 169
#endif
#ifdef CPP_ELEMENTAL
       CASE (diag_elemental)
          IF (it==1) THEN !switch between direct solver and iterative solver
170
             CALL elemental(lapw%nmat,dimension%nbasfcn/n_size,SUB_COMM,eig,ne_found,1,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c)
171
          ELSE
172
             CALL elemental(lapw%nmat,,dimension%nbasfcn/n_size,SUB_COMM,eig,ne_found,0,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c)
173 174 175 176 177
          ENDIF

#endif
#ifdef CPP_SCALAPACK
       CASE (diag_scalapack)
178
          CALL chani(lapw%nmat,dimension%nbasfcn/n_size,ndim, n_rank,n_size,SUB_COMM,mpi%mpi_comm,eig,ne_found,hamOvlp,zMat)
179 180 181
#endif
#ifdef CPP_MAGMA
       CASE (diag_magma)
182
          CALL magma_diag(lapw%nmat,eig,ne_found,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c)
183 184 185 186
#endif
       CASE (diag_lapack2)
          if (noco%l_ss) call juDFT_error("zsymsecloc not tested with noco%l_ss")
          if (input%gw>1) call juDFT_error("zsymsecloc not tested with input%gw>1")
187 188
          IF (l_real) THEN
          CALL zsymsecloc(jsp,input,lapw,bkpt,atoms,kveclo, sym,cell, dimension,matsize,ndim,&
189
                jij,matind,nred,eig,ne_found,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r)
190
       else
191
          CALL zsymsecloc(jsp,input,lapw,bkpt,atoms,kveclo, sym,cell, dimension,matsize,ndim,&
192
               jij,matind,nred,eig,ne_found,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c)
193
       endif
194
       CASE (diag_lapack)
195 196
          CALL franza(dimension%nbasfcn,ndim, lapw%nmat,(sym%l_zref.AND.(atoms%nlotot.EQ.0)),&
                      jij%l_j,matind,nred,input%gw,eig,ne_found,hamOvlp,zMat)
197 198 199 200 201 202 203 204 205 206
       CASE DEFAULT
          !This should only happen if you select a solver by hand which was not compiled against
          print*, "You selected a diagonalization scheme without compiling for it"
          CALL priv_solver_error(priv_select_solver(parallel),parallel)
       END SELECT
       CALL timestop("Diagonalization")
       !
    ELSE
       call timestart("aline")
       CALL aline(eig_id,nk,atoms,dimension,sym,cell,input,jsp,el,&
207
            ud,lapw,td,noco,oneD,bkpt,eig,ne_found,zMat,hamOvlp%a_r,hamOvlp%b_r,hamOvlp%a_c,hamOvlp%b_c)
208 209 210 211 212 213 214 215
       call timestop("aline")
    ENDIF
    !--->         SECOND VARIATION STEP
    IF (input%secvar) THEN
       !--->           compute and diagonalize the second variation
       !--->           hamiltonian
       call timestart("second variation diagonalization")

216
       CALL aline_muff(atoms,dimension,sym, cell, jsp,ne_found, ud,td, bkpt,lapw,eig,zMat%z_r,zMat%z_c,realdata)
217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233
       call timestop("second variation diagonalization")
    END IF
  END SUBROUTINE eigen_diag

  FUNCTION priv_select_solver(parallel) result(diag_solver)
    LOGICAL,INTENT(IN):: parallel
    INTEGER           :: diag_solver
    diag_solver=-99

    !Determine the default solver
    IF (parallel) THEN
#ifdef CPP_ELPA
       diag_solver=diag_elpa
#else
       diag_solver=diag_scalapack
#endif
    ELSE
234 235 236
#ifdef CPP_MAGMA
       diag_solver=diag_magma
#else
237
       diag_solver=diag_lapack
238
#endif
239 240 241 242 243 244 245 246
    ENDIF

    !check if a special solver was requested
    IF (juDFT_was_argument("-elpa")) diag_solver=diag_elpa
    IF (juDFT_was_argument("-scalapack")) diag_solver=diag_scalapack
    IF (juDFT_was_argument("-elemental")) diag_solver=diag_elemental
    IF (juDFT_was_argument("-lapack")) diag_solver=diag_lapack
    IF (juDFT_was_argument("-lapack2")) diag_solver=diag_lapack2
247
    IF (juDFT_was_argument("-magma")) diag_solver=diag_magma
248 249 250 251 252 253 254 255 256 257 258
    !Check if solver is possible
    if (diag_solver<0) call priv_solver_error(diag_solver,parallel)
    if (any((/diag_lapack,diag_lapack2/)==diag_solver).and.parallel) call priv_solver_error(diag_solver,parallel)
    if (any((/diag_elpa,diag_elemental,diag_scalapack/)==diag_solver).and..not.parallel) call priv_solver_error(diag_solver,parallel)

  END FUNCTION priv_select_solver

  SUBROUTINE priv_solver_error(diag_solver,parallel)
    IMPLICIT NONE
    INTEGER,INTENT(IN):: diag_solver
    LOGICAL,INTENT(IN)::parallel
259
    print *,diag_solver,parallel
260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286
    SELECT CASE(diag_solver)
    CASE (diag_elpa)
       IF (parallel) THEN
          CALL juDFT_error("You did not compile with the ELPA solver and hence can not use it")
       ELSE
          CALL juDFT_error("The ELPA solver can not be used in serial")
       ENDIF
    CASE (diag_elemental)
       IF (parallel) THEN
          CALL juDFT_error("You did not compile with the ELEMENTAL solver and hence can not use it")
       ELSE
          CALL juDFT_error("The ELEMENTAL solver can not be used in serial")
       ENDIF
    CASE (diag_scalapack)
       IF (parallel) THEN
          CALL juDFT_error("You did not compile with the SCALAPACK solver and hence can not use it")
       ELSE
          CALL juDFT_error("The SCALAPACK solver can not be used in serial")
       ENDIF
    CASE (diag_lapack)
       IF (parallel) THEN
          CALL juDFT_error("The LAPACK solver can not be used in parallel")
       ENDIF
    CASE (diag_lapack2)
       IF (parallel) THEN
          CALL juDFT_error("The LAPACK2 solver can not be used in parallel")
       ENDIF  
287 288
    CASE (diag_magma)
       CALL juDFT_error("You have not compiled with MAGMA support")
289 290 291 292 293 294 295
    CASE DEFAULT
       CALL judft_error("You have selected an unkown eigensolver")
    END SELECT
  END SUBROUTINE priv_solver_error


END MODULE m_eigen_diag