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