1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief used for collecting diagonalization schemes available for cp_cfm_type 8!> \note 9!> first version : only one routine right now 10!> \author Joost VandeVondele (2003-09) 11! ************************************************************************************************** 12MODULE cp_cfm_diag 13 USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose, & 14 cp_cfm_gemm, & 15 cp_cfm_column_scale, & 16 cp_cfm_scale, & 17 cp_cfm_triangular_invert, & 18 cp_cfm_triangular_multiply 19 USE cp_cfm_types, ONLY: cp_cfm_get_info, & 20 cp_cfm_set_element, & 21 cp_cfm_to_cfm, & 22 cp_cfm_type 23 USE kinds, ONLY: dp 24#if defined (__HAS_IEEE_EXCEPTIONS) 25 USE ieee_exceptions, ONLY: ieee_get_halting_mode, & 26 ieee_set_halting_mode, & 27 ieee_all 28#endif 29 30#include "../base/base_uses.f90" 31 32 IMPLICIT NONE 33 PRIVATE 34 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_diag' 35 36 PUBLIC :: cp_cfm_heevd, cp_cfm_geeig, cp_cfm_geeig_canon 37 38CONTAINS 39 40! ************************************************************************************************** 41!> \brief Perform a diagonalisation of a complex matrix 42!> \param matrix ... 43!> \param eigenvectors ... 44!> \param eigenvalues ... 45!> \par History 46!> - (De)Allocation checks updated (15.02.2011,MK) 47!> \author Joost VandeVondele 48! ************************************************************************************************** 49 SUBROUTINE cp_cfm_heevd(matrix, eigenvectors, eigenvalues) 50 51 TYPE(cp_cfm_type), POINTER :: matrix, eigenvectors 52 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues 53 54 CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd', & 55 routineP = moduleN//':'//routineN 56 57 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work 58 COMPLEX(KIND=dp), DIMENSION(:, :), & 59 POINTER :: m 60 INTEGER :: handle, info, liwork, & 61 lrwork, lwork, n 62 INTEGER, DIMENSION(:), POINTER :: iwork 63 REAL(KIND=dp), DIMENSION(:), POINTER :: rwork 64#if defined(__SCALAPACK) 65 INTEGER, DIMENSION(9) :: descm, descv 66 COMPLEX(KIND=dp), DIMENSION(:, :), & 67 POINTER :: v 68#if defined (__HAS_IEEE_EXCEPTIONS) 69 LOGICAL, DIMENSION(5) :: halt 70#endif 71#endif 72 73 CALL timeset(routineN, handle) 74 75 CPASSERT(ASSOCIATED(matrix)) 76 CPASSERT(ASSOCIATED(eigenvectors)) 77 78 n = matrix%matrix_struct%nrow_global 79 m => matrix%local_data 80 ALLOCATE (iwork(1), rwork(1), work(1)) 81 ! work space query 82 lwork = -1 83 lrwork = -1 84 liwork = -1 85 86#if defined(__SCALAPACK) 87 v => eigenvectors%local_data 88 descm(:) = matrix%matrix_struct%descriptor(:) 89 descv(:) = eigenvectors%matrix_struct%descriptor(:) 90 CALL PZHEEVD('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, & 91 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info) 92 lwork = CEILING(REAL(work(1), KIND=dp)) 93 ! needed to correct for a bug in scalapack, unclear how much the right number is 94 lrwork = CEILING(REAL(rwork(1), KIND=dp)) + 1000000 95 liwork = iwork(1) 96#else 97 CALL ZHEEVD('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), & 98 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info) 99 lwork = CEILING(REAL(work(1), KIND=dp)) 100 lrwork = CEILING(REAL(rwork(1), KIND=dp)) 101 liwork = iwork(1) 102#endif 103 104 DEALLOCATE (iwork, rwork, work) 105 106 ALLOCATE (iwork(liwork)) 107 iwork(:) = 0 108 ALLOCATE (rwork(lrwork)) 109 rwork(:) = 0.0_dp 110 ALLOCATE (work(lwork)) 111 work(:) = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 112 113#if defined(__SCALAPACK) 114! Scalapack takes advantage of IEEE754 exceptions for speedup. 115! Therefore, we disable floating point traps temporarily. 116#if defined (__HAS_IEEE_EXCEPTIONS) 117 CALL ieee_get_halting_mode(IEEE_ALL, halt) 118 CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.) 119#endif 120 121 CALL PZHEEVD('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, & 122 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info) 123 124#if defined (__HAS_IEEE_EXCEPTIONS) 125 CALL ieee_set_halting_mode(IEEE_ALL, halt) 126#endif 127#else 128 CALL ZHEEVD('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), & 129 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info) 130 eigenvectors%local_data = matrix%local_data 131#endif 132 133 IF (info /= 0) & 134 CPABORT("Diagonalisation of a complex matrix failed") 135 DEALLOCATE (iwork, rwork, work) 136 137 CALL timestop(handle) 138 139 END SUBROUTINE cp_cfm_heevd 140 141! ************************************************************************************************** 142!> \brief General Eigenvalue Problem AX = BXE 143!> Single option version: Cholesky decomposition of B 144!> \param amatrix ... 145!> \param bmatrix ... 146!> \param eigenvectors ... 147!> \param eigenvalues ... 148!> \param work ... 149! ************************************************************************************************** 150 SUBROUTINE cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work) 151 152 TYPE(cp_cfm_type), POINTER :: amatrix, bmatrix, eigenvectors 153 REAL(KIND=dp), DIMENSION(:) :: eigenvalues 154 TYPE(cp_cfm_type), POINTER :: work 155 156 CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig', routineP = moduleN//':'//routineN 157 158 INTEGER :: handle, nao, nmo 159 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals 160 161 CALL timeset(routineN, handle) 162 163 CALL cp_cfm_get_info(amatrix, nrow_global=nao) 164 ALLOCATE (evals(nao)) 165 ! Cholesky decompose S=U(T)U 166 CALL cp_cfm_cholesky_decompose(bmatrix) 167 ! Invert to get U^(-1) 168 CALL cp_cfm_triangular_invert(bmatrix) 169 ! Reduce to get U^(-T) * H * U^(-1) 170 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R") 171 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C") 172 ! Diagonalize 173 CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals) 174 ! Restore vectors C = U^(-1) * C* 175 CALL cp_cfm_triangular_multiply(bmatrix, work) 176 nmo = SIZE(eigenvalues) 177 CALL cp_cfm_to_cfm(work, eigenvectors, nmo) 178 eigenvalues(1:nmo) = evals(1:nmo) 179 180 DEALLOCATE (evals) 181 182 CALL timestop(handle) 183 184 END SUBROUTINE cp_cfm_geeig 185 186! ************************************************************************************************** 187!> \brief General Eigenvalue Problem AX = BXE 188!> Use canonical orthogonalization 189!> \param amatrix ... 190!> \param bmatrix ... 191!> \param eigenvectors ... 192!> \param eigenvalues ... 193!> \param work ... 194!> \param epseig ... 195! ************************************************************************************************** 196 SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig) 197 198 TYPE(cp_cfm_type), POINTER :: amatrix, bmatrix, eigenvectors 199 REAL(KIND=dp), DIMENSION(:) :: eigenvalues 200 TYPE(cp_cfm_type), POINTER :: work 201 REAL(KIND=dp), INTENT(IN) :: epseig 202 203 CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig_canon', & 204 routineP = moduleN//':'//routineN 205 COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), & 206 czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 207 208 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals 209 INTEGER :: handle, i, icol, irow, nao, nc, ncol, & 210 nmo, nx 211 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals 212 213 CALL timeset(routineN, handle) 214 215 ! Test sizees 216 CALL cp_cfm_get_info(amatrix, nrow_global=nao) 217 nmo = SIZE(eigenvalues) 218 ALLOCATE (evals(nao), cevals(nao)) 219 220 ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum 221 CALL cp_cfm_scale(-cone, bmatrix) 222 CALL cp_cfm_heevd(bmatrix, work, evals) 223 evals(:) = -evals(:) 224 nc = nao 225 DO i = 1, nao 226 IF (evals(i) < epseig) THEN 227 nc = i - 1 228 EXIT 229 END IF 230 END DO 231 CPASSERT(nc /= 0) 232 233 IF (nc /= nao) THEN 234 IF (nc < nmo) THEN 235 ! Copy NULL space definition to last vectors of eigenvectors (if needed) 236 ncol = nmo - nc 237 CALL cp_cfm_to_cfm(work, eigenvectors, ncol, nc + 1, nc + 1) 238 END IF 239 ! Set NULL space in eigenvector matrix of S to zero 240 DO icol = nc + 1, nao 241 DO irow = 1, nao 242 CALL cp_cfm_set_element(work, irow, icol, czero) 243 END DO 244 END DO 245 ! Set small eigenvalues to a dummy save value 246 evals(nc + 1:nao) = 1.0_dp 247 END IF 248 ! calculate U*s**(-1/2) 249 cevals(:) = CMPLX(1.0_dp/SQRT(evals(:)), 0.0_dp, KIND=dp) 250 CALL cp_cfm_column_scale(work, cevals) 251 ! Reduce to get U^(-C) * H * U^(-1) 252 CALL cp_cfm_gemm("C", "N", nao, nao, nao, cone, work, amatrix, czero, bmatrix) 253 CALL cp_cfm_gemm("N", "N", nao, nao, nao, cone, bmatrix, work, czero, amatrix) 254 IF (nc /= nao) THEN 255 ! set diagonal values to save large value 256 DO icol = nc + 1, nao 257 CALL cp_cfm_set_element(amatrix, icol, icol, CMPLX(10000.0_dp, 0.0_dp, KIND=dp)) 258 END DO 259 END IF 260 ! Diagonalize 261 CALL cp_cfm_heevd(amatrix, bmatrix, evals) 262 eigenvalues(1:nmo) = evals(1:nmo) 263 nx = MIN(nc, nmo) 264 ! Restore vectors C = U^(-1) * C* 265 CALL cp_cfm_gemm("N", "N", nao, nx, nc, cone, work, bmatrix, czero, eigenvectors) 266 267 DEALLOCATE (evals) 268 269 CALL timestop(handle) 270 271 END SUBROUTINE cp_cfm_geeig_canon 272 273END MODULE cp_cfm_diag 274