Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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_blacs_env, ONLY: cp_blacs_env_type
16 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
17 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm, &
18 : cp_cfm_column_scale, &
19 : cp_cfm_scale, &
20 : cp_cfm_triangular_invert, &
21 : cp_cfm_triangular_multiply
22 : USE cp_cfm_types, ONLY: cp_cfm_create, &
23 : cp_cfm_get_info, &
24 : cp_cfm_release, &
25 : cp_cfm_set_element, &
26 : cp_cfm_to_cfm, &
27 : cp_cfm_type
28 : USE cp_fm_diag, ONLY: diag_check_requested, &
29 : diag_check_warning_threshold, &
30 : diag_lib_explicit, &
31 : diag_type, &
32 : direct_generalized_diagonalization, &
33 : cusolver_n_min, &
34 : elpa_neigvec_min, &
35 : FM_DIAG_TYPE_CUSOLVER, &
36 : FM_DIAG_TYPE_ELPA, &
37 : FM_DIAG_TYPE_SCALAPACK, &
38 : set_removed_eigval_to
39 : USE cp_cfm_elpa, ONLY: cp_cfm_diag_elpa, &
40 : is_elpa_c_broken
41 : USE cp_fm_cusolver_api, ONLY: cp_cfm_general_cusolver
42 : #if defined(__DLAF)
43 : USE cp_cfm_dlaf_api, ONLY: cp_cfm_diag_gen_dlaf, &
44 : cp_cfm_diag_dlaf
45 : USE cp_dlaf_utils_api, ONLY: cp_dlaf_initialize, cp_dlaf_create_grid
46 : USE cp_fm_diag, ONLY: dlaf_neigvec_min, FM_DIAG_TYPE_DLAF
47 : #endif
48 : USE cp_log_handling, ONLY: cp_to_string
49 : USE kinds, ONLY: default_string_length, &
50 : dp
51 : USE machine, ONLY: default_output_unit
52 : USE mathconstants, ONLY: z_one, &
53 : z_zero
54 : #if defined (__HAS_IEEE_EXCEPTIONS)
55 : USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
56 : ieee_set_halting_mode, &
57 : IEEE_ALL
58 : #endif
59 : #include "../base/base_uses.f90"
60 :
61 : IMPLICIT NONE
62 : PRIVATE
63 :
64 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_diag'
65 :
66 : PUBLIC :: cp_cfm_heevd, cp_cfm_geeig, cp_cfm_geeig_canon, &
67 : cp_cfm_geeig_local, cp_cfm_geeig_canon_local
68 :
69 : ! Caller-owned LAPACK scratch: no module SAVE state and no cross-worker sharing.
70 : TYPE, PUBLIC :: cp_cfm_diag_workspace_type
71 : INTEGER :: n = 0
72 : LOGICAL :: generalized = .FALSE.
73 : COMPLEX(KIND=dp), ALLOCATABLE :: work(:)
74 : REAL(KIND=dp), ALLOCATABLE :: evals(:), rwork(:)
75 : INTEGER, ALLOCATABLE :: iwork(:)
76 : END TYPE cp_cfm_diag_workspace_type
77 : PUBLIC :: cp_cfm_local_workspace_prepare
78 :
79 : CONTAINS
80 :
81 : ! **************************************************************************************************
82 : !> \brief Perform a diagonalisation of a complex matrix
83 : !> \param matrix ...
84 : !> \param eigenvectors ...
85 : !> \param eigenvalues ...
86 : !> \par History
87 : !> 12.2024 Added DLA-Future support [Rocco Meli]
88 : !> 08.2026 Added ELPA support
89 : !> \author Joost VandeVondele
90 : ! **************************************************************************************************
91 64797 : SUBROUTINE cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
92 :
93 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
94 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
95 :
96 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd'
97 :
98 : INTEGER :: handle
99 :
100 64797 : CALL timeset(routineN, handle)
101 :
102 : #if defined(__DLAF)
103 : IF (diag_type == FM_DIAG_TYPE_DLAF .AND. matrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN
104 : ! Initialize DLA-Future on-demand; if already initialized, does nothing
105 : CALL cp_dlaf_initialize()
106 :
107 : ! Create DLAF grid from BLACS context; if already present, does nothing
108 : CALL cp_dlaf_create_grid(matrix%matrix_struct%context%get_handle())
109 :
110 : CALL cp_cfm_diag_dlaf(matrix, eigenvectors, eigenvalues)
111 : ELSE
112 : #endif
113 : ! We don't trust ELPA with very small matrices and use it for complex matrices
114 : ! only when the diagonalization library was requested explicitly.
115 : ! A runtime correctness check may have disabled ELPA for mis-compiled BLOCK2 kernels.
116 : IF (diag_type == FM_DIAG_TYPE_ELPA .AND. diag_lib_explicit .AND. &
117 64797 : .NOT. is_elpa_c_broken() .AND. &
118 : matrix%matrix_struct%nrow_global >= elpa_neigvec_min) THEN
119 72 : CALL cp_cfm_diag_elpa(matrix, eigenvectors, eigenvalues)
120 : ELSE
121 64725 : CALL cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
122 : END IF
123 : #if defined(__DLAF)
124 : END IF
125 : #endif
126 :
127 64797 : CALL timestop(handle)
128 :
129 64797 : END SUBROUTINE cp_cfm_heevd
130 :
131 : ! **************************************************************************************************
132 : !> \brief Perform a diagonalisation of a complex matrix
133 : !> \param matrix ...
134 : !> \param eigenvectors ...
135 : !> \param eigenvalues ...
136 : !> \par History
137 : !> - (De)Allocation checks updated (15.02.2011,MK)
138 : !> \author Joost VandeVondele
139 : ! **************************************************************************************************
140 64725 : SUBROUTINE cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
141 :
142 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
143 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
144 :
145 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd_base'
146 :
147 64725 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work
148 : COMPLEX(KIND=dp), DIMENSION(:, :), &
149 64725 : POINTER :: m
150 : INTEGER :: handle, info, liwork, &
151 : lrwork, lwork, n
152 64725 : INTEGER, DIMENSION(:), POINTER :: iwork
153 64725 : REAL(KIND=dp), DIMENSION(:), POINTER :: rwork
154 : #if defined(__parallel)
155 : INTEGER, DIMENSION(9) :: descm, descv
156 : COMPLEX(KIND=dp), DIMENSION(:, :), &
157 64725 : POINTER :: v
158 : #endif
159 : #if defined (__HAS_IEEE_EXCEPTIONS)
160 : LOGICAL, DIMENSION(5) :: halt
161 : #endif
162 :
163 64725 : CALL timeset(routineN, handle)
164 :
165 64725 : n = matrix%matrix_struct%nrow_global
166 64725 : m => matrix%local_data
167 64725 : ALLOCATE (iwork(1), rwork(1), work(1))
168 : ! work space query
169 64725 : lwork = -1
170 64725 : lrwork = -1
171 64725 : liwork = -1
172 :
173 : #if defined(__parallel)
174 64725 : v => eigenvectors%local_data
175 647250 : descm(:) = matrix%matrix_struct%descriptor(:)
176 647250 : descv(:) = eigenvectors%matrix_struct%descriptor(:)
177 : CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
178 64725 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
179 : ! The work space query for lwork does not return always sufficiently large values.
180 : ! Let's add some margin to avoid crashes.
181 64725 : lwork = CEILING(REAL(work(1), KIND=dp)) + 1000
182 : ! needed to correct for a bug in scalapack, unclear how much the right number is
183 64725 : lrwork = CEILING(rwork(1)) + 1000000
184 64725 : liwork = iwork(1)
185 : #else
186 : CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
187 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
188 : lwork = CEILING(REAL(work(1), KIND=dp))
189 : lrwork = CEILING(rwork(1))
190 : liwork = iwork(1)
191 : #endif
192 :
193 64725 : DEALLOCATE (iwork, rwork, work)
194 453075 : ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
195 :
196 : ! (Sca-)LAPACK takes advantage of IEEE754 exceptions for speedup.
197 : ! Therefore, we disable floating point traps temporarily.
198 : #if defined (__HAS_IEEE_EXCEPTIONS)
199 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
200 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
201 : #endif
202 : #if defined(__parallel)
203 : CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
204 64725 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
205 : #else
206 : CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
207 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
208 : eigenvectors%local_data = matrix%local_data
209 : #endif
210 : #if defined (__HAS_IEEE_EXCEPTIONS)
211 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
212 : #endif
213 :
214 64725 : DEALLOCATE (iwork, rwork, work)
215 64725 : IF (info /= 0) CPABORT("Diagonalisation of a complex matrix failed")
216 :
217 64725 : CALL timestop(handle)
218 :
219 64725 : END SUBROUTINE cp_cfm_heevd_base
220 :
221 : ! **************************************************************************************************
222 : !> \brief Check C^H*S*C = I for a generalized complex eigenvalue problem.
223 : !> \param overlap original overlap matrix S; used as work matrix and overwritten
224 : !> \param eigenvectors eigenvectors C to be checked
225 : !> \param scratch work matrix
226 : !> \param nvec ...
227 : ! **************************************************************************************************
228 16 : SUBROUTINE check_generalized_diag(overlap, eigenvectors, scratch, nvec)
229 :
230 : TYPE(cp_cfm_type), INTENT(IN) :: eigenvectors
231 : TYPE(cp_cfm_type), INTENT(INOUT) :: overlap, scratch
232 : INTEGER, INTENT(IN) :: nvec
233 :
234 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_generalized_diag'
235 :
236 : CHARACTER(LEN=default_string_length) :: diag_type_name
237 : COMPLEX(KIND=dp) :: gold, test
238 : INTEGER :: handle, i, j, ncol, nrow, output_unit
239 : REAL(KIND=dp) :: eps, eps_abort, eps_warning
240 : #if defined(__parallel)
241 : TYPE(cp_blacs_env_type), POINTER :: context
242 : INTEGER :: il, jl, ipcol, iprow, &
243 : mypcol, myprow, npcol, nprow
244 : INTEGER, DIMENSION(9) :: desca
245 : #endif
246 :
247 16 : CALL timeset(routineN, handle)
248 :
249 16 : IF (.NOT. diag_check_requested()) THEN
250 0 : CALL timestop(handle)
251 0 : RETURN
252 : END IF
253 :
254 16 : output_unit = default_output_unit
255 16 : eps_warning = diag_check_warning_threshold()
256 16 : eps_abort = 10.0_dp*eps_warning
257 :
258 16 : nrow = eigenvectors%matrix_struct%nrow_global
259 16 : ncol = MIN(eigenvectors%matrix_struct%ncol_global, nvec)
260 :
261 16 : CALL cp_cfm_gemm("N", "N", nrow, ncol, nrow, z_one, overlap, eigenvectors, z_zero, scratch)
262 16 : CALL cp_cfm_gemm("C", "N", ncol, ncol, nrow, z_one, eigenvectors, scratch, z_zero, overlap)
263 :
264 16 : gold = z_zero
265 16 : test = z_zero
266 16 : eps = 0.0_dp
267 :
268 : #if defined(__parallel)
269 16 : context => overlap%matrix_struct%context
270 16 : myprow = context%mepos(1)
271 16 : mypcol = context%mepos(2)
272 16 : nprow = context%num_pe(1)
273 16 : npcol = context%num_pe(2)
274 160 : desca(:) = overlap%matrix_struct%descriptor(:)
275 160 : outer: DO j = 1, ncol
276 1456 : DO i = 1, ncol
277 1296 : CALL infog2l(i, j, desca, nprow, npcol, myprow, mypcol, il, jl, iprow, ipcol)
278 1440 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
279 648 : gold = MERGE(z_zero, z_one, i /= j)
280 648 : test = overlap%local_data(il, jl)
281 648 : eps = ABS(test - gold)
282 648 : IF (eps > eps_warning) EXIT outer
283 : END IF
284 : END DO
285 : END DO outer
286 : #else
287 : outer: DO j = 1, ncol
288 : DO i = 1, ncol
289 : gold = MERGE(z_zero, z_one, i /= j)
290 : test = overlap%local_data(i, j)
291 : eps = ABS(test - gold)
292 : IF (eps > eps_warning) EXIT outer
293 : END DO
294 : END DO outer
295 : #endif
296 :
297 16 : IF (eps > eps_warning) THEN
298 0 : IF (diag_type == FM_DIAG_TYPE_SCALAPACK) THEN
299 0 : diag_type_name = "HEGVX"
300 0 : ELSE IF (diag_type == FM_DIAG_TYPE_CUSOLVER) THEN
301 0 : diag_type_name = "CUSOLVER"
302 0 : ELSE IF (diag_type == FM_DIAG_TYPE_ELPA .AND. diag_lib_explicit) THEN
303 0 : diag_type_name = "ELPA"
304 : #if defined(__DLAF)
305 : ELSE IF (diag_type == FM_DIAG_TYPE_DLAF) THEN
306 : diag_type_name = "DLAF"
307 : #endif
308 : ELSE
309 0 : diag_type_name = "generalized eigensolver"
310 : END IF
311 : WRITE (UNIT=output_unit, FMT="(/,T2,A,/,T2,A,I0,A,I0,A,ES10.3,/,T2,A,F0.0,A,ES10.3)") &
312 0 : "The generalized eigenvectors returned by "//TRIM(diag_type_name)//" are not S-orthonormal", &
313 0 : "Absolute deviation of matrix element (", i, ", ", j, ") is ", eps, &
314 0 : "The deviation from the expected value ", REAL(gold, KIND=dp), " is", eps
315 0 : IF (eps > eps_abort) THEN
316 : CALL cp_abort(__LOCATION__, &
317 0 : "ERROR in "//routineN//": Check of generalized matrix diagonalization failed")
318 : ELSE
319 0 : CPWARN("Check of generalized matrix diagonalization failed in routine "//routineN)
320 : END IF
321 : END IF
322 :
323 16 : CALL timestop(handle)
324 :
325 : END SUBROUTINE check_generalized_diag
326 :
327 : ! **************************************************************************************************
328 : !> \brief General Eigenvalue Problem AX = BXE
329 : !> Single option version: Cholesky decomposition of B
330 : !> \param amatrix ...
331 : !> \param bmatrix ...
332 : !> \param eigenvectors ...
333 : !> \param eigenvalues ...
334 : !> \param work ...
335 : !> \param lowest_subset compute only the requested lowest eigenpairs with ScaLAPACK when available
336 : !> \par History
337 : !> 12.2024 Added DLA-Future support [Rocco Meli]
338 : ! **************************************************************************************************
339 36870 : SUBROUTINE cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work, lowest_subset)
340 :
341 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
342 : REAL(KIND=dp), DIMENSION(:) :: eigenvalues
343 : TYPE(cp_cfm_type), INTENT(IN) :: work
344 : LOGICAL, INTENT(IN), OPTIONAL :: lowest_subset
345 :
346 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig'
347 :
348 : INTEGER :: handle, nao, nmo
349 : LOGICAL :: check_eigenvectors, use_lowest_subset
350 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
351 : TYPE(cp_cfm_type) :: overlap_check, scratch_check
352 :
353 36870 : CALL timeset(routineN, handle)
354 :
355 36870 : CALL cp_cfm_get_info(amatrix, nrow_global=nao)
356 110610 : ALLOCATE (evals(nao))
357 36870 : nmo = SIZE(eigenvalues)
358 36870 : check_eigenvectors = diag_check_requested()
359 36870 : use_lowest_subset = .FALSE.
360 36870 : IF (PRESENT(lowest_subset)) use_lowest_subset = lowest_subset .AND. nmo < nao
361 : #if !defined(__parallel)
362 : use_lowest_subset = .FALSE.
363 : #endif
364 :
365 : IF (use_lowest_subset) THEN
366 : #if defined(__parallel)
367 176 : IF (check_eigenvectors) THEN
368 0 : CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
369 0 : CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
370 0 : CALL cp_cfm_to_cfm(bmatrix, overlap_check)
371 : END IF
372 176 : CALL cp_cfm_geeig_scalapack(amatrix, bmatrix, work, evals(1:nmo))
373 176 : IF (check_eigenvectors) THEN
374 0 : CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
375 0 : CALL cp_cfm_release(scratch_check)
376 0 : CALL cp_cfm_release(overlap_check)
377 : END IF
378 : #endif
379 36694 : ELSE IF (diag_type == FM_DIAG_TYPE_CUSOLVER .AND. direct_generalized_diagonalization .AND. &
380 : nao >= cusolver_n_min) THEN
381 : ! Use cuSolverMP generalized eigenvalue solver without a CP2K-side
382 : ! Cholesky reduction.
383 0 : IF (check_eigenvectors) THEN
384 0 : CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
385 0 : CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
386 0 : CALL cp_cfm_to_cfm(bmatrix, overlap_check)
387 : END IF
388 0 : CALL cp_cfm_general_cusolver(amatrix, bmatrix, work, evals)
389 0 : IF (check_eigenvectors) THEN
390 0 : CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
391 0 : CALL cp_cfm_release(scratch_check)
392 0 : CALL cp_cfm_release(overlap_check)
393 : END IF
394 : #if defined(__DLAF)
395 : ELSE IF (diag_type == FM_DIAG_TYPE_DLAF .AND. direct_generalized_diagonalization .AND. &
396 : nao >= dlaf_neigvec_min) THEN
397 : ! Initialize DLA-Future on-demand; if already initialized, does nothing
398 : CALL cp_dlaf_initialize()
399 :
400 : ! Create DLAF grid from BLACS context; if already present, does nothing
401 : CALL cp_dlaf_create_grid(amatrix%matrix_struct%context%get_handle())
402 : CALL cp_dlaf_create_grid(bmatrix%matrix_struct%context%get_handle())
403 : CALL cp_dlaf_create_grid(eigenvectors%matrix_struct%context%get_handle())
404 :
405 : ! Use DLA-Future generalized eigenvalue solver for large matrices
406 : IF (check_eigenvectors) THEN
407 : CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
408 : CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
409 : CALL cp_cfm_to_cfm(bmatrix, overlap_check)
410 : END IF
411 : CALL cp_cfm_diag_gen_dlaf(amatrix, bmatrix, work, evals)
412 : IF (check_eigenvectors) THEN
413 : CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
414 : CALL cp_cfm_release(scratch_check)
415 : CALL cp_cfm_release(overlap_check)
416 : END IF
417 : #endif
418 : #if defined(__parallel)
419 36694 : ELSE IF (diag_type == FM_DIAG_TYPE_SCALAPACK .AND. direct_generalized_diagonalization) THEN
420 : ! Use ScaLAPACK generalized eigenvalue solver without a CP2K-side
421 : ! Cholesky reduction.
422 16 : IF (check_eigenvectors) THEN
423 16 : CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
424 16 : CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
425 16 : CALL cp_cfm_to_cfm(bmatrix, overlap_check)
426 : END IF
427 16 : CALL cp_cfm_geeig_scalapack(amatrix, bmatrix, work, evals)
428 16 : IF (check_eigenvectors) THEN
429 16 : CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
430 16 : CALL cp_cfm_release(scratch_check)
431 16 : CALL cp_cfm_release(overlap_check)
432 : END IF
433 : #endif
434 : ELSE
435 : ! Cholesky decompose S=U(T)U
436 36678 : CALL cp_cfm_cholesky_decompose(bmatrix)
437 : ! Invert to get U^(-1)
438 36678 : CALL cp_cfm_triangular_invert(bmatrix)
439 : ! Reduce to get U^(-T) * H * U^(-1)
440 36678 : CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
441 36678 : CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
442 : ! Diagonalize
443 36678 : CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
444 : ! Restore vectors C = U^(-1) * C*
445 36678 : CALL cp_cfm_triangular_multiply(bmatrix, work)
446 : END IF
447 :
448 36870 : CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
449 796076 : eigenvalues(1:nmo) = evals(1:nmo)
450 :
451 36870 : DEALLOCATE (evals)
452 :
453 36870 : CALL timestop(handle)
454 :
455 36870 : END SUBROUTINE cp_cfm_geeig
456 :
457 : ! **************************************************************************************************
458 : !> \brief General Eigenvalue Problem AX = BXE using ScaLAPACK PZHEGVX.
459 : !> \param amatrix ...
460 : !> \param bmatrix ...
461 : !> \param eigenvectors ...
462 : !> \param eigenvalues ...
463 : ! **************************************************************************************************
464 192 : SUBROUTINE cp_cfm_geeig_scalapack(amatrix, bmatrix, eigenvectors, eigenvalues)
465 :
466 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
467 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
468 :
469 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig_scalapack'
470 :
471 : #if defined(__parallel)
472 : REAL(KIND=dp), PARAMETER :: orfac = -1.0_dp, &
473 : vl = 0.0_dp, &
474 : vu = 0.0_dp
475 :
476 192 : COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
477 192 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a, b, z
478 : INTEGER :: handle, info, liwork, lwork, lrwork, &
479 : m, n, nb, neig, npcol, nprow, nz
480 : INTEGER, DIMENSION(9) :: desca, descb, descz
481 192 : INTEGER, DIMENSION(:), ALLOCATABLE :: iclustr, ifail, iwork
482 : REAL(KIND=dp) :: abstol
483 192 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: gap, rwork, w
484 :
485 : INTEGER :: mq0, nn, np0, npe
486 : INTEGER, EXTERNAL :: iceil, numroc
487 : REAL(KIND=dp), EXTERNAL :: dlamch
488 : #if defined (__HAS_IEEE_EXCEPTIONS)
489 : LOGICAL, DIMENSION(5) :: halt
490 : #endif
491 : #else
492 : INTEGER :: handle
493 : #endif
494 :
495 192 : CALL timeset(routineN, handle)
496 :
497 : #if defined(__parallel)
498 192 : n = amatrix%matrix_struct%nrow_global
499 192 : neig = MIN(SIZE(eigenvalues), n)
500 :
501 192 : IF (neig == 0) THEN
502 0 : CALL timestop(handle)
503 0 : RETURN
504 : END IF
505 :
506 192 : IF (amatrix%matrix_struct%nrow_block /= amatrix%matrix_struct%ncol_block) THEN
507 0 : CPABORT("ERROR in "//routineN//": Invalid blocksize (no square blocks) found")
508 : END IF
509 :
510 192 : a => amatrix%local_data
511 192 : b => bmatrix%local_data
512 192 : z => eigenvectors%local_data
513 1920 : desca(:) = amatrix%matrix_struct%descriptor(:)
514 1920 : descb(:) = bmatrix%matrix_struct%descriptor(:)
515 1920 : descz(:) = eigenvectors%matrix_struct%descriptor(:)
516 :
517 192 : nprow = amatrix%matrix_struct%context%num_pe(1)
518 192 : npcol = amatrix%matrix_struct%context%num_pe(2)
519 192 : npe = nprow*npcol
520 192 : nb = amatrix%matrix_struct%nrow_block
521 192 : nn = MAX(n, nb, 2)
522 192 : np0 = numroc(nn, nb, 0, 0, nprow)
523 192 : mq0 = MAX(numroc(nn, nb, 0, 0, npcol), nb)
524 :
525 192 : lwork = n + (np0 + mq0 + nb)*nb
526 192 : lrwork = 4*n + MAX(5*nn, np0*mq0) + iceil(neig, npe)*nn + MAX(0, neig - 1)*n
527 192 : liwork = 6*MAX(n, npe + 1, 4)
528 :
529 576 : ALLOCATE (gap(npe))
530 192 : gap = 0.0_dp
531 576 : ALLOCATE (iclustr(2*npe))
532 192 : iclustr = 0
533 576 : ALLOCATE (ifail(n))
534 192 : ifail = 0
535 576 : ALLOCATE (iwork(liwork))
536 576 : ALLOCATE (rwork(lrwork))
537 576 : ALLOCATE (w(n))
538 576 : ALLOCATE (work(lwork))
539 :
540 192 : abstol = 2.0_dp*dlamch("S")
541 :
542 : #if defined (__HAS_IEEE_EXCEPTIONS)
543 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
544 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
545 : #endif
546 : CALL pzhegvx(1, "V", "I", "U", n, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb, &
547 : vl, vu, 1, neig, abstol, m, nz, w(1), orfac, z(1, 1), 1, 1, descz, &
548 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, ifail(1), &
549 192 : iclustr(1), gap(1), info)
550 : #if defined (__HAS_IEEE_EXCEPTIONS)
551 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
552 : #endif
553 :
554 192 : IF (info /= 0 .OR. m < neig .OR. nz < neig) THEN
555 0 : CPABORT("ERROR in PZHEGVX (ScaLAPACK), info="//TRIM(cp_to_string(info)))
556 : END IF
557 :
558 1209 : eigenvalues(:) = 0.0_dp
559 1209 : eigenvalues(1:neig) = w(1:neig)
560 :
561 192 : DEALLOCATE (gap, iclustr, ifail, iwork, rwork, w, work)
562 : #else
563 : MARK_USED(amatrix)
564 : MARK_USED(bmatrix)
565 : MARK_USED(eigenvectors)
566 : MARK_USED(eigenvalues)
567 : CPABORT("ERROR in "//routineN//": PZHEGVX requested without ScaLAPACK support")
568 : #endif
569 :
570 192 : CALL timestop(handle)
571 :
572 192 : END SUBROUTINE cp_cfm_geeig_scalapack
573 :
574 : ! **************************************************************************************************
575 : !> \brief General Eigenvalue Problem AX = BXE
576 : !> Use canonical orthogonalization
577 : !> \param amatrix ...
578 : !> \param bmatrix ...
579 : !> \param eigenvectors ...
580 : !> \param eigenvalues ...
581 : !> \param work ...
582 : !> \param epseig ...
583 : !> \param nmo_retained ...
584 : ! **************************************************************************************************
585 2596 : SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig, &
586 : nmo_retained)
587 :
588 : TYPE(cp_cfm_type), INTENT(INOUT) :: amatrix, bmatrix, eigenvectors
589 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
590 : TYPE(cp_cfm_type), INTENT(INOUT) :: work
591 : REAL(KIND=dp), INTENT(IN) :: epseig
592 : INTEGER, INTENT(OUT), OPTIONAL :: nmo_retained
593 :
594 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig_canon'
595 :
596 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
597 : INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
598 : nmo, nx
599 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
600 :
601 2596 : CALL timeset(routineN, handle)
602 :
603 : ! Test sizes
604 2596 : CALL cp_cfm_get_info(amatrix, nrow_global=nao)
605 2596 : nmo = SIZE(eigenvalues)
606 12980 : ALLOCATE (evals(nao), cevals(nao))
607 :
608 : ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
609 2596 : CALL cp_cfm_scale(-z_one, bmatrix)
610 2596 : CALL cp_cfm_heevd(bmatrix, work, evals)
611 81048 : evals(:) = -evals(:)
612 2596 : nc = nao
613 80760 : DO i = 1, nao
614 80760 : IF (evals(i) < epseig) THEN
615 16 : nc = i - 1
616 16 : EXIT
617 : END IF
618 : END DO
619 2596 : CPASSERT(nc /= 0)
620 :
621 2596 : IF (nc /= nao) THEN
622 16 : IF (nc < nmo) THEN
623 : ! Copy NULL space definition to last vectors of eigenvectors (if needed)
624 0 : ncol = nmo - nc
625 0 : CALL cp_cfm_to_cfm(work, eigenvectors, ncol, nc + 1, nc + 1)
626 : END IF
627 : ! Set NULL space in eigenvector matrix of S to zero
628 304 : DO icol = nc + 1, nao
629 50992 : DO irow = 1, nao
630 50976 : CALL cp_cfm_set_element(work, irow, icol, z_zero)
631 : END DO
632 : END DO
633 : ! Set small eigenvalues to a dummy save value
634 304 : evals(nc + 1:nao) = 1.0_dp
635 : END IF
636 : ! Calculate U*s**(-1/2)
637 81048 : cevals(:) = CMPLX(1.0_dp/SQRT(evals(:)), 0.0_dp, KIND=dp)
638 2596 : CALL cp_cfm_column_scale(work, cevals)
639 : ! Reduce to get U^(-C) * H * U^(-1)
640 2596 : CALL cp_cfm_gemm("C", "N", nao, nao, nao, z_one, work, amatrix, z_zero, bmatrix)
641 2596 : CALL cp_cfm_gemm("N", "N", nao, nao, nao, z_one, bmatrix, work, z_zero, amatrix)
642 2596 : IF (nc /= nao) THEN
643 : ! set diagonal values to save large value
644 304 : DO icol = nc + 1, nao
645 : CALL cp_cfm_set_element(amatrix, icol, icol, &
646 304 : CMPLX(set_removed_eigval_to, 0.0_dp, KIND=dp))
647 : END DO
648 : END IF
649 : ! Diagonalize
650 2596 : CALL cp_cfm_heevd(amatrix, bmatrix, evals)
651 31910 : eigenvalues(1:nmo) = evals(1:nmo)
652 2596 : nx = MIN(nc, nmo)
653 : ! Restore vectors C = U^(-1) * C*
654 2596 : CALL cp_cfm_gemm("N", "N", nao, nx, nc, z_one, work, bmatrix, z_zero, eigenvectors)
655 :
656 : ! Number of basis modes that survived the linear-dependency filter. The remaining
657 : ! nao - nc entries of eigenvalues(:) are the placeholders set above.
658 2596 : IF (PRESENT(nmo_retained)) nmo_retained = nc
659 :
660 2596 : DEALLOCATE (evals)
661 :
662 2596 : CALL timestop(handle)
663 :
664 5192 : END SUBROUTINE cp_cfm_geeig_canon
665 :
666 : ! **************************************************************************************************
667 : !> \brief Solve a generalized complex eigenproblem using the local LAPACK backend.
668 : !> This routine is restricted to a one-rank BLACS grid. It deliberately
669 : !> avoids ScaLAPACK so independent k-points can be evaluated concurrently
670 : !> without making overlapping MPI calls from OpenMP worker threads.
671 : !> \param amatrix Hamiltonian, overwritten
672 : !> \param bmatrix overlap matrix, overwritten
673 : !> \param eigenvectors eigenvectors
674 : !> \param eigenvalues eigenvalues
675 : !> \param ws reusable private LAPACK scratch
676 : ! **************************************************************************************************
677 42325 : SUBROUTINE cp_cfm_geeig_local(amatrix, bmatrix, eigenvectors, eigenvalues, ws)
678 :
679 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
680 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
681 : TYPE(cp_cfm_diag_workspace_type), INTENT(INOUT) :: ws
682 :
683 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig_local'
684 :
685 42325 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a, b, z
686 : INTEGER :: handle, info, n, nmo
687 : #if defined (__HAS_IEEE_EXCEPTIONS)
688 : LOGICAL, DIMENSION(5) :: halt
689 : #endif
690 :
691 42325 : CALL timeset(routineN, handle)
692 :
693 126975 : CPASSERT(PRODUCT(amatrix%matrix_struct%context%num_pe) == 1)
694 42325 : n = amatrix%matrix_struct%nrow_global
695 42325 : nmo = MIN(n, SIZE(eigenvalues))
696 42325 : a => amatrix%local_data
697 42325 : b => bmatrix%local_data
698 42325 : z => eigenvectors%local_data
699 :
700 42325 : CALL cp_cfm_local_workspace_prepare(ws, n, a, b)
701 :
702 : #if defined (__HAS_IEEE_EXCEPTIONS)
703 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
704 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
705 : #endif
706 : CALL zhegvd(1, 'V', 'U', n, a(1, 1), SIZE(a, 1), b(1, 1), SIZE(b, 1), ws%evals(1), &
707 42325 : ws%work(1), SIZE(ws%work), ws%rwork(1), SIZE(ws%rwork), ws%iwork(1), SIZE(ws%iwork), info)
708 : #if defined (__HAS_IEEE_EXCEPTIONS)
709 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
710 : #endif
711 42325 : IF (info /= 0) THEN
712 0 : CALL cp_abort(__LOCATION__, "Local ZHEGVD failed, info="//TRIM(cp_to_string(info)))
713 : END IF
714 :
715 875424 : eigenvalues(1:nmo) = ws%evals(1:nmo)
716 74654112 : z(1:n, 1:nmo) = a(1:n, 1:nmo)
717 :
718 42325 : CALL timestop(handle)
719 :
720 42325 : END SUBROUTINE cp_cfm_geeig_local
721 :
722 : ! **************************************************************************************************
723 : !> \brief Canonical generalized complex diagonalization on a one-rank BLACS grid.
724 : !> \param amatrix Hamiltonian, overwritten
725 : !> \param bmatrix overlap matrix, overwritten and used as work storage
726 : !> \param eigenvectors eigenvectors
727 : !> \param eigenvalues eigenvalues
728 : !> \param work work matrix
729 : !> \param epseig overlap eigenvalue threshold
730 : !> \param ws reusable private LAPACK scratch
731 : ! **************************************************************************************************
732 1400 : SUBROUTINE cp_cfm_geeig_canon_local(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig, ws)
733 :
734 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
735 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
736 : TYPE(cp_cfm_type), INTENT(IN) :: work
737 : REAL(KIND=dp), INTENT(IN) :: epseig
738 : TYPE(cp_cfm_diag_workspace_type), INTENT(INOUT) :: ws
739 :
740 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig_canon_local'
741 :
742 1400 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a, b, u, z
743 : INTEGER :: handle, i, info, n, nc, nmo, nx
744 :
745 1400 : CALL timeset(routineN, handle)
746 :
747 4200 : CPASSERT(PRODUCT(amatrix%matrix_struct%context%num_pe) == 1)
748 1400 : n = amatrix%matrix_struct%nrow_global
749 1400 : nmo = MIN(n, SIZE(eigenvalues))
750 1400 : a => amatrix%local_data
751 1400 : b => bmatrix%local_data
752 1400 : u => work%local_data
753 1400 : z => eigenvectors%local_data
754 1400 : CALL cp_cfm_local_workspace_prepare(ws, n, b)
755 :
756 1024034 : b(1:n, 1:n) = -b(1:n, 1:n)
757 1400 : CALL cp_cfm_heevd_local(b, u, info, ws)
758 1400 : IF (info /= 0) THEN
759 0 : CALL cp_abort(__LOCATION__, "Local overlap ZHEEVD failed, info="//TRIM(cp_to_string(info)))
760 : END IF
761 29494 : ws%evals(:) = -ws%evals(:)
762 1400 : nc = n
763 29356 : DO i = 1, n
764 29356 : IF (ws%evals(i) < epseig) THEN
765 28 : nc = i - 1
766 28 : EXIT
767 : END IF
768 : END DO
769 1400 : CPASSERT(nc /= 0)
770 :
771 1400 : IF (nc < n) THEN
772 28 : IF (nc < nmo) z(1:n, nc + 1:nmo) = u(1:n, nc + 1:nmo)
773 14518 : u(1:n, nc + 1:n) = z_zero
774 166 : ws%evals(nc + 1:n) = 1.0_dp
775 : END IF
776 29494 : DO i = 1, n
777 1024034 : u(1:n, i) = u(1:n, i)/SQRT(ws%evals(i))
778 : END DO
779 :
780 : CALL zgemm('C', 'N', n, n, n, z_one, u(1, 1), SIZE(u, 1), a(1, 1), SIZE(a, 1), &
781 1400 : z_zero, b(1, 1), SIZE(b, 1))
782 : CALL zgemm('N', 'N', n, n, n, z_one, b(1, 1), SIZE(b, 1), u(1, 1), SIZE(u, 1), &
783 1400 : z_zero, a(1, 1), SIZE(a, 1))
784 1400 : IF (nc < n) THEN
785 166 : DO i = nc + 1, n
786 166 : a(i, i) = CMPLX(10000.0_dp, 0.0_dp, KIND=dp)
787 : END DO
788 : END IF
789 :
790 1400 : CALL cp_cfm_heevd_local(a, b, info, ws)
791 1400 : IF (info /= 0) THEN
792 0 : CALL cp_abort(__LOCATION__, "Local Hamiltonian ZHEEVD failed, info="//TRIM(cp_to_string(info)))
793 : END IF
794 11135 : eigenvalues(1:nmo) = ws%evals(1:nmo)
795 1400 : nx = MIN(nc, nmo)
796 : CALL zgemm('N', 'N', n, nx, nc, z_one, u(1, 1), SIZE(u, 1), b(1, 1), SIZE(b, 1), &
797 1400 : z_zero, z(1, 1), SIZE(z, 1))
798 :
799 1400 : CALL timestop(handle)
800 :
801 1400 : END SUBROUTINE cp_cfm_geeig_canon_local
802 :
803 : ! **************************************************************************************************
804 : !> \brief Local LAPACK ZHEEVD helper. The eigenvectors are copied to vectors.
805 : !> \param matrix ...
806 : !> \param vectors ...
807 : !> \param info ...
808 : !> \param ws prepared private LAPACK scratch
809 : ! **************************************************************************************************
810 2800 : SUBROUTINE cp_cfm_heevd_local(matrix, vectors, info, ws)
811 :
812 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix, vectors
813 : INTEGER, INTENT(OUT) :: info
814 :
815 : TYPE(cp_cfm_diag_workspace_type), INTENT(INOUT) :: ws
816 : INTEGER :: n
817 : #if defined (__HAS_IEEE_EXCEPTIONS)
818 : LOGICAL, DIMENSION(5) :: halt
819 : #endif
820 :
821 2800 : n = ws%n
822 : #if defined (__HAS_IEEE_EXCEPTIONS)
823 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
824 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
825 : #endif
826 : CALL zheevd('V', 'U', n, matrix(1, 1), SIZE(matrix, 1), ws%evals(1), &
827 2800 : ws%work(1), SIZE(ws%work), ws%rwork(1), SIZE(ws%rwork), ws%iwork(1), SIZE(ws%iwork), info)
828 : #if defined (__HAS_IEEE_EXCEPTIONS)
829 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
830 : #endif
831 2048068 : vectors(1:n, 1:n) = matrix(1:n, 1:n)
832 :
833 2800 : END SUBROUTINE cp_cfm_heevd_local
834 :
835 : ! **************************************************************************************************
836 : !> \brief Query LAPACK once per matrix order/solver and retain private scratch across k points.
837 : !> \param ws caller-owned workspace
838 : !> \param n matrix order
839 : !> \param a query matrix; workspace queries do not read its values
840 : !> \param b optional overlap selects ZHEGVD; otherwise query ZHEEVD
841 : ! **************************************************************************************************
842 67151 : SUBROUTINE cp_cfm_local_workspace_prepare(ws, n, a, b)
843 : TYPE(cp_cfm_diag_workspace_type), INTENT(INOUT) :: ws
844 : INTEGER, INTENT(IN) :: n
845 : COMPLEX(KIND=dp), INTENT(INOUT) :: a(:, :)
846 : COMPLEX(KIND=dp), INTENT(INOUT), OPTIONAL :: b(:, :)
847 :
848 : COMPLEX(KIND=dp) :: work(1)
849 : INTEGER :: info, iwork(1)
850 : REAL(KIND=dp) :: rwork(1)
851 :
852 67151 : IF (ws%n == n .AND. (ws%generalized .EQV. PRESENT(b))) RETURN
853 23426 : IF (ALLOCATED(ws%work)) DEALLOCATE (ws%work, ws%rwork, ws%iwork, ws%evals)
854 70278 : ALLOCATE (ws%evals(n))
855 23426 : IF (PRESENT(b)) THEN
856 : CALL zhegvd(1, 'V', 'U', n, a(1, 1), SIZE(a, 1), b(1, 1), SIZE(b, 1), ws%evals(1), &
857 22830 : work(1), -1, rwork(1), -1, iwork(1), -1, info)
858 : ELSE
859 596 : CALL zheevd('V', 'U', n, a(1, 1), SIZE(a, 1), ws%evals(1), work(1), -1, rwork(1), -1, iwork(1), -1, info)
860 : END IF
861 23426 : IF (info /= 0) THEN
862 : CALL cp_abort(__LOCATION__, &
863 0 : "Local LAPACK workspace query failed, info="//TRIM(cp_to_string(info)))
864 : END IF
865 0 : ALLOCATE (ws%work(MAX(1, CEILING(REAL(work(1), dp)))), &
866 163982 : ws%rwork(MAX(1, CEILING(rwork(1)))), ws%iwork(MAX(1, iwork(1))))
867 23426 : ws%n = n
868 23426 : ws%generalized = PRESENT(b)
869 : END SUBROUTINE cp_cfm_local_workspace_prepare
870 :
871 0 : END MODULE cp_cfm_diag
|
