(git:48c3be8)
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! **************************************************************************************************
22 USE cp_cfm_types, ONLY: cp_cfm_create, &
31 diag_type, &
39 USE cp_cfm_elpa, ONLY: cp_cfm_diag_elpa, &
42#if defined(__DLAF)
47#endif
49 USE kinds, ONLY: default_string_length, &
50 dp
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
68
69CONTAINS
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 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 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
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.
107 .NOT. is_elpa_c_broken() .AND. &
108 matrix%matrix_struct%nrow_global >= elpa_neigvec_min) THEN
109 CALL cp_cfm_diag_elpa(matrix, eigenvectors, eigenvalues)
110 ELSE
111 CALL cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
112 END IF
113#if defined(__DLAF)
114 END IF
115#endif
116
117 CALL timestop(handle)
118
119 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 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 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work
138 COMPLEX(KIND=dp), DIMENSION(:, :), &
139 POINTER :: m
140 INTEGER :: handle, info, liwork, &
141 lrwork, lwork, n
142 INTEGER, DIMENSION(:), POINTER :: iwork
143 REAL(kind=dp), DIMENSION(:), POINTER :: rwork
144#if defined(__parallel)
145 INTEGER, DIMENSION(9) :: descm, descv
146 COMPLEX(KIND=dp), DIMENSION(:, :), &
147 POINTER :: v
148#endif
149#if defined (__HAS_IEEE_EXCEPTIONS)
150 LOGICAL, DIMENSION(5) :: halt
151#endif
152
153 CALL timeset(routinen, handle)
154
155 n = matrix%matrix_struct%nrow_global
156 m => matrix%local_data
157 ALLOCATE (iwork(1), rwork(1), work(1))
158 ! work space query
159 lwork = -1
160 lrwork = -1
161 liwork = -1
162
163#if defined(__parallel)
164 v => eigenvectors%local_data
165 descm(:) = matrix%matrix_struct%descriptor(:)
166 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 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 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 lrwork = ceiling(rwork(1)) + 1000000
174 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 DEALLOCATE (iwork, rwork, work)
184 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 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 DEALLOCATE (iwork, rwork, work)
205 IF (info /= 0) cpabort("Diagonalisation of a complex matrix failed")
206
207 CALL timestop(handle)
208
209 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 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 CALL timeset(routinen, handle)
238
239 IF (.NOT. diag_check_requested()) THEN
240 CALL timestop(handle)
241 RETURN
242 END IF
243
244 output_unit = default_output_unit
245 eps_warning = diag_check_warning_threshold()
246 eps_abort = 10.0_dp*eps_warning
247
248 nrow = eigenvectors%matrix_struct%nrow_global
249 ncol = min(eigenvectors%matrix_struct%ncol_global, nvec)
250
251 CALL cp_cfm_gemm("N", "N", nrow, ncol, nrow, z_one, overlap, eigenvectors, z_zero, scratch)
252 CALL cp_cfm_gemm("C", "N", ncol, ncol, nrow, z_one, eigenvectors, scratch, z_zero, overlap)
253
254 gold = z_zero
255 test = z_zero
256 eps = 0.0_dp
257
258#if defined(__parallel)
259 context => overlap%matrix_struct%context
260 myprow = context%mepos(1)
261 mypcol = context%mepos(2)
262 nprow = context%num_pe(1)
263 npcol = context%num_pe(2)
264 desca(:) = overlap%matrix_struct%descriptor(:)
265 outer: DO j = 1, ncol
266 DO i = 1, ncol
267 CALL infog2l(i, j, desca, nprow, npcol, myprow, mypcol, il, jl, iprow, ipcol)
268 IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
269 gold = merge(z_zero, z_one, i /= j)
270 test = overlap%local_data(il, jl)
271 eps = abs(test - gold)
272 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 IF (eps > eps_warning) THEN
289 diag_type_name = "HEGVX"
290 ELSE IF (diag_type == fm_diag_type_cusolver) THEN
291 diag_type_name = "CUSOLVER"
292 ELSE IF (diag_type == fm_diag_type_elpa .AND. diag_lib_explicit) THEN
293 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 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 "The generalized eigenvectors returned by "//trim(diag_type_name)//" are not S-orthonormal", &
303 "Absolute deviation of matrix element (", i, ", ", j, ") is ", eps, &
304 "The deviation from the expected value ", real(gold, kind=dp), " is", eps
305 IF (eps > eps_abort) THEN
306 cpabort("ERROR in "//routinen//": Check of generalized matrix diagonalization failed")
307 ELSE
308 cpwarn("Check of generalized matrix diagonalization failed in routine "//routinen)
309 END IF
310 END IF
311
312 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 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 CALL timeset(routinen, handle)
343
344 CALL cp_cfm_get_info(amatrix, nrow_global=nao)
345 ALLOCATE (evals(nao))
346 nmo = SIZE(eigenvalues)
347 check_eigenvectors = diag_check_requested()
348 use_lowest_subset = .false.
349 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 IF (check_eigenvectors) THEN
357 CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
358 CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
359 CALL cp_cfm_to_cfm(bmatrix, overlap_check)
360 END IF
361 CALL cp_cfm_geeig_scalapack(amatrix, bmatrix, work, evals(1:nmo))
362 IF (check_eigenvectors) THEN
363 CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
364 CALL cp_cfm_release(scratch_check)
365 CALL cp_cfm_release(overlap_check)
366 END IF
367#endif
369 nao >= cusolver_n_min) THEN
370 ! Use cuSolverMP generalized eigenvalue solver without a CP2K-side
371 ! Cholesky reduction.
372 IF (check_eigenvectors) THEN
373 CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
374 CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
375 CALL cp_cfm_to_cfm(bmatrix, overlap_check)
376 END IF
377 CALL cp_cfm_general_cusolver(amatrix, bmatrix, work, evals)
378 IF (check_eigenvectors) THEN
379 CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
380 CALL cp_cfm_release(scratch_check)
381 CALL cp_cfm_release(overlap_check)
382 END IF
383#if defined(__DLAF)
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)
409 ! Use ScaLAPACK generalized eigenvalue solver without a CP2K-side
410 ! Cholesky reduction.
411 IF (check_eigenvectors) THEN
412 CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
413 CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
414 CALL cp_cfm_to_cfm(bmatrix, overlap_check)
415 END IF
416 CALL cp_cfm_geeig_scalapack(amatrix, bmatrix, work, evals)
417 IF (check_eigenvectors) THEN
418 CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
419 CALL cp_cfm_release(scratch_check)
420 CALL cp_cfm_release(overlap_check)
421 END IF
422#endif
423 ELSE
424 ! Cholesky decompose S=U(T)U
425 CALL cp_cfm_cholesky_decompose(bmatrix)
426 ! Invert to get U^(-1)
427 CALL cp_cfm_triangular_invert(bmatrix)
428 ! Reduce to get U^(-T) * H * U^(-1)
429 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
430 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
431 ! Diagonalize
432 CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
433 ! Restore vectors C = U^(-1) * C*
434 CALL cp_cfm_triangular_multiply(bmatrix, work)
435 END IF
436
437 CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
438 eigenvalues(1:nmo) = evals(1:nmo)
439
440 DEALLOCATE (evals)
441
442 CALL timestop(handle)
443
444 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 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 COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
466 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 INTEGER, DIMENSION(:), ALLOCATABLE :: iclustr, ifail, iwork
471 REAL(kind=dp) :: abstol
472 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 CALL timeset(routinen, handle)
485
486#if defined(__parallel)
487 n = amatrix%matrix_struct%nrow_global
488 neig = min(SIZE(eigenvalues), n)
489
490 IF (neig == 0) THEN
491 CALL timestop(handle)
492 RETURN
493 END IF
494
495 IF (amatrix%matrix_struct%nrow_block /= amatrix%matrix_struct%ncol_block) THEN
496 cpabort("ERROR in "//routinen//": Invalid blocksize (no square blocks) found")
497 END IF
498
499 a => amatrix%local_data
500 b => bmatrix%local_data
501 z => eigenvectors%local_data
502 desca(:) = amatrix%matrix_struct%descriptor(:)
503 descb(:) = bmatrix%matrix_struct%descriptor(:)
504 descz(:) = eigenvectors%matrix_struct%descriptor(:)
505
506 nprow = amatrix%matrix_struct%context%num_pe(1)
507 npcol = amatrix%matrix_struct%context%num_pe(2)
508 npe = nprow*npcol
509 nb = amatrix%matrix_struct%nrow_block
510 nn = max(n, nb, 2)
511 np0 = numroc(nn, nb, 0, 0, nprow)
512 mq0 = max(numroc(nn, nb, 0, 0, npcol), nb)
513
514 lwork = n + (np0 + mq0 + nb)*nb
515 lrwork = 4*n + max(5*nn, np0*mq0) + iceil(neig, npe)*nn + max(0, neig - 1)*n
516 liwork = 6*max(n, npe + 1, 4)
517
518 ALLOCATE (gap(npe))
519 gap = 0.0_dp
520 ALLOCATE (iclustr(2*npe))
521 iclustr = 0
522 ALLOCATE (ifail(n))
523 ifail = 0
524 ALLOCATE (iwork(liwork))
525 ALLOCATE (rwork(lrwork))
526 ALLOCATE (w(n))
527 ALLOCATE (work(lwork))
528
529 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 iclustr(1), gap(1), info)
539#if defined (__HAS_IEEE_EXCEPTIONS)
540 CALL ieee_set_halting_mode(ieee_all, halt)
541#endif
542
543 IF (info /= 0 .OR. m < neig .OR. nz < neig) THEN
544 cpabort("ERROR in PZHEGVX (ScaLAPACK), info="//trim(cp_to_string(info)))
545 END IF
546
547 eigenvalues(:) = 0.0_dp
548 eigenvalues(1:neig) = w(1:neig)
549
550 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 CALL timestop(handle)
560
561 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 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 CALL timeset(routinen, handle)
591
592 ! Test sizes
593 CALL cp_cfm_get_info(amatrix, nrow_global=nao)
594 nmo = SIZE(eigenvalues)
595 ALLOCATE (evals(nao), cevals(nao))
596
597 ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
598 CALL cp_cfm_scale(-z_one, bmatrix)
599 CALL cp_cfm_heevd(bmatrix, work, evals)
600 evals(:) = -evals(:)
601 nc = nao
602 DO i = 1, nao
603 IF (evals(i) < epseig) THEN
604 nc = i - 1
605 EXIT
606 END IF
607 END DO
608 cpassert(nc /= 0)
609
610 IF (nc /= nao) THEN
611 IF (nc < nmo) THEN
612 ! Copy NULL space definition to last vectors of eigenvectors (if needed)
613 ncol = nmo - nc
614 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 DO icol = nc + 1, nao
618 DO irow = 1, nao
619 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 evals(nc + 1:nao) = 1.0_dp
624 END IF
625 ! Calculate U*s**(-1/2)
626 cevals(:) = cmplx(1.0_dp/sqrt(evals(:)), 0.0_dp, kind=dp)
627 CALL cp_cfm_column_scale(work, cevals)
628 ! Reduce to get U^(-C) * H * U^(-1)
629 CALL cp_cfm_gemm("C", "N", nao, nao, nao, z_one, work, amatrix, z_zero, bmatrix)
630 CALL cp_cfm_gemm("N", "N", nao, nao, nao, z_one, bmatrix, work, z_zero, amatrix)
631 IF (nc /= nao) THEN
632 ! set diagonal values to save large value
633 DO icol = nc + 1, nao
634 CALL cp_cfm_set_element(amatrix, icol, icol, &
635 cmplx(set_removed_eigval_to, 0.0_dp, kind=dp))
636 END DO
637 END IF
638 ! Diagonalize
639 CALL cp_cfm_heevd(amatrix, bmatrix, evals)
640 eigenvalues(1:nmo) = evals(1:nmo)
641 nx = min(nc, nmo)
642 ! Restore vectors C = U^(-1) * C*
643 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 IF (PRESENT(nmo_retained)) nmo_retained = nc
648
649 DEALLOCATE (evals)
650
651 CALL timestop(handle)
652
653 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 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 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
673 COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a, b, z
674 INTEGER :: handle, info, liwork, lrwork, lwork, &
675 n, nmo
676 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
677 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals, rwork
678#if defined (__HAS_IEEE_EXCEPTIONS)
679 LOGICAL, DIMENSION(5) :: halt
680#endif
681
682 CALL timeset(routinen, handle)
683
684 cpassert(product(amatrix%matrix_struct%context%num_pe) == 1)
685 n = amatrix%matrix_struct%nrow_global
686 nmo = min(n, SIZE(eigenvalues))
687 a => amatrix%local_data
688 b => bmatrix%local_data
689 z => eigenvectors%local_data
690
691 ALLOCATE (evals(n), iwork(1), rwork(1), work(1))
692 lwork = -1
693 lrwork = -1
694 liwork = -1
695 CALL zhegvd(1, 'V', 'U', n, a(1, 1), SIZE(a, 1), b(1, 1), SIZE(b, 1), evals(1), &
696 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
697 IF (info /= 0) cpabort("Local ZHEGVD workspace query failed, info="//trim(cp_to_string(info)))
698 lwork = max(1, ceiling(real(work(1), kind=dp)))
699 lrwork = max(1, ceiling(rwork(1)))
700 liwork = max(1, iwork(1))
701 DEALLOCATE (iwork, rwork, work)
702 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 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 IF (info /= 0) cpabort("Local ZHEGVD failed, info="//trim(cp_to_string(info)))
714
715 eigenvalues(1:nmo) = evals(1:nmo)
716 z(1:n, 1:nmo) = a(1:n, 1:nmo)
717
718 DEALLOCATE (evals, iwork, rwork, work)
719 CALL timestop(handle)
720
721 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 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 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 CALL timeset(routinen, handle)
746
747 cpassert(product(amatrix%matrix_struct%context%num_pe) == 1)
748 n = amatrix%matrix_struct%nrow_global
749 nmo = min(n, SIZE(eigenvalues))
750 a => amatrix%local_data
751 b => bmatrix%local_data
752 u => work%local_data
753 z => eigenvectors%local_data
754 ALLOCATE (evals(n))
755
756 b(1:n, 1:n) = -b(1:n, 1:n)
757 CALL cp_cfm_heevd_local(b, u, evals, info)
758 IF (info /= 0) cpabort("Local overlap ZHEEVD failed, info="//trim(cp_to_string(info)))
759 evals(:) = -evals(:)
760 nc = n
761 DO i = 1, n
762 IF (evals(i) < epseig) THEN
763 nc = i - 1
764 EXIT
765 END IF
766 END DO
767 cpassert(nc /= 0)
768
769 IF (nc < n) THEN
770 IF (nc < nmo) z(1:n, nc + 1:nmo) = u(1:n, nc + 1:nmo)
771 u(1:n, nc + 1:n) = z_zero
772 evals(nc + 1:n) = 1.0_dp
773 END IF
774 DO i = 1, n
775 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 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 z_zero, a(1, 1), SIZE(a, 1))
782 IF (nc < n) THEN
783 DO i = nc + 1, n
784 a(i, i) = cmplx(10000.0_dp, 0.0_dp, kind=dp)
785 END DO
786 END IF
787
788 CALL cp_cfm_heevd_local(a, b, evals, info)
789 IF (info /= 0) cpabort("Local Hamiltonian ZHEEVD failed, info="//trim(cp_to_string(info)))
790 eigenvalues(1:nmo) = evals(1:nmo)
791 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 z_zero, z(1, 1), SIZE(z, 1))
794
795 DEALLOCATE (evals)
796 CALL timestop(handle)
797
798 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 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 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
814 INTEGER :: liwork, lrwork, lwork, n
815 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
816 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rwork
817#if defined (__HAS_IEEE_EXCEPTIONS)
818 LOGICAL, DIMENSION(5) :: halt
819#endif
820
821 n = SIZE(eigenvalues)
822 ALLOCATE (iwork(1), rwork(1), work(1))
823 lwork = -1
824 lrwork = -1
825 liwork = -1
826 CALL zheevd('V', 'U', n, matrix(1, 1), SIZE(matrix, 1), eigenvalues(1), &
827 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
828 IF (info /= 0) cpabort("Local ZHEEVD workspace query failed, info="//trim(cp_to_string(info)))
829 lwork = max(1, ceiling(real(work(1), kind=dp)))
830 lrwork = max(1, ceiling(rwork(1)))
831 liwork = max(1, iwork(1))
832 DEALLOCATE (iwork, rwork, work)
833 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 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 vectors(1:n, 1:n) = matrix(1:n, 1:n)
844 DEALLOCATE (iwork, rwork, work)
845
846 END SUBROUTINE cp_cfm_heevd_local
847
848END MODULE cp_cfm_diag
methods related to the blacs parallel environment
Basic linear algebra operations for complex full matrices.
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 ) + ...
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...
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
subroutine, public cp_cfm_triangular_invert(matrix_a, uplo, info_out)
Inverts a triangular matrix.
various cholesky decomposition related routines
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...
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
subroutine, public cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig, nmo_retained)
General Eigenvalue Problem AX = BXE Use canonical orthogonalization.
subroutine, public cp_cfm_geeig_canon_local(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
Canonical generalized complex diagonalization on a one-rank BLACS grid.
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition cp_cfm_diag.F:82
subroutine, public cp_cfm_geeig_local(amatrix, bmatrix, eigenvectors, eigenvalues)
Solve a generalized complex eigenproblem using the local LAPACK backend. This routine is restricted t...
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work, lowest_subset)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
subroutine, public cp_cfm_diag_dlaf(matrix, eigenvectors, eigenvalues)
DLA-Future eigensolver for complex Hermitian matrices.
subroutine, public cp_cfm_diag_gen_dlaf(amatrix, bmatrix, eigenvectors, eigenvalues)
DLA-Future generalized eigensolver for complex Hermitian matrices.
Wrapper for ELPA (complex matrices, i.e. cp_cfm_type).
Definition cp_cfm_elpa.F:11
subroutine, public cp_cfm_diag_elpa(matrix, eigenvectors, eigenvalues)
Driver routine to diagonalize a CFM matrix with the ELPA library.
logical function, public is_elpa_c_broken()
Returns .TRUE. if the runtime check detected broken BLOCK2 kernels.
Represents a complex full matrix distributed on many processors.
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
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.
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
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.
subroutine, public cp_dlaf_create_grid(blacs_context)
Create DLA-Future grid from BLACS context.
subroutine, public cp_dlaf_initialize()
Initialize DLA-Future and pika runtime.
Wrapper for cuSOLVERMp.
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.
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
integer, parameter, public fm_diag_type_cusolver
Definition cp_fm_diag.F:112
real(kind=dp), parameter, public set_removed_eigval_to
Definition cp_fm_diag.F:89
integer, parameter, public fm_diag_type_dlaf
Definition cp_fm_diag.F:112
integer, parameter, public fm_diag_type_scalapack
Definition cp_fm_diag.F:112
logical, save, public direct_generalized_diagonalization
Definition cp_fm_diag.F:106
integer, save, public elpa_neigvec_min
Definition cp_fm_diag.F:99
real(kind=dp) function, public diag_check_warning_threshold()
Return the warning threshold for diagonalization checks.
Definition cp_fm_diag.F:331
integer, save, public diag_type
Definition cp_fm_diag.F:93
integer, save, public dlaf_neigvec_min
Definition cp_fm_diag.F:105
integer, parameter, public fm_diag_type_elpa
Definition cp_fm_diag.F:112
integer, parameter, public cusolver_n_min
Definition cp_fm_diag.F:102
logical function, public diag_check_requested()
Return whether diagonalization checks should be performed.
Definition cp_fm_diag.F:316
logical, save, public diag_lib_explicit
Definition cp_fm_diag.F:96
various routines to log and control the output. The idea is that decisions about where to log should ...
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
integer, parameter, public default_output_unit
Definition machine.F:46
Definition of mathematical constants and functions.
complex(kind=dp), parameter, public z_one
complex(kind=dp), parameter, public z_zero
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Represent a complex full matrix.