eigen_diag.F90 8.44 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
    IMPLICIT NONE
Daniel Wortmann's avatar
Daniel Wortmann committed
11
    PRIVATE
12
#ifdef CPP_ELPA
13 14 15 16 17 18 19 20 21 22 23 24 25
  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
26 27 28 29 30
#endif
#ifdef CPP_MAGMA
  INTEGER,PARAMETER:: diag_magma=6
#else
  INTEGER,PARAMETER:: diag_magma=-6
31
#endif
Daniel Wortmann's avatar
Daniel Wortmann committed
32
#ifdef CPP_CHASE
33
  INTEGER,PARAMETER:: diag_chase=7
Daniel Wortmann's avatar
Daniel Wortmann committed
34 35 36
#else
  INTEGER,PARAMETER:: diag_chase=-7
#endif
37
#ifdef CPP_CUSOLVER
Daniel Wortmann's avatar
Daniel Wortmann committed
38 39 40 41 42 43
  INTEGER,PARAMETER:: diag_cusolver=8
#else
  INTEGER,PARAMETER:: diag_cusolver=-8
#endif
  
  INTEGER,PARAMETER:: diag_lapack=4
44
#ifdef CPP_ELPA_ONENODE
45
  INTEGER,PARAMETER:: diag_elpa_1node=14
46 47 48
#else
  INTEGER,PARAMETER:: diag_elpa_1node=-14
#endif 
Daniel Wortmann's avatar
Daniel Wortmann committed
49
  INTEGER,PARAMETER:: diag_debugout=99
Daniel Wortmann's avatar
Daniel Wortmann committed
50
  PUBLIC eigen_diag,parallel_solver_available
51
CONTAINS
Daniel Wortmann's avatar
Daniel Wortmann committed
52 53 54 55 56

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

57
  SUBROUTINE eigen_diag(mpi,hmat,smat,ikpt,jsp,iter,ne,eig,ev)
Daniel Wortmann's avatar
Daniel Wortmann committed
58
    USE m_lapack_diag
59
    USE m_magma
60
    USE m_elpa
61
    USE m_elpa_onenode
Daniel Wortmann's avatar
Daniel Wortmann committed
62
    USE m_scalapack
63
    USE m_elemental
64 65
    USE m_chase_diag
    USE m_types_mpimat
Daniel Wortmann's avatar
Daniel Wortmann committed
66
    USE m_types_gpumat
67
    USE m_matrix_copy
Daniel Wortmann's avatar
Daniel Wortmann committed
68
    USE m_cusolver_diag
69
    IMPLICIT NONE
70
#ifdef CPP_MPI
71 72
    include 'mpif.h'
#endif
73
    TYPE(t_mpi),               INTENT(IN)    :: mpi
74 75 76 77
    CLASS(t_mat),              INTENT(INOUT) :: smat,hmat
    CLASS(t_mat), ALLOCATABLE, INTENT(OUT)   :: ev
    INTEGER,                   INTENT(IN)    :: ikpt
    INTEGER,                   INTENT(IN)    :: jsp
78
    INTEGER,                   INTENT(IN)    :: iter
79 80
    INTEGER,                   INTENT(INOUT) :: ne
    REAL,                      INTENT(OUT)   :: eig(:)
81

82
    !Locals  
83
    LOGICAL :: parallel
84 85
    !For check-mode
    TYPE(t_mat) :: s_store,h_store
86
    
87
         
Daniel Wortmann's avatar
Daniel Wortmann committed
88 89 90 91 92 93
    SELECT TYPE(smat)
    CLASS IS (t_mpimat)
       parallel=.TRUE.
    CLASS default
       parallel=.FALSE.
    END SELECT
94 95

    !Create a copy of the matrix if in test mode
96
    IF (TRIM(judft_string_for_argument("-diag"))=="test") THEN
97 98 99 100 101 102 103 104 105 106 107 108
       SELECT TYPE(hmat)
        CLASS IS (t_mpimat)
          CALL s_store%init(hmat%l_real,hmat%global_size1,hmat%global_size2)
          CALL h_store%init(hmat%l_real,hmat%global_size1,hmat%global_size2)
       CLASS default
          CALL s_store%init(hmat%l_real,hmat%matsize1,hmat%matsize2)
          CALL h_store%init(hmat%l_real,hmat%matsize1,hmat%matsize2)
       END SELECT
       CALL matrix_copy(smat,s_store)
       CALL matrix_copy(hmat,h_store)
    END IF
       
Daniel Wortmann's avatar
Daniel Wortmann committed
109 110 111 112
    CALL timestart("Diagonalization")
    !Select the solver
    SELECT CASE (priv_select_solver(parallel))
    CASE (diag_elpa)
Uliana Alekseeva's avatar
Uliana Alekseeva committed
113
       CALL elpa_diag(hmat,smat,ne,eig,ev)
114 115
    CASE (diag_elpa_1node)
       CALL elpa_diag_onenode(hmat,smat,ne,eig,ev)
Daniel Wortmann's avatar
Daniel Wortmann committed
116 117 118 119 120 121
    CASE (diag_elemental)
       !CALL ELEMENTAL(hmat,smat,ne,eig,ev)
    CASE (diag_scalapack)
       CALL scalapack(hmat,smat,ne,eig,ev)
    CASE (diag_magma)
       !CALL magma_diag(hmat,smat,ne,eig,ev)
Daniel Wortmann's avatar
Daniel Wortmann committed
122 123
    CASE (diag_cusolver)
       CALL cusolver_diag(hmat,smat,ne,eig,ev)
Daniel Wortmann's avatar
Daniel Wortmann committed
124 125
    CASE (diag_lapack)
       CALL lapack_diag(hmat,smat,ne,eig,ev)
126
    CASE (diag_chase)
127
       CALL chase_diag(hmat,smat,ikpt,jsp,iter,ne,eig,ev)
Daniel Wortmann's avatar
Daniel Wortmann committed
128 129 130
    CASE (diag_debugout)
       CALL priv_debug_out(smat,hmat)
    END SELECT
131 132

    !Create test the solution
133
    IF (TRIM(judft_string_for_argument("-diag"))=="test") THEN
134 135 136
       CALL priv_test_solution(mpi,ne,s_store,h_store,eig,ev)
       CALL judft_error("Diagonalization tested")
    END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
137
    CALL timestop("Diagonalization")
138 139 140
    !
  END SUBROUTINE eigen_diag

Daniel Wortmann's avatar
Daniel Wortmann committed
141

142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185
  SUBROUTINE priv_test_solution(mpi,ne,s_store,h_store,eig1,ev)
    USE m_types
    USE m_lapack_diag
    USE m_matrix_copy
    IMPLICIT NONE
    TYPE(t_mpi),INTENT(in):: mpi
    INTEGER,INTENT(INOUT) :: ne
    TYPE(t_mat)           :: s_store,h_store
    REAL,INTENT(in)       :: eig1(:)
    CLASS(t_mat)          :: ev

    
    REAL,ALLOCATABLE::eig2(:)
    TYPE(t_mat)     ::ev1
    CLASS(t_mat),ALLOCATABLE    ::ev2
    INTEGER         :: i,irank
    

    ALLOCATE(eig2(ne))
    CALL ev1%init(s_store%l_real,s_store%matsize1,s_store%matsize2)
    CALL matrix_copy(ev,ev1)
    
    IF (mpi%irank==0) THEN
       CALL lapack_diag(h_store,s_store,ne,eig2,ev2)

       OPEN(99,file="diag.compare")
       WRITE(99,*) "Eigenvalues"
       DO i=1,ne
          WRITE(99,*) i,eig1(i),eig2(i)
       ENDDO
       WRITE(99,*) "Eigenvectors"
       DO i=1,ne
          IF (ev1%l_real) THEN
             WRITE(99,"(i0,20(1x,f10.5))") i,ev1%data_r(1:10,i)
             WRITE(99,"(i0,20(1x,f10.5))") i,ev2%data_r(1:10,i)
          ELSE
             WRITE(99,"(i0,20(1x,f10.5))") i,ev1%data_c(1:10,i)
             WRITE(99,"(i0,20(1x,f10.5))") i,ev2%data_c(1:10,i)
          END IF
       ENDDO
       CLOSE(99)
    END IF
  END SUBROUTINE priv_test_solution
  
Daniel Wortmann's avatar
Daniel Wortmann committed
186
  SUBROUTINE priv_debug_out(smat,hmat)
187 188 189
    USE m_types
    use m_judft
    IMPLICIT NONE
Daniel Wortmann's avatar
Daniel Wortmann committed
190
    CLASS(t_mat),INTENT(INOUT) :: hmat,smat
191
    !small subroutine that does only wite the matrix to a file
Daniel Wortmann's avatar
Daniel Wortmann committed
192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212
    INTEGER:: i,ii,irank,ierr
    CHARACTER(len=20)::filename
#ifdef CPP_MPI
    INCLUDE 'mpif.h'
    CALL MPI_COMM_RANK(MPI_COMM_WORLD,irank,ierr)
#else
    irank=0
#endif
    WRITE(filename,"(a,i0)") "hmat",irank
    OPEN(999,file=TRIM(filename))
    WRITE(filename,"(a,i0)") "smat",irank
    OPEN(998,file=TRIM(filename))
    DO i=1,hmat%matsize2
       DO ii=1,hmat%matsize1
          IF (hmat%l_real) THEN
             WRITE(999,"(2i6,f15.6)") ii,i,hmat%data_r(ii,i)
             WRITE(998,"(2i6,f15.6)") ii,i,smat%data_r(ii,i)
          ELSE
             WRITE(999,"(2i6,2f15.6)") ii,i,hmat%data_c(ii,i)
             WRITE(998,"(2i6,2f15.6)") ii,i,smat%data_c(ii,i)
          ENDIF
213
       END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
214
    ENDDO
215 216
    CLOSE(999)
    CLOSE(998)
Daniel Wortmann's avatar
Daniel Wortmann committed
217 218 219
#ifdef CPP_MPI
    CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
#endif
220 221 222
    CALL judft_error("STOP in eigen_diag:debug_diag")
  END SUBROUTINE priv_debug_out

Daniel Wortmann's avatar
Daniel Wortmann committed
223
  FUNCTION priv_select_solver(parallel) RESULT(diag_solver)
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239
    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
       diag_solver=diag_lapack
    ENDIF

    !check if a special solver was requested
240 241 242 243 244 245 246 247 248
    IF (TRIM(juDFT_string_for_argument("-diag"))=="elpa")       diag_solver=diag_elpa
    IF (TRIM(juDFT_string_for_argument("-diag"))=="elpa_1node") diag_solver=diag_elpa_1node
    IF (trim(juDFT_string_for_argument("-diag"))=="scalapack")  diag_solver=diag_scalapack
    IF (trim(juDFT_string_for_argument("-diag"))=="elemental")  diag_solver=diag_elemental
    IF (trim(juDFT_string_for_argument("-diag"))=="lapack")     diag_solver=diag_lapack
    IF (trim(juDFT_string_for_argument("-diag"))=="magma")      diag_solver=diag_magma
    IF (trim(juDFT_string_for_argument("-diag"))=="chase")      diag_solver=diag_chase
    IF (trim(juDFT_string_for_argument("-diag"))=="cusolver")   diag_solver=diag_cusolver
    IF (trim(juDFT_string_for_argument("-diag"))=="debugout")   diag_solver=diag_debugout
Daniel Wortmann's avatar
Daniel Wortmann committed
249
    
250
    !Check if solver is possible
Daniel Wortmann's avatar
Daniel Wortmann committed
251 252 253
    IF (diag_solver<0)  CALL juDFT_error("You selected a solver for the eigenvalue problem that is not available",hint="You most probably did not provide the appropriate libraries for compilation/linking")
    IF (ANY((/diag_lapack,diag_magma,diag_cusolver/)==diag_solver).AND.parallel) CALL judft_error("You selected an eigensolver that does not support distributed memory parallism",hint="Try scalapack,elpa or another supported solver for parallel matrices")
    IF (ANY((/diag_elpa,diag_elemental,diag_scalapack/)==diag_solver).AND..NOT.parallel) CALL judft_error("You selected an eigensolver for matrices that are memory distributed",hint="Try lapack, cusolver or another supported solver for non-distributed matrices")
254

255

256 257
  END FUNCTION priv_select_solver

Daniel Wortmann's avatar
Daniel Wortmann committed
258

Daniel Wortmann's avatar
Daniel Wortmann committed
259
 
260 261

END MODULE m_eigen_diag