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