eigen_diag.F90 11.5 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
43
   INTEGER,PARAMETER:: diag_debugout=99
Daniel Wortmann's avatar
Daniel Wortmann committed
44
  PUBLIC eigen_diag,parallel_solver_available
45
CONTAINS
Daniel Wortmann's avatar
Daniel Wortmann committed
46 47 48 49 50

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

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

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

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

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

    l_real=hamOvlp%l_real
108
    if (present(realdata)) l_real=realdata
109 110 111 112

#if 1==2
    !This is only needed for debugging
    print *,n_rank,lapw%nmat
113 114
    print *,"SR:",n_size,n_rank
    print *,mpi
115 116 117 118 119 120 121 122
    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
123 124
          write(99,'(2(i10,1x),4(f15.8,1x))') n,nn,hamOvlp%a_c(i)
          write(98,'(2(i10,1x),4(f15.8,1x))') n,nn,hamOvlp%b_c(i)
125 126 127 128 129 130 131
       ENDDO
    ENDDO
    CALL MPI_BARRIER(MPI_COMM_WORLD,err)
    close(99)
     close(98)
#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 173 174
          CALL MPI_COMM_RANK(sub_comm,n,err)
          write(*,*) "DIAG:",mpi%irank,sub_comm,n
          write(*,*) "ELPA:",mpi%irank,lapw%nmat,ne_found,size(eig),ndim,ndim1
175
          IF (hamovlp%l_real) THEN
176
              CALL elpa_diag(lapw%nmat,SUB_COMM,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r,eig,ne_found)
177
          ELSE
178
              CALL elpa_diag(lapw%nmat,SUB_COMM,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c,eig,ne_found)
179
          ENDIF
180 181 182 183
#endif
#ifdef CPP_ELEMENTAL
       CASE (diag_elemental)
          IF (it==1) THEN !switch between direct solver and iterative solver
184
             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)
185
          ELSE
186
             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)
187 188 189 190 191
          ENDIF

#endif
#ifdef CPP_SCALAPACK
       CASE (diag_scalapack)
192
          CALL chani(lapw%nmat,ndim, n_rank,n_size,SUB_COMM,mpi%mpi_comm,eig,ne_found,hamOvlp,zMat)
193 194 195
#endif
#ifdef CPP_MAGMA
       CASE (diag_magma)
196 197 198 199 200 201
          if (l_real) THEN
             call juDFT_error("REAL diagonalization not implemented in magma.F90")
          else
             print *,"Start magma_diag"
             CALL magma_diag(lapw%nmat,eig,ne_found,a_c=hamOvlp%a_c,b_c=hamOvlp%b_c,z_c=zMat%z_c)
          endif
202 203 204 205
#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")
206
          IF (l_real) THEN
207
#ifndef __PGI
Daniel Wortmann's avatar
Daniel Wortmann committed
208
          CALL zsymsecloc(jsp,input,lapw,bkpt,atoms,kveclo, sym,l_zref,cell, dimension,matsize,ndim,&
209
                jij,matind,nred,eig,ne_found,hamOvlp%a_r,hamOvlp%b_r,zMat%z_r)
210
       else
Daniel Wortmann's avatar
Daniel Wortmann committed
211
          CALL zsymsecloc(jsp,input,lapw,bkpt,atoms,kveclo, sym,l_zref,cell, dimension,matsize,ndim,&
212
               jij,matind,nred,eig,ne_found,hamOvlp%a_c,hamOvlp%b_c,zMat%z_c)
213
#endif
214
       endif
215
       CASE (diag_lapack)
Daniel Wortmann's avatar
Daniel Wortmann committed
216
          CALL franza(dimension%nbasfcn,ndim, lapw%nmat,(l_zref.AND.(atoms%nlotot.EQ.0)),&
217 218 219
               jij%l_j,matind,nred,input%gw,eig,ne_found,hamOvlp,zMat)
           CASE (diag_debugout)
          CALL priv_debug_out(hamovlp)
220 221 222 223 224 225 226 227 228 229
       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,&
230
            ud,lapw,td,noco,oneD,bkpt,eig,ne_found,zMat,hamOvlp%a_r,hamOvlp%b_r,hamOvlp%a_c,hamOvlp%b_c)
231 232 233 234 235 236 237 238
       call timestop("aline")
    ENDIF
    !--->         SECOND VARIATION STEP
    IF (input%secvar) THEN
       !--->           compute and diagonalize the second variation
       !--->           hamiltonian
       call timestart("second variation diagonalization")

239
       CALL aline_muff(atoms,dimension,sym, cell, jsp,ne_found, ud,td, bkpt,lapw,eig,zMat%z_r,zMat%z_c,realdata)
240 241 242 243
       call timestop("second variation diagonalization")
    END IF
  END SUBROUTINE eigen_diag

244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268
  SUBROUTINE priv_debug_out(hamovlp)
    USE m_types
    use m_judft
    IMPLICIT NONE
    TYPE(t_hamovlp),INTENT(INOUT) :: hamovlp
    !small subroutine that does only wite the matrix to a file
    INTEGER:: i,ii
    OPEN(999,file="hmat")
    OPEN(998,file="smat")
    IF (hamovlp%l_real) THEN
       DO i=1,SIZE(hamovlp%a_r)
          WRITE(999,"(i6,f15.6)")  i,hamovlp%a_r(i)
          WRITE(998,"(i6,f15.6)")  i,hamovlp%b_r(i)
       END DO
       ELSE
       DO i=1,SIZE(hamovlp%a_c)
          WRITE(999,"(i6,2f15.6)") i,hamovlp%a_c(i)
          WRITE(998,"(i6,2f15.6)") i,hamovlp%b_c(i)
       ENDDO
    ENDIF
    CLOSE(999)
    CLOSE(998)
    CALL judft_error("STOP in eigen_diag:debug_diag")
  END SUBROUTINE priv_debug_out

269 270 271 272 273 274 275 276 277 278 279 280 281
  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
282 283 284
#ifdef CPP_MAGMA
       diag_solver=diag_magma
#else
285
       diag_solver=diag_lapack
286
#endif
287 288 289 290 291 292 293 294
    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
295
    IF (juDFT_was_argument("-magma")) diag_solver=diag_magma
296
     IF (juDFT_was_argument("-debugdiag")) diag_solver=diag_debugout
297 298 299 300 301 302 303
    !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
304

305 306 307 308
  SUBROUTINE priv_solver_error(diag_solver,parallel)
    IMPLICIT NONE
    INTEGER,INTENT(IN):: diag_solver
    LOGICAL,INTENT(IN)::parallel
309
    print *,diag_solver,parallel
310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336
    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  
337 338
    CASE (diag_magma)
       CALL juDFT_error("You have not compiled with MAGMA support")
339 340 341 342 343 344 345
    CASE DEFAULT
       CALL judft_error("You have selected an unkown eigensolver")
    END SELECT
  END SUBROUTINE priv_solver_error


END MODULE m_eigen_diag