eigen_diag.F90 10.3 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
#endif
    IMPLICIT NONE
Daniel Wortmann's avatar
Daniel Wortmann committed
20
    PRIVATE
21
#ifdef CPP_ELPA
22 23 24 25 26 27 28 29 30 31 32 33 34
  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
35 36 37 38 39
#endif
#ifdef CPP_MAGMA
  INTEGER,PARAMETER:: diag_magma=6
#else
  INTEGER,PARAMETER:: diag_magma=-6
40 41 42
#endif
  INTEGER,PARAMETER:: diag_lapack=4
  INTEGER,PARAMETER:: diag_lapack2=5
Daniel Wortmann's avatar
Daniel Wortmann committed
43
  PUBLIC eigen_diag,parallel_solver_available
44
CONTAINS
Daniel Wortmann's avatar
Daniel Wortmann committed
45 46 47 48 49

  LOGICAL FUNCTION parallel_solver_available()
    parallel_solver_available=any((/diag_elpa,diag_elemental,diag_scalapack/)>0)
  END FUNCTION parallel_solver_available

50
  SUBROUTINE eigen_diag(jsp,eig_id,it,atoms,dimension,matsize,mpi, n_rank,n_size,ne,nk,lapw,input,nred,sub_comm,&
Daniel Wortmann's avatar
Daniel Wortmann committed
51
       sym,l_zref,matind,kveclo, noco,cell,bkpt,el,jij,l_wu,oneD,td,ud, eig,ne_found,hamOvlp,zMat,realdata)
52 53 54 55 56
    USE m_zsymsecloc
    USE m_aline
    USE m_alinemuff
    USE m_types
    USE m_franza
57 58 59
#ifdef CPP_MAGMA
    USE m_magma
#endif
60 61 62 63 64 65
#ifdef CPP_ELPA
    USE m_elpa
#endif
#ifdef CPP_SCALAPACK
    USE m_chani
#endif
66
#ifdef CPP_ELEMENTAL
67 68
    USE m_elemental
#endif
69
    IMPLICIT NONE
70 71 72
#ifdef CPP_MPI    
    include 'mpif.h'
#endif
73 74 75 76 77 78 79 80 81 82 83 84
    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
85 86 87
    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)
88 89
    INTEGER,INTENT(IN)  :: ne
    INTEGER,INTENT(OUT) :: ne_found
90
    REAL,INTENT(IN)     :: el(:,:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
91
    LOGICAL, INTENT(IN) :: l_wu,l_zref
92 93 94 95 96 97
    REAL,INTENT(INOUT)  :: bkpt(3)
    TYPE(t_tlmplm),INTENT(IN) :: td
    TYPE(t_usdus),INTENT(IN)  :: ud

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

98
    LOGICAL,OPTIONAL,INTENT(IN) :: realdata
99 100 101 102

    !Locals
    REAL :: time1
    REAL :: time2
103
    INTEGER :: ndim,err,n,nn,i,ndim1
104
    LOGICAL :: parallel
105
    CHARACTER(len=20)::f
106
    LOGICAL :: l_real
107 108

    l_real=hamOvlp%l_real
109
    if (present(realdata)) l_real=realdata
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131

#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
    
132 133 134 135 136
    !
    !  The array z will contain the eigenvectors and is allocated now
    !
    IF (n_size.NE.1) THEN
       ndim = CEILING(real(dimension%neigd)/n_size)
137 138
       ndim1=lapw%nmat
     ELSE
139
       ndim = dimension%neigd
140
       ndim1=dimension%nbasfcn   
141
    ENDIF
142
    if (l_real) THEN
143
       ALLOCATE ( zMat%z_r(ndim1,ndim), STAT = err )
144
    ELSE
145
       ALLOCATE ( zMat%z_c(ndim1,ndim), STAT = err )
146
    END if
147 148 149 150

    zMat%nbasfcn = ndim1
    zMat%nbands = ndim
    zMat%l_real = l_real
151
    
152 153 154 155 156
    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
157
    if (l_real) THEN
158
       zMat%z_r=0.0
159
    else
160
       zMat%z_c=0.0
161
    endif
162
   
163 164 165 166 167
    !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)
168
       ne_found=ne
169 170 171
       SELECT CASE (priv_select_solver(parallel))
#ifdef CPP_ELPA
       CASE (diag_elpa)
172
          IF (hamovlp%l_real) THEN
173
              CALL elpa_diag(lapw%nmat,n,SUB_COMM,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r,eig,ne_found)
174
          ELSE
175
              CALL elpa_diag(lapw%nmat,n,SUB_COMM,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c,eig,ne_found)
176
          ENDIF
177 178 179 180
#endif
#ifdef CPP_ELEMENTAL
       CASE (diag_elemental)
          IF (it==1) THEN !switch between direct solver and iterative solver
181
             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)
182
          ELSE
183
             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)
184 185 186 187 188
          ENDIF

#endif
#ifdef CPP_SCALAPACK
       CASE (diag_scalapack)
189
          CALL chani(lapw%nmat,dimension%nbasfcn/n_size,ndim, n_rank,n_size,SUB_COMM,mpi%mpi_comm,eig,ne_found,hamOvlp,zMat)
190 191 192
#endif
#ifdef CPP_MAGMA
       CASE (diag_magma)
193
          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)
194 195 196 197
#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")
198
          IF (l_real) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
199
          CALL zsymsecloc(jsp,input,lapw,bkpt,atoms,kveclo, sym,l_zref,cell, dimension,matsize,ndim,&
200
                jij,matind,nred,eig,ne_found,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r)
201
       else
Daniel Wortmann's avatar
Daniel Wortmann committed
202
          CALL zsymsecloc(jsp,input,lapw,bkpt,atoms,kveclo, sym,l_zref,cell, dimension,matsize,ndim,&
203
               jij,matind,nred,eig,ne_found,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c)
204
       endif
205
       CASE (diag_lapack)
Daniel Wortmann's avatar
Daniel Wortmann committed
206
          CALL franza(dimension%nbasfcn,ndim, lapw%nmat,(l_zref.AND.(atoms%nlotot.EQ.0)),&
207
                      jij%l_j,matind,nred,input%gw,eig,ne_found,hamOvlp,zMat)
208 209 210 211 212 213 214 215 216 217
       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,&
218
            ud,lapw,td,noco,oneD,bkpt,eig,ne_found,zMat,hamOvlp%a_r,hamOvlp%b_r,hamOvlp%a_c,hamOvlp%b_c)
219 220 221 222 223 224 225 226
       call timestop("aline")
    ENDIF
    !--->         SECOND VARIATION STEP
    IF (input%secvar) THEN
       !--->           compute and diagonalize the second variation
       !--->           hamiltonian
       call timestart("second variation diagonalization")

227
       CALL aline_muff(atoms,dimension,sym, cell, jsp,ne_found, ud,td, bkpt,lapw,eig,zMat%z_r,zMat%z_c,realdata)
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244
       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
245 246 247
#ifdef CPP_MAGMA
       diag_solver=diag_magma
#else
248
       diag_solver=diag_lapack
249
#endif
250 251 252 253 254 255 256 257
    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
258
    IF (juDFT_was_argument("-magma")) diag_solver=diag_magma
259 260 261 262 263 264 265
    !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

Daniel Wortmann's avatar
Daniel Wortmann committed
266

267 268 269 270
  SUBROUTINE priv_solver_error(diag_solver,parallel)
    IMPLICIT NONE
    INTEGER,INTENT(IN):: diag_solver
    LOGICAL,INTENT(IN)::parallel
271
    print *,diag_solver,parallel
272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298
    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  
299 300
    CASE (diag_magma)
       CALL juDFT_error("You have not compiled with MAGMA support")
301 302 303 304 305 306 307
    CASE DEFAULT
       CALL judft_error("You have selected an unkown eigensolver")
    END SELECT
  END SUBROUTINE priv_solver_error


END MODULE m_eigen_diag