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