Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief used for collecting diagonalization schemes available for cp_cfm_type
10 : !> \note
11 : !> first version : only one routine right now
12 : !> \author Joost VandeVondele (2003-09)
13 : ! **************************************************************************************************
14 : MODULE cp_cfm_diag
15 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
16 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm, &
17 : cp_cfm_column_scale, &
18 : cp_cfm_scale, &
19 : cp_cfm_triangular_invert, &
20 : cp_cfm_triangular_multiply
21 : USE cp_cfm_types, ONLY: cp_cfm_get_info, &
22 : cp_cfm_set_element, &
23 : cp_cfm_to_cfm, &
24 : cp_cfm_type
25 : #if defined(__DLAF)
26 : USE cp_cfm_dlaf_api, ONLY: cp_cfm_diag_gen_dlaf, &
27 : cp_cfm_diag_dlaf
28 : USE cp_dlaf_utils_api, ONLY: cp_dlaf_initialize, cp_dlaf_create_grid
29 : USE cp_fm_diag, ONLY: diag_type, dlaf_neigvec_min, FM_DIAG_TYPE_DLAF
30 : #endif
31 : USE kinds, ONLY: dp
32 : #if defined (__HAS_IEEE_EXCEPTIONS)
33 : USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
34 : ieee_set_halting_mode, &
35 : IEEE_ALL
36 : #endif
37 :
38 : #include "../base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 : PRIVATE
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_diag'
44 :
45 : PUBLIC :: cp_cfm_heevd, cp_cfm_geeig, cp_cfm_geeig_canon
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief Perform a diagonalisation of a complex matrix
51 : !> \param matrix ...
52 : !> \param eigenvectors ...
53 : !> \param eigenvalues ...
54 : !> \par History
55 : !> 12.2024 Added DLA-Future support [Rocco Meli]
56 : !> \author Joost VandeVondele
57 : ! **************************************************************************************************
58 26256 : SUBROUTINE cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
59 :
60 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
61 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
62 :
63 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd'
64 :
65 : INTEGER :: handle
66 :
67 26256 : CALL timeset(routineN, handle)
68 :
69 : #if defined(__DLAF)
70 : IF (diag_type == FM_DIAG_TYPE_DLAF .AND. matrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN
71 : ! Initialize DLA-Future on-demand; if already initialized, does nothing
72 : CALL cp_dlaf_initialize()
73 :
74 : ! Create DLAF grid from BLACS context; if already present, does nothing
75 : CALL cp_dlaf_create_grid(matrix%matrix_struct%context%get_handle())
76 :
77 : CALL cp_cfm_diag_dlaf(matrix, eigenvectors, eigenvalues)
78 : ELSE
79 : #endif
80 26256 : CALL cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
81 : #if defined(__DLAF)
82 : END IF
83 : #endif
84 :
85 26256 : CALL timestop(handle)
86 :
87 26256 : END SUBROUTINE cp_cfm_heevd
88 :
89 : ! **************************************************************************************************
90 : !> \brief Perform a diagonalisation of a complex matrix
91 : !> \param matrix ...
92 : !> \param eigenvectors ...
93 : !> \param eigenvalues ...
94 : !> \par History
95 : !> - (De)Allocation checks updated (15.02.2011,MK)
96 : !> \author Joost VandeVondele
97 : ! **************************************************************************************************
98 26256 : SUBROUTINE cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
99 :
100 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
101 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
102 :
103 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd_base'
104 :
105 26256 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work
106 : COMPLEX(KIND=dp), DIMENSION(:, :), &
107 26256 : POINTER :: m
108 : INTEGER :: handle, info, liwork, &
109 : lrwork, lwork, n
110 26256 : INTEGER, DIMENSION(:), POINTER :: iwork
111 26256 : REAL(KIND=dp), DIMENSION(:), POINTER :: rwork
112 : #if defined(__parallel)
113 : INTEGER, DIMENSION(9) :: descm, descv
114 : COMPLEX(KIND=dp), DIMENSION(:, :), &
115 26256 : POINTER :: v
116 : #if defined (__HAS_IEEE_EXCEPTIONS)
117 : LOGICAL, DIMENSION(5) :: halt
118 : #endif
119 : #endif
120 :
121 26256 : CALL timeset(routineN, handle)
122 :
123 26256 : n = matrix%matrix_struct%nrow_global
124 26256 : m => matrix%local_data
125 26256 : ALLOCATE (iwork(1), rwork(1), work(1))
126 : ! work space query
127 26256 : lwork = -1
128 26256 : lrwork = -1
129 26256 : liwork = -1
130 :
131 : #if defined(__parallel)
132 26256 : v => eigenvectors%local_data
133 262560 : descm(:) = matrix%matrix_struct%descriptor(:)
134 262560 : descv(:) = eigenvectors%matrix_struct%descriptor(:)
135 : CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
136 26256 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
137 : ! The work space query for lwork does not return always sufficiently large values.
138 : ! Let's add some margin to avoid crashes.
139 26256 : lwork = CEILING(REAL(work(1), KIND=dp)) + 1000
140 : ! needed to correct for a bug in scalapack, unclear how much the right number is
141 26256 : lrwork = CEILING(rwork(1)) + 1000000
142 26256 : liwork = iwork(1)
143 : #else
144 : CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
145 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
146 : lwork = CEILING(REAL(work(1), KIND=dp))
147 : lrwork = CEILING(rwork(1))
148 : liwork = iwork(1)
149 : #endif
150 :
151 26256 : DEALLOCATE (iwork, rwork, work)
152 183792 : ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
153 :
154 : #if defined(__parallel)
155 : ! Scalapack takes advantage of IEEE754 exceptions for speedup.
156 : ! Therefore, we disable floating point traps temporarily.
157 : #if defined (__HAS_IEEE_EXCEPTIONS)
158 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
159 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
160 : #endif
161 :
162 : CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
163 26256 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
164 :
165 : #if defined (__HAS_IEEE_EXCEPTIONS)
166 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
167 : #endif
168 : #else
169 : CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
170 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
171 : eigenvectors%local_data = matrix%local_data
172 : #endif
173 :
174 26256 : DEALLOCATE (iwork, rwork, work)
175 26256 : IF (info /= 0) &
176 0 : CPABORT("Diagonalisation of a complex matrix failed")
177 :
178 26256 : CALL timestop(handle)
179 :
180 26256 : END SUBROUTINE cp_cfm_heevd_base
181 :
182 : ! **************************************************************************************************
183 : !> \brief General Eigenvalue Problem AX = BXE
184 : !> Single option version: Cholesky decomposition of B
185 : !> \param amatrix ...
186 : !> \param bmatrix ...
187 : !> \param eigenvectors ...
188 : !> \param eigenvalues ...
189 : !> \param work ...
190 : !> \par History
191 : !> 12.2024 Added DLA-Future support [Rocco Meli]
192 : ! **************************************************************************************************
193 16710 : SUBROUTINE cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
194 :
195 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
196 : REAL(KIND=dp), DIMENSION(:) :: eigenvalues
197 : TYPE(cp_cfm_type), INTENT(IN) :: work
198 :
199 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig'
200 :
201 : INTEGER :: handle, nao, nmo
202 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
203 :
204 16710 : CALL timeset(routineN, handle)
205 :
206 16710 : CALL cp_cfm_get_info(amatrix, nrow_global=nao)
207 50130 : ALLOCATE (evals(nao))
208 16710 : nmo = SIZE(eigenvalues)
209 :
210 : #if defined(__DLAF)
211 : IF (diag_type == FM_DIAG_TYPE_DLAF .AND. amatrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN
212 : ! Initialize DLA-Future on-demand; if already initialized, does nothing
213 : CALL cp_dlaf_initialize()
214 :
215 : ! Create DLAF grid from BLACS context; if already present, does nothing
216 : CALL cp_dlaf_create_grid(amatrix%matrix_struct%context%get_handle())
217 : CALL cp_dlaf_create_grid(bmatrix%matrix_struct%context%get_handle())
218 : CALL cp_dlaf_create_grid(eigenvectors%matrix_struct%context%get_handle())
219 :
220 : ! Use DLA-Future generalized eigenvalue solver for large matrices
221 : CALL cp_cfm_diag_gen_dlaf(amatrix, bmatrix, work, evals)
222 : ELSE
223 : #endif
224 : ! Cholesky decompose S=U(T)U
225 16710 : CALL cp_cfm_cholesky_decompose(bmatrix)
226 : ! Invert to get U^(-1)
227 16710 : CALL cp_cfm_triangular_invert(bmatrix)
228 : ! Reduce to get U^(-T) * H * U^(-1)
229 16710 : CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
230 16710 : CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
231 : ! Diagonalize
232 16710 : CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
233 : ! Restore vectors C = U^(-1) * C*
234 16710 : CALL cp_cfm_triangular_multiply(bmatrix, work)
235 : #if defined(__DLAF)
236 : END IF
237 : #endif
238 :
239 16710 : CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
240 333414 : eigenvalues(1:nmo) = evals(1:nmo)
241 :
242 16710 : DEALLOCATE (evals)
243 :
244 16710 : CALL timestop(handle)
245 :
246 16710 : END SUBROUTINE cp_cfm_geeig
247 :
248 : ! **************************************************************************************************
249 : !> \brief General Eigenvalue Problem AX = BXE
250 : !> Use canonical orthogonalization
251 : !> \param amatrix ...
252 : !> \param bmatrix ...
253 : !> \param eigenvectors ...
254 : !> \param eigenvalues ...
255 : !> \param work ...
256 : !> \param epseig ...
257 : ! **************************************************************************************************
258 432 : SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
259 :
260 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
261 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
262 : TYPE(cp_cfm_type), INTENT(IN) :: work
263 : REAL(KIND=dp), INTENT(IN) :: epseig
264 :
265 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig_canon'
266 : COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
267 : czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
268 :
269 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
270 : INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
271 : nmo, nx
272 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
273 :
274 432 : CALL timeset(routineN, handle)
275 :
276 : ! Test sizes
277 432 : CALL cp_cfm_get_info(amatrix, nrow_global=nao)
278 432 : nmo = SIZE(eigenvalues)
279 2160 : ALLOCATE (evals(nao), cevals(nao))
280 :
281 : ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
282 432 : CALL cp_cfm_scale(-cone, bmatrix)
283 432 : CALL cp_cfm_heevd(bmatrix, work, evals)
284 19502 : evals(:) = -evals(:)
285 432 : nc = nao
286 18938 : DO i = 1, nao
287 18938 : IF (evals(i) < epseig) THEN
288 72 : nc = i - 1
289 72 : EXIT
290 : END IF
291 : END DO
292 432 : CPASSERT(nc /= 0)
293 :
294 432 : IF (nc /= nao) THEN
295 72 : IF (nc < nmo) THEN
296 : ! Copy NULL space definition to last vectors of eigenvectors (if needed)
297 0 : ncol = nmo - nc
298 0 : CALL cp_cfm_to_cfm(work, eigenvectors, ncol, nc + 1, nc + 1)
299 : END IF
300 : ! Set NULL space in eigenvector matrix of S to zero
301 636 : DO icol = nc + 1, nao
302 80028 : DO irow = 1, nao
303 79956 : CALL cp_cfm_set_element(work, irow, icol, czero)
304 : END DO
305 : END DO
306 : ! Set small eigenvalues to a dummy save value
307 636 : evals(nc + 1:nao) = 1.0_dp
308 : END IF
309 : ! Calculate U*s**(-1/2)
310 19502 : cevals(:) = CMPLX(1.0_dp/SQRT(evals(:)), 0.0_dp, KIND=dp)
311 432 : CALL cp_cfm_column_scale(work, cevals)
312 : ! Reduce to get U^(-C) * H * U^(-1)
313 432 : CALL cp_cfm_gemm("C", "N", nao, nao, nao, cone, work, amatrix, czero, bmatrix)
314 432 : CALL cp_cfm_gemm("N", "N", nao, nao, nao, cone, bmatrix, work, czero, amatrix)
315 432 : IF (nc /= nao) THEN
316 : ! set diagonal values to save large value
317 636 : DO icol = nc + 1, nao
318 636 : CALL cp_cfm_set_element(amatrix, icol, icol, CMPLX(10000.0_dp, 0.0_dp, KIND=dp))
319 : END DO
320 : END IF
321 : ! Diagonalize
322 432 : CALL cp_cfm_heevd(amatrix, bmatrix, evals)
323 6678 : eigenvalues(1:nmo) = evals(1:nmo)
324 432 : nx = MIN(nc, nmo)
325 : ! Restore vectors C = U^(-1) * C*
326 432 : CALL cp_cfm_gemm("N", "N", nao, nx, nc, cone, work, bmatrix, czero, eigenvectors)
327 :
328 432 : DEALLOCATE (evals)
329 :
330 432 : CALL timestop(handle)
331 :
332 864 : END SUBROUTINE cp_cfm_geeig_canon
333 :
334 : END MODULE cp_cfm_diag
|