(git:f824184)
Loading...
Searching...
No Matches
cp_cfm_diag.F
Go to the documentation of this file.
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! **************************************************************************************************
14MODULE 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_heevd_local, cp_cfm_local_workspace_prepare
80
81CONTAINS
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 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 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 .NOT. is_elpa_c_broken() .AND. &
120 matrix%matrix_struct%nrow_global >= elpa_neigvec_min) THEN
121 CALL cp_cfm_diag_elpa(matrix, eigenvectors, eigenvalues)
122 ELSE
123 CALL cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
124 END IF
125#if defined(__DLAF)
126 END IF
127#endif
128
129 CALL timestop(handle)
130
131 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 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 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work
150 COMPLEX(KIND=dp), DIMENSION(:, :), &
151 POINTER :: m
152 INTEGER :: handle, info, liwork, &
153 lrwork, lwork, n
154 INTEGER, DIMENSION(:), POINTER :: iwork
155 REAL(kind=dp), DIMENSION(:), POINTER :: rwork
156#if defined(__parallel)
157 INTEGER, DIMENSION(9) :: descm, descv
158 COMPLEX(KIND=dp), DIMENSION(:, :), &
159 POINTER :: v
160#endif
161#if defined (__HAS_IEEE_EXCEPTIONS)
162 LOGICAL, DIMENSION(5) :: halt
163#endif
164
165 CALL timeset(routinen, handle)
166
167 n = matrix%matrix_struct%nrow_global
168 m => matrix%local_data
169 ALLOCATE (iwork(1), rwork(1), work(1))
170 ! work space query
171 lwork = -1
172 lrwork = -1
173 liwork = -1
174
175#if defined(__parallel)
176 v => eigenvectors%local_data
177 descm(:) = matrix%matrix_struct%descriptor(:)
178 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 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 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 lrwork = ceiling(rwork(1)) + 1000000
186 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 DEALLOCATE (iwork, rwork, work)
196 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 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 DEALLOCATE (iwork, rwork, work)
217 IF (info /= 0) cpabort("Diagonalisation of a complex matrix failed")
218
219 CALL timestop(handle)
220
221 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
253 TYPE(cp_cfm_type), POINTER :: cfm_work
254 TYPE(cp_cfm_type), TARGET :: local_work
255
256 CALL timeset(routinen, handle)
257
258 CALL cp_cfm_get_info(matrix, nrow_global=nrow_global, ncol_global=ncol_global)
259 cpassert(nrow_global == ncol_global)
260 cpassert(threshold >= 0.0_dp)
261 IF (PRESENT(work)) THEN
262 cfm_work => work
263 ELSE
264 CALL cp_cfm_create(local_work, matrix%matrix_struct, set_zero=.true.)
265 cfm_work => local_work
266 END IF
267 ALLOCATE (eigenvalues(nrow_global), eigenvalue_scales(nrow_global))
268
269 CALL cp_cfm_heevd(matrix, cfm_work, eigenvalues)
270 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 DO i = 1, nrow_global
276 ! Keep zero modes quenched, including inverse powers with threshold=0.
277 IF (eigenvalues(i) < threshold .OR. eigenvalues(i) <= 0.0_dp) THEN
278 sqrt_power_eval = 0.0_dp
279 ndep = ndep + 1
280 ELSE
281 sqrt_power_eval = eigenvalues(i)**(0.5_dp*exponent)
282 END IF
283 eigenvalue_scales(i) = cmplx(sqrt_power_eval, 0.0_dp, kind=dp)
284 END DO
285 IF (PRESENT(n_dependent)) n_dependent = ndep
286
287 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 cfm_work, cfm_work, z_zero, matrix)
292
293 my_verbose = .false.
294 IF (PRESENT(verbose)) my_verbose = verbose
295 IF (matrix%matrix_struct%para_env%is_source() .AND. my_verbose) THEN
296 condition_number = huge(0.0_dp)
297 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 "CP_CFM_POWER: smallest eigenvalue:", eigenvalues(1), &
300 "CP_CFM_POWER: largest eigenvalue: ", eigenvalues(nrow_global), &
301 "CP_CFM_POWER: condition number: ", condition_number
302 IF (eigenvalues(1) <= 0.0_dp) THEN
303 WRITE (unit=cp_logger_get_default_unit_nr(), fmt="(/,T2,A)") &
304 "WARNING: matrix has a non-positive eigenvalue, tighten EPS_DEFAULT"
305 END IF
306 IF (condition_number > 1.0e12_dp) THEN
307 WRITE (unit=cp_logger_get_default_unit_nr(), fmt="(/,T2,A)") &
308 "WARNING: high condition number => possibly ill-conditioned matrix"
309 END IF
310 END IF
311
312 DEALLOCATE (eigenvalues, eigenvalue_scales)
313
314 IF (.NOT. PRESENT(work)) CALL cp_cfm_release(local_work)
315
316 CALL timestop(handle)
317
318 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 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 CALL timeset(routinen, handle)
347
348 IF (.NOT. diag_check_requested()) THEN
349 CALL timestop(handle)
350 RETURN
351 END IF
352
353 output_unit = default_output_unit
354 eps_warning = diag_check_warning_threshold()
355 eps_abort = 10.0_dp*eps_warning
356
357 nrow = eigenvectors%matrix_struct%nrow_global
358 ncol = min(eigenvectors%matrix_struct%ncol_global, nvec)
359
360 CALL cp_cfm_gemm("N", "N", nrow, ncol, nrow, z_one, overlap, eigenvectors, z_zero, scratch)
361 CALL cp_cfm_gemm("C", "N", ncol, ncol, nrow, z_one, eigenvectors, scratch, z_zero, overlap)
362
363 gold = z_zero
364 test = z_zero
365 eps = 0.0_dp
366
367#if defined(__parallel)
368 context => overlap%matrix_struct%context
369 myprow = context%mepos(1)
370 mypcol = context%mepos(2)
371 nprow = context%num_pe(1)
372 npcol = context%num_pe(2)
373 desca(:) = overlap%matrix_struct%descriptor(:)
374 outer: DO j = 1, ncol
375 DO i = 1, ncol
376 CALL infog2l(i, j, desca, nprow, npcol, myprow, mypcol, il, jl, iprow, ipcol)
377 IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
378 gold = merge(z_zero, z_one, i /= j)
379 test = overlap%local_data(il, jl)
380 eps = abs(test - gold)
381 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 IF (eps > eps_warning) THEN
397 IF (diag_type == fm_diag_type_scalapack) THEN
398 diag_type_name = "HEGVX"
399 ELSE IF (diag_type == fm_diag_type_cusolver) THEN
400 diag_type_name = "CUSOLVER"
401 ELSE IF (diag_type == fm_diag_type_elpa .AND. diag_lib_explicit) THEN
402 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 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 "The generalized eigenvectors returned by "//trim(diag_type_name)//" are not S-orthonormal", &
412 "Absolute deviation of matrix element (", i, ", ", j, ") is ", eps, &
413 "The deviation from the expected value ", real(gold, kind=dp), " is", eps
414 IF (eps > eps_abort) THEN
415 CALL cp_abort(__location__, &
416 "ERROR in "//routinen//": Check of generalized matrix diagonalization failed")
417 ELSE
418 cpwarn("Check of generalized matrix diagonalization failed in routine "//routinen)
419 END IF
420 END IF
421
422 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 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 CALL timeset(routinen, handle)
453
454 CALL cp_cfm_get_info(amatrix, nrow_global=nao)
455 ALLOCATE (evals(nao))
456 nmo = SIZE(eigenvalues)
457 check_eigenvectors = diag_check_requested()
458 use_lowest_subset = .false.
459 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 IF (check_eigenvectors) THEN
467 CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
468 CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
469 CALL cp_cfm_to_cfm(bmatrix, overlap_check)
470 END IF
471 CALL cp_cfm_geeig_scalapack(amatrix, bmatrix, work, evals(1:nmo))
472 IF (check_eigenvectors) THEN
473 CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
474 CALL cp_cfm_release(scratch_check)
475 CALL cp_cfm_release(overlap_check)
476 END IF
477#endif
478 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 IF (check_eigenvectors) THEN
483 CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
484 CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
485 CALL cp_cfm_to_cfm(bmatrix, overlap_check)
486 END IF
487 CALL cp_cfm_general_cusolver(amatrix, bmatrix, work, evals)
488 IF (check_eigenvectors) THEN
489 CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
490 CALL cp_cfm_release(scratch_check)
491 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 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 IF (check_eigenvectors) THEN
522 CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
523 CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
524 CALL cp_cfm_to_cfm(bmatrix, overlap_check)
525 END IF
526 CALL cp_cfm_geeig_scalapack(amatrix, bmatrix, work, evals)
527 IF (check_eigenvectors) THEN
528 CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
529 CALL cp_cfm_release(scratch_check)
530 CALL cp_cfm_release(overlap_check)
531 END IF
532#endif
533 ELSE
534 ! Cholesky decompose S=U(T)U
535 CALL cp_cfm_cholesky_decompose(bmatrix)
536 ! Invert to get U^(-1)
537 CALL cp_cfm_triangular_invert(bmatrix)
538 ! Reduce to get U^(-T) * H * U^(-1)
539 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
540 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
541 ! Diagonalize
542 CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
543 ! Restore vectors C = U^(-1) * C*
544 CALL cp_cfm_triangular_multiply(bmatrix, work)
545 END IF
546
547 CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
548 eigenvalues(1:nmo) = evals(1:nmo)
549
550 DEALLOCATE (evals)
551
552 CALL timestop(handle)
553
554 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 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 COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
576 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 INTEGER, DIMENSION(:), ALLOCATABLE :: iclustr, ifail, iwork
581 REAL(kind=dp) :: abstol
582 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 CALL timeset(routinen, handle)
595
596#if defined(__parallel)
597 n = amatrix%matrix_struct%nrow_global
598 neig = min(SIZE(eigenvalues), n)
599
600 IF (neig == 0) THEN
601 CALL timestop(handle)
602 RETURN
603 END IF
604
605 IF (amatrix%matrix_struct%nrow_block /= amatrix%matrix_struct%ncol_block) THEN
606 cpabort("ERROR in "//routinen//": Invalid blocksize (no square blocks) found")
607 END IF
608
609 a => amatrix%local_data
610 b => bmatrix%local_data
611 z => eigenvectors%local_data
612 desca(:) = amatrix%matrix_struct%descriptor(:)
613 descb(:) = bmatrix%matrix_struct%descriptor(:)
614 descz(:) = eigenvectors%matrix_struct%descriptor(:)
615
616 nprow = amatrix%matrix_struct%context%num_pe(1)
617 npcol = amatrix%matrix_struct%context%num_pe(2)
618 npe = nprow*npcol
619 nb = amatrix%matrix_struct%nrow_block
620 nn = max(n, nb, 2)
621 np0 = numroc(nn, nb, 0, 0, nprow)
622 mq0 = max(numroc(nn, nb, 0, 0, npcol), nb)
623
624 lwork = n + (np0 + mq0 + nb)*nb
625 lrwork = 4*n + max(5*nn, np0*mq0) + iceil(neig, npe)*nn + max(0, neig - 1)*n
626 liwork = 6*max(n, npe + 1, 4)
627
628 ALLOCATE (gap(npe))
629 gap = 0.0_dp
630 ALLOCATE (iclustr(2*npe))
631 iclustr = 0
632 ALLOCATE (ifail(n))
633 ifail = 0
634 ALLOCATE (iwork(liwork))
635 ALLOCATE (rwork(lrwork))
636 ALLOCATE (w(n))
637 ALLOCATE (work(lwork))
638
639 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 iclustr(1), gap(1), info)
649#if defined (__HAS_IEEE_EXCEPTIONS)
650 CALL ieee_set_halting_mode(ieee_all, halt)
651#endif
652
653 IF (info /= 0 .OR. m < neig .OR. nz < neig) THEN
654 cpabort("ERROR in PZHEGVX (ScaLAPACK), info="//trim(cp_to_string(info)))
655 END IF
656
657 eigenvalues(:) = 0.0_dp
658 eigenvalues(1:neig) = w(1:neig)
659
660 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 CALL timestop(handle)
670
671 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 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 CALL timeset(routinen, handle)
701
702 ! Test sizes
703 CALL cp_cfm_get_info(amatrix, nrow_global=nao)
704 nmo = SIZE(eigenvalues)
705 ALLOCATE (evals(nao), cevals(nao))
706
707 ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
708 CALL cp_cfm_scale(-z_one, bmatrix)
709 CALL cp_cfm_heevd(bmatrix, work, evals)
710 evals(:) = -evals(:)
711 nc = nao
712 DO i = 1, nao
713 IF (evals(i) < epseig) THEN
714 nc = i - 1
715 EXIT
716 END IF
717 END DO
718 cpassert(nc /= 0)
719
720 IF (nc /= nao) THEN
721 IF (nc < nmo) THEN
722 ! Copy NULL space definition to last vectors of eigenvectors (if needed)
723 ncol = nmo - nc
724 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 DO icol = nc + 1, nao
728 DO irow = 1, nao
729 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 evals(nc + 1:nao) = 1.0_dp
734 END IF
735 ! Calculate U*s**(-1/2)
736 cevals(:) = cmplx(1.0_dp/sqrt(evals(:)), 0.0_dp, kind=dp)
737 CALL cp_cfm_column_scale(work, cevals)
738 ! Reduce to get U^(-C) * H * U^(-1)
739 CALL cp_cfm_gemm("C", "N", nao, nao, nao, z_one, work, amatrix, z_zero, bmatrix)
740 CALL cp_cfm_gemm("N", "N", nao, nao, nao, z_one, bmatrix, work, z_zero, amatrix)
741 IF (nc /= nao) THEN
742 ! set diagonal values to save large value
743 DO icol = nc + 1, nao
744 CALL cp_cfm_set_element(amatrix, icol, icol, &
745 cmplx(set_removed_eigval_to, 0.0_dp, kind=dp))
746 END DO
747 END IF
748 ! Diagonalize
749 CALL cp_cfm_heevd(amatrix, bmatrix, evals)
750 eigenvalues(1:nmo) = evals(1:nmo)
751 nx = min(nc, nmo)
752 ! Restore vectors C = U^(-1) * C*
753 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 IF (PRESENT(nmo_retained)) nmo_retained = nc
758
759 DEALLOCATE (evals)
760
761 CALL timestop(handle)
762
763 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 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 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 CALL timeset(routinen, handle)
791
792 cpassert(product(amatrix%matrix_struct%context%num_pe) == 1)
793 n = amatrix%matrix_struct%nrow_global
794 nmo = min(n, SIZE(eigenvalues))
795 a => amatrix%local_data
796 b => bmatrix%local_data
797 z => eigenvectors%local_data
798
799 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 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 IF (info /= 0) THEN
811 CALL cp_abort(__location__, "Local ZHEGVD failed, info="//trim(cp_to_string(info)))
812 END IF
813
814 eigenvalues(1:nmo) = ws%evals(1:nmo)
815 z(1:n, 1:nmo) = a(1:n, 1:nmo)
816
817 CALL timestop(handle)
818
819 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 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 COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a, b, u, z
842 INTEGER :: handle, i, info, n, nc, nmo, nx
843
844 CALL timeset(routinen, handle)
845
846 cpassert(product(amatrix%matrix_struct%context%num_pe) == 1)
847 n = amatrix%matrix_struct%nrow_global
848 nmo = min(n, SIZE(eigenvalues))
849 a => amatrix%local_data
850 b => bmatrix%local_data
851 u => work%local_data
852 z => eigenvectors%local_data
853 CALL cp_cfm_local_workspace_prepare(ws, n, b)
854
855 b(1:n, 1:n) = -b(1:n, 1:n)
856 CALL cp_cfm_heevd_local(b, u, info, ws)
857 IF (info /= 0) THEN
858 CALL cp_abort(__location__, "Local overlap ZHEEVD failed, info="//trim(cp_to_string(info)))
859 END IF
860 ws%evals(:) = -ws%evals(:)
861 nc = n
862 DO i = 1, n
863 IF (ws%evals(i) < epseig) THEN
864 nc = i - 1
865 EXIT
866 END IF
867 END DO
868 cpassert(nc /= 0)
869
870 IF (nc < n) THEN
871 IF (nc < nmo) z(1:n, nc + 1:nmo) = u(1:n, nc + 1:nmo)
872 u(1:n, nc + 1:n) = z_zero
873 ws%evals(nc + 1:n) = 1.0_dp
874 END IF
875 DO i = 1, n
876 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 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 z_zero, a(1, 1), SIZE(a, 1))
883 IF (nc < n) THEN
884 DO i = nc + 1, n
885 a(i, i) = cmplx(10000.0_dp, 0.0_dp, kind=dp)
886 END DO
887 END IF
888
889 CALL cp_cfm_heevd_local(a, b, info, ws)
890 IF (info /= 0) THEN
891 CALL cp_abort(__location__, "Local Hamiltonian ZHEEVD failed, info="//trim(cp_to_string(info)))
892 END IF
893 eigenvalues(1:nmo) = ws%evals(1:nmo)
894 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 z_zero, z(1, 1), SIZE(z, 1))
897
898 CALL timestop(handle)
899
900 END SUBROUTINE cp_cfm_geeig_canon_local
901
902! **************************************************************************************************
903!> \brief Local LAPACK ZHEEVD with prepared scratch, without allocation or communication.
904!> \param matrix input matrix, overwritten by eigenvectors
905!> \param vectors optional copy of the eigenvectors
906!> \param info ...
907!> \param ws prepared private LAPACK scratch
908! **************************************************************************************************
909 SUBROUTINE cp_cfm_heevd_local(matrix, vectors, info, ws)
910
911 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix
912 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: vectors
913 INTEGER, INTENT(OUT) :: info
914
915 TYPE(cp_cfm_diag_workspace_type), INTENT(INOUT) :: ws
916 INTEGER :: n
917#if defined (__HAS_IEEE_EXCEPTIONS)
918 LOGICAL, DIMENSION(5) :: halt
919#endif
920
921 n = ws%n
922#if defined (__HAS_IEEE_EXCEPTIONS)
923 CALL ieee_get_halting_mode(ieee_all, halt)
924 CALL ieee_set_halting_mode(ieee_all, .false.)
925#endif
926 CALL zheevd('V', 'U', n, matrix(1, 1), SIZE(matrix, 1), ws%evals(1), &
927 ws%work(1), SIZE(ws%work), ws%rwork(1), SIZE(ws%rwork), ws%iwork(1), SIZE(ws%iwork), info)
928#if defined (__HAS_IEEE_EXCEPTIONS)
929 CALL ieee_set_halting_mode(ieee_all, halt)
930#endif
931 IF (PRESENT(vectors)) vectors(1:n, 1:n) = matrix(1:n, 1:n)
932
933 END SUBROUTINE cp_cfm_heevd_local
934
935! **************************************************************************************************
936!> \brief Query LAPACK once per matrix order/solver and retain private scratch across k points.
937!> \param ws caller-owned workspace
938!> \param n matrix order
939!> \param a query matrix; workspace queries do not read its values
940!> \param b optional overlap selects ZHEGVD; otherwise query ZHEEVD
941! **************************************************************************************************
942 SUBROUTINE cp_cfm_local_workspace_prepare(ws, n, a, b)
943 TYPE(cp_cfm_diag_workspace_type), INTENT(INOUT) :: ws
944 INTEGER, INTENT(IN) :: n
945 COMPLEX(KIND=dp), INTENT(INOUT) :: a(:, :)
946 COMPLEX(KIND=dp), INTENT(INOUT), OPTIONAL :: b(:, :)
947
948 COMPLEX(KIND=dp) :: work(1)
949 INTEGER :: info, iwork(1)
950 REAL(kind=dp) :: rwork(1)
951
952 IF (ws%n == n .AND. (ws%generalized .EQV. PRESENT(b))) RETURN
953 IF (ALLOCATED(ws%work)) DEALLOCATE (ws%work, ws%rwork, ws%iwork, ws%evals)
954 ALLOCATE (ws%evals(n))
955 IF (PRESENT(b)) THEN
956 CALL zhegvd(1, 'V', 'U', n, a(1, 1), SIZE(a, 1), b(1, 1), SIZE(b, 1), ws%evals(1), &
957 work(1), -1, rwork(1), -1, iwork(1), -1, info)
958 ELSE
959 CALL zheevd('V', 'U', n, a(1, 1), SIZE(a, 1), ws%evals(1), work(1), -1, rwork(1), -1, iwork(1), -1, info)
960 END IF
961 IF (info /= 0) THEN
962 CALL cp_abort(__location__, &
963 "Local LAPACK workspace query failed, info="//trim(cp_to_string(info)))
964 END IF
965 ALLOCATE (ws%work(max(1, ceiling(real(work(1), dp)))), &
966 ws%rwork(max(1, ceiling(rwork(1)))), ws%iwork(max(1, iwork(1))))
967 ws%n = n
968 ws%generalized = PRESENT(b)
969 END SUBROUTINE cp_cfm_local_workspace_prepare
970
971END MODULE cp_cfm_diag
cp_cfm_basic_linalg::cp_cfm_scale
Definition cp_cfm_basic_linalg.F:65
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:57
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_gemm
subroutine, public cp_cfm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
Performs one of the matrix-matrix operations: matrix_c = alpha * op1( matrix_a ) * op2( matrix_b ) + ...
Definition cp_cfm_basic_linalg.F:564
cp_cfm_basic_linalg::cp_cfm_triangular_multiply
subroutine, public cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
Multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
Definition cp_cfm_basic_linalg.F:983
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling, first_col, source)
Scales columns of the full matrix by corresponding factors.
Definition cp_cfm_basic_linalg.F:636
cp_cfm_basic_linalg::cp_cfm_triangular_invert
subroutine, public cp_cfm_triangular_invert(matrix_a, uplo, info_out)
Inverts a triangular matrix.
Definition cp_cfm_basic_linalg.F:1065
cp_cfm_cholesky
various cholesky decomposition related routines
Definition cp_cfm_cholesky.F:13
cp_cfm_cholesky::cp_cfm_cholesky_decompose
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Definition cp_cfm_cholesky.F:53
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_geeig_canon
subroutine, public cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig, nmo_retained)
General Eigenvalue Problem AX = BXE Use canonical orthogonalization.
Definition cp_cfm_diag.F:686
cp_cfm_diag::cp_cfm_heevd_local
subroutine, public cp_cfm_heevd_local(matrix, vectors, info, ws)
Local LAPACK ZHEEVD with prepared scratch, without allocation or communication.
Definition cp_cfm_diag.F:910
cp_cfm_diag::cp_cfm_geeig_local
subroutine, public cp_cfm_geeig_local(amatrix, bmatrix, eigenvectors, eigenvalues, ws)
Solve a generalized complex eigenproblem using the local LAPACK backend. This routine is restricted t...
Definition cp_cfm_diag.F:777
cp_cfm_diag::cp_cfm_local_workspace_prepare
subroutine, public cp_cfm_local_workspace_prepare(ws, n, a, b)
Query LAPACK once per matrix order/solver and retain private scratch across k points.
Definition cp_cfm_diag.F:943
cp_cfm_diag::cp_cfm_geeig_canon_local
subroutine, public cp_cfm_geeig_canon_local(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig, ws)
Canonical generalized complex diagonalization on a one-rank BLACS grid.
Definition cp_cfm_diag.F:832
cp_cfm_diag::cp_cfm_heevd
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition cp_cfm_diag.F:94
cp_cfm_diag::cp_cfm_power
subroutine, public cp_cfm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
Raise the positive spectral part of a Hermitian matrix to a real power. Eigenmodes with eigenvalues b...
Definition cp_cfm_diag.F:239
cp_cfm_diag::cp_cfm_geeig
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work, lowest_subset)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_cfm_diag.F:439
cp_cfm_dlaf_api
Definition cp_cfm_dlaf_api.F:8
cp_cfm_dlaf_api::cp_cfm_diag_dlaf
subroutine, public cp_cfm_diag_dlaf(matrix, eigenvectors, eigenvalues)
DLA-Future eigensolver for complex Hermitian matrices.
Definition cp_cfm_dlaf_api.F:119
cp_cfm_dlaf_api::cp_cfm_diag_gen_dlaf
subroutine, public cp_cfm_diag_gen_dlaf(amatrix, bmatrix, eigenvectors, eigenvalues)
DLA-Future generalized eigensolver for complex Hermitian matrices.
Definition cp_cfm_dlaf_api.F:157
cp_cfm_elpa
Wrapper for ELPA (complex matrices, i.e. cp_cfm_type).
Definition cp_cfm_elpa.F:11
cp_cfm_elpa::cp_cfm_diag_elpa
subroutine, public cp_cfm_diag_elpa(matrix, eigenvectors, eigenvalues)
Driver routine to diagonalize a CFM matrix with the ELPA library.
Definition cp_cfm_elpa.F:391
cp_cfm_elpa::is_elpa_c_broken
logical function, public is_elpa_c_broken()
Returns .TRUE. if the runtime check detected broken BLOCK2 kernels.
Definition cp_cfm_elpa.F:379
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:183
cp_cfm_types::cp_cfm_set_element
subroutine, public cp_cfm_set_element(matrix, irow_global, icol_global, alpha)
Set the matrix element (irow_global,icol_global) of the full matrix to alpha.
Definition cp_cfm_types.F:303
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:126
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:632
cp_dlaf_utils_api
Definition cp_dlaf_utils_api.F:8
cp_dlaf_utils_api::cp_dlaf_create_grid
subroutine, public cp_dlaf_create_grid(blacs_context)
Create DLA-Future grid from BLACS context.
Definition cp_dlaf_utils_api.F:75
cp_dlaf_utils_api::cp_dlaf_initialize
subroutine, public cp_dlaf_initialize()
Initialize DLA-Future and pika runtime.
Definition cp_dlaf_utils_api.F:38
cp_fm_cusolver_api
Wrapper for cuSOLVERMp.
Definition cp_fm_cusolver_api.F:12
cp_fm_cusolver_api::cp_cfm_general_cusolver
subroutine, public cp_cfm_general_cusolver(amatrix, bmatrix, eigenvectors, eigenvalues)
Driver routine to solve generalized complex eigenvalue problem A*x = lambda*B*x with cuSOLVERMp.
Definition cp_fm_cusolver_api.F:193
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::fm_diag_type_cusolver
integer, parameter, public fm_diag_type_cusolver
Definition cp_fm_diag.F:112
cp_fm_diag::set_removed_eigval_to
real(kind=dp), parameter, public set_removed_eigval_to
Definition cp_fm_diag.F:89
cp_fm_diag::fm_diag_type_dlaf
integer, parameter, public fm_diag_type_dlaf
Definition cp_fm_diag.F:112
cp_fm_diag::fm_diag_type_scalapack
integer, parameter, public fm_diag_type_scalapack
Definition cp_fm_diag.F:112
cp_fm_diag::direct_generalized_diagonalization
logical, save, public direct_generalized_diagonalization
Definition cp_fm_diag.F:106
cp_fm_diag::elpa_neigvec_min
integer, save, public elpa_neigvec_min
Definition cp_fm_diag.F:99
cp_fm_diag::diag_check_warning_threshold
real(kind=dp) function, public diag_check_warning_threshold()
Return the warning threshold for diagonalization checks.
Definition cp_fm_diag.F:331
cp_fm_diag::diag_type
integer, save, public diag_type
Definition cp_fm_diag.F:93
cp_fm_diag::dlaf_neigvec_min
integer, save, public dlaf_neigvec_min
Definition cp_fm_diag.F:105
cp_fm_diag::fm_diag_type_elpa
integer, parameter, public fm_diag_type_elpa
Definition cp_fm_diag.F:112
cp_fm_diag::cusolver_n_min
integer, parameter, public cusolver_n_min
Definition cp_fm_diag.F:102
cp_fm_diag::diag_check_requested
logical function, public diag_check_requested()
Return whether diagonalization checks should be performed.
Definition cp_fm_diag.F:316
cp_fm_diag::diag_lib_explicit
logical, save, public diag_lib_explicit
Definition cp_fm_diag.F:96
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:583
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::default_output_unit
integer, parameter, public default_output_unit
Definition machine.F:48
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:144
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_diag::cp_cfm_diag_workspace_type
Definition cp_cfm_diag.F:72
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:70