(git:85b8a9b)
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_type, &
31 direct_generalized_diagonalization, &
32 cusolver_n_min, &
33 fm_diag_type_cusolver, &
34 fm_diag_type_scalapack
35 USE cp_fm_cusolver_api, ONLY: cp_cfm_general_cusolver
36#if defined(__DLAF)
37 USE cp_cfm_dlaf_api, ONLY: cp_cfm_diag_gen_dlaf, &
38 cp_cfm_diag_dlaf
39 USE cp_dlaf_utils_api, ONLY: cp_dlaf_initialize, cp_dlaf_create_grid
40 USE cp_fm_diag, ONLY: dlaf_neigvec_min, fm_diag_type_dlaf
41#endif
42 USE cp_log_handling, ONLY: cp_to_string
43 USE kinds, ONLY: default_string_length, &
44 dp
45 USE machine, ONLY: default_output_unit
46 USE mathconstants, ONLY: z_one, &
47 z_zero
48#if defined (__HAS_IEEE_EXCEPTIONS)
49 USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
50 ieee_set_halting_mode, &
51 ieee_all
52#endif
53#include "../base/base_uses.f90"
54
55 IMPLICIT NONE
56 PRIVATE
57
58 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_diag'
59
60 PUBLIC :: cp_cfm_heevd, cp_cfm_geeig, cp_cfm_geeig_canon
61
62CONTAINS
63
64! **************************************************************************************************
65!> \brief Perform a diagonalisation of a complex matrix
66!> \param matrix ...
67!> \param eigenvectors ...
68!> \param eigenvalues ...
69!> \par History
70!> 12.2024 Added DLA-Future support [Rocco Meli]
71!> \author Joost VandeVondele
72! **************************************************************************************************
73 SUBROUTINE cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
74
75 TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
76 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
77
78 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_heevd'
79
80 INTEGER :: handle
81
82 CALL timeset(routinen, handle)
83
84#if defined(__DLAF)
85 IF (diag_type == fm_diag_type_dlaf .AND. matrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN
86 ! Initialize DLA-Future on-demand; if already initialized, does nothing
87 CALL cp_dlaf_initialize()
88
89 ! Create DLAF grid from BLACS context; if already present, does nothing
90 CALL cp_dlaf_create_grid(matrix%matrix_struct%context%get_handle())
91
92 CALL cp_cfm_diag_dlaf(matrix, eigenvectors, eigenvalues)
93 ELSE
94#endif
95 CALL cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
96#if defined(__DLAF)
97 END IF
98#endif
99
100 CALL timestop(handle)
101
102 END SUBROUTINE cp_cfm_heevd
103
104! **************************************************************************************************
105!> \brief Perform a diagonalisation of a complex matrix
106!> \param matrix ...
107!> \param eigenvectors ...
108!> \param eigenvalues ...
109!> \par History
110!> - (De)Allocation checks updated (15.02.2011,MK)
111!> \author Joost VandeVondele
112! **************************************************************************************************
113 SUBROUTINE cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
114
115 TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
116 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
117
118 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_heevd_base'
119
120 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work
121 COMPLEX(KIND=dp), DIMENSION(:, :), &
122 POINTER :: m
123 INTEGER :: handle, info, liwork, &
124 lrwork, lwork, n
125 INTEGER, DIMENSION(:), POINTER :: iwork
126 REAL(kind=dp), DIMENSION(:), POINTER :: rwork
127#if defined(__parallel)
128 INTEGER, DIMENSION(9) :: descm, descv
129 COMPLEX(KIND=dp), DIMENSION(:, :), &
130 POINTER :: v
131#endif
132#if defined (__HAS_IEEE_EXCEPTIONS)
133 LOGICAL, DIMENSION(5) :: halt
134#endif
135
136 CALL timeset(routinen, handle)
137
138 n = matrix%matrix_struct%nrow_global
139 m => matrix%local_data
140 ALLOCATE (iwork(1), rwork(1), work(1))
141 ! work space query
142 lwork = -1
143 lrwork = -1
144 liwork = -1
145
146#if defined(__parallel)
147 v => eigenvectors%local_data
148 descm(:) = matrix%matrix_struct%descriptor(:)
149 descv(:) = eigenvectors%matrix_struct%descriptor(:)
150 CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
151 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
152 ! The work space query for lwork does not return always sufficiently large values.
153 ! Let's add some margin to avoid crashes.
154 lwork = ceiling(real(work(1), kind=dp)) + 1000
155 ! needed to correct for a bug in scalapack, unclear how much the right number is
156 lrwork = ceiling(rwork(1)) + 1000000
157 liwork = iwork(1)
158#else
159 CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
160 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
161 lwork = ceiling(real(work(1), kind=dp))
162 lrwork = ceiling(rwork(1))
163 liwork = iwork(1)
164#endif
165
166 DEALLOCATE (iwork, rwork, work)
167 ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
168
169! (Sca-)LAPACK takes advantage of IEEE754 exceptions for speedup.
170! Therefore, we disable floating point traps temporarily.
171#if defined (__HAS_IEEE_EXCEPTIONS)
172 CALL ieee_get_halting_mode(ieee_all, halt)
173 CALL ieee_set_halting_mode(ieee_all, .false.)
174#endif
175#if defined(__parallel)
176 CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
177 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
178#else
179 CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
180 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
181 eigenvectors%local_data = matrix%local_data
182#endif
183#if defined (__HAS_IEEE_EXCEPTIONS)
184 CALL ieee_set_halting_mode(ieee_all, halt)
185#endif
186
187 DEALLOCATE (iwork, rwork, work)
188 IF (info /= 0) cpabort("Diagonalisation of a complex matrix failed")
189
190 CALL timestop(handle)
191
192 END SUBROUTINE cp_cfm_heevd_base
193
194! **************************************************************************************************
195!> \brief Check C^H*S*C = I for a generalized complex eigenvalue problem.
196!> \param overlap original overlap matrix S; used as work matrix and overwritten
197!> \param eigenvectors eigenvectors C to be checked
198!> \param scratch work matrix
199!> \param nvec ...
200! **************************************************************************************************
201 SUBROUTINE check_generalized_diag(overlap, eigenvectors, scratch, nvec)
202
203 TYPE(cp_cfm_type), INTENT(IN) :: eigenvectors
204 TYPE(cp_cfm_type), INTENT(INOUT) :: overlap, scratch
205 INTEGER, INTENT(IN) :: nvec
206
207 CHARACTER(LEN=*), PARAMETER :: routinen = 'check_generalized_diag'
208
209 CHARACTER(LEN=default_string_length) :: diag_type_name
210 COMPLEX(KIND=dp) :: gold, test
211 INTEGER :: handle, i, j, ncol, nrow, output_unit
212 REAL(kind=dp) :: eps, eps_abort, eps_warning
213#if defined(__parallel)
214 TYPE(cp_blacs_env_type), POINTER :: context
215 INTEGER :: il, jl, ipcol, iprow, &
216 mypcol, myprow, npcol, nprow
217 INTEGER, DIMENSION(9) :: desca
218#endif
219
220 CALL timeset(routinen, handle)
221
222 IF (.NOT. diag_check_requested()) THEN
223 CALL timestop(handle)
224 RETURN
225 END IF
226
227 output_unit = default_output_unit
228 eps_warning = diag_check_warning_threshold()
229 eps_abort = 10.0_dp*eps_warning
230
231 nrow = eigenvectors%matrix_struct%nrow_global
232 ncol = min(eigenvectors%matrix_struct%ncol_global, nvec)
233
234 CALL cp_cfm_gemm("N", "N", nrow, ncol, nrow, z_one, overlap, eigenvectors, z_zero, scratch)
235 CALL cp_cfm_gemm("C", "N", ncol, ncol, nrow, z_one, eigenvectors, scratch, z_zero, overlap)
236
237 gold = z_zero
238 test = z_zero
239 eps = 0.0_dp
240
241#if defined(__parallel)
242 context => overlap%matrix_struct%context
243 myprow = context%mepos(1)
244 mypcol = context%mepos(2)
245 nprow = context%num_pe(1)
246 npcol = context%num_pe(2)
247 desca(:) = overlap%matrix_struct%descriptor(:)
248 outer: DO j = 1, ncol
249 DO i = 1, ncol
250 CALL infog2l(i, j, desca, nprow, npcol, myprow, mypcol, il, jl, iprow, ipcol)
251 IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
252 gold = merge(z_zero, z_one, i /= j)
253 test = overlap%local_data(il, jl)
254 eps = abs(test - gold)
255 IF (eps > eps_warning) EXIT outer
256 END IF
257 END DO
258 END DO outer
259#else
260 outer: DO j = 1, ncol
261 DO i = 1, ncol
262 gold = merge(z_zero, z_one, i /= j)
263 test = overlap%local_data(i, j)
264 eps = abs(test - gold)
265 IF (eps > eps_warning) EXIT outer
266 END DO
267 END DO outer
268#endif
269
270 IF (eps > eps_warning) THEN
271 IF (diag_type == fm_diag_type_scalapack) THEN
272 diag_type_name = "HEGVX"
273 ELSE IF (diag_type == fm_diag_type_cusolver) THEN
274 diag_type_name = "CUSOLVER"
275#if defined(__DLAF)
276 ELSE IF (diag_type == fm_diag_type_dlaf) THEN
277 diag_type_name = "DLAF"
278#endif
279 ELSE
280 diag_type_name = "generalized eigensolver"
281 END IF
282 WRITE (unit=output_unit, fmt="(/,T2,A,/,T2,A,I0,A,I0,A,ES10.3,/,T2,A,F0.0,A,ES10.3)") &
283 "The generalized eigenvectors returned by "//trim(diag_type_name)//" are not S-orthonormal", &
284 "Absolute deviation of matrix element (", i, ", ", j, ") is ", eps, &
285 "The deviation from the expected value ", real(gold, kind=dp), " is", eps
286 IF (eps > eps_abort) THEN
287 cpabort("ERROR in "//routinen//": Check of generalized matrix diagonalization failed")
288 ELSE
289 cpwarn("Check of generalized matrix diagonalization failed in routine "//routinen)
290 END IF
291 END IF
292
293 CALL timestop(handle)
294
295 END SUBROUTINE check_generalized_diag
296
297! **************************************************************************************************
298!> \brief General Eigenvalue Problem AX = BXE
299!> Single option version: Cholesky decomposition of B
300!> \param amatrix ...
301!> \param bmatrix ...
302!> \param eigenvectors ...
303!> \param eigenvalues ...
304!> \param work ...
305!> \par History
306!> 12.2024 Added DLA-Future support [Rocco Meli]
307! **************************************************************************************************
308 SUBROUTINE cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
309
310 TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
311 REAL(kind=dp), DIMENSION(:) :: eigenvalues
312 TYPE(cp_cfm_type), INTENT(IN) :: work
313
314 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_geeig'
315
316 INTEGER :: handle, nao, nmo
317 LOGICAL :: check_eigenvectors
318 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
319 TYPE(cp_cfm_type) :: overlap_check, scratch_check
320
321 CALL timeset(routinen, handle)
322
323 CALL cp_cfm_get_info(amatrix, nrow_global=nao)
324 ALLOCATE (evals(nao))
325 nmo = SIZE(eigenvalues)
326 check_eigenvectors = diag_check_requested()
327
328 IF (diag_type == fm_diag_type_cusolver .AND. direct_generalized_diagonalization .AND. &
329 nao >= cusolver_n_min) THEN
330 ! Use cuSolverMP generalized eigenvalue solver without a CP2K-side
331 ! Cholesky reduction.
332 IF (check_eigenvectors) THEN
333 CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
334 CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
335 CALL cp_cfm_to_cfm(bmatrix, overlap_check)
336 END IF
337 CALL cp_cfm_general_cusolver(amatrix, bmatrix, work, evals)
338 IF (check_eigenvectors) THEN
339 CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
340 CALL cp_cfm_release(scratch_check)
341 CALL cp_cfm_release(overlap_check)
342 END IF
343#if defined(__DLAF)
344 ELSE IF (diag_type == fm_diag_type_dlaf .AND. direct_generalized_diagonalization .AND. &
345 nao >= dlaf_neigvec_min) THEN
346 ! Initialize DLA-Future on-demand; if already initialized, does nothing
347 CALL cp_dlaf_initialize()
348
349 ! Create DLAF grid from BLACS context; if already present, does nothing
350 CALL cp_dlaf_create_grid(amatrix%matrix_struct%context%get_handle())
351 CALL cp_dlaf_create_grid(bmatrix%matrix_struct%context%get_handle())
352 CALL cp_dlaf_create_grid(eigenvectors%matrix_struct%context%get_handle())
353
354 ! Use DLA-Future generalized eigenvalue solver for large matrices
355 IF (check_eigenvectors) THEN
356 CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
357 CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
358 CALL cp_cfm_to_cfm(bmatrix, overlap_check)
359 END IF
360 CALL cp_cfm_diag_gen_dlaf(amatrix, bmatrix, work, evals)
361 IF (check_eigenvectors) THEN
362 CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
363 CALL cp_cfm_release(scratch_check)
364 CALL cp_cfm_release(overlap_check)
365 END IF
366#endif
367#if defined(__parallel)
368 ELSE IF (diag_type == fm_diag_type_scalapack .AND. direct_generalized_diagonalization) THEN
369 ! Use ScaLAPACK generalized eigenvalue solver without a CP2K-side
370 ! Cholesky reduction.
371 IF (check_eigenvectors) THEN
372 CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
373 CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
374 CALL cp_cfm_to_cfm(bmatrix, overlap_check)
375 END IF
376 CALL cp_cfm_geeig_scalapack(amatrix, bmatrix, work, evals)
377 IF (check_eigenvectors) THEN
378 CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
379 CALL cp_cfm_release(scratch_check)
380 CALL cp_cfm_release(overlap_check)
381 END IF
382#endif
383 ELSE
384 ! Cholesky decompose S=U(T)U
385 CALL cp_cfm_cholesky_decompose(bmatrix)
386 ! Invert to get U^(-1)
387 CALL cp_cfm_triangular_invert(bmatrix)
388 ! Reduce to get U^(-T) * H * U^(-1)
389 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
390 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
391 ! Diagonalize
392 CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
393 ! Restore vectors C = U^(-1) * C*
394 CALL cp_cfm_triangular_multiply(bmatrix, work)
395 END IF
396
397 CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
398 eigenvalues(1:nmo) = evals(1:nmo)
399
400 DEALLOCATE (evals)
401
402 CALL timestop(handle)
403
404 END SUBROUTINE cp_cfm_geeig
405
406! **************************************************************************************************
407!> \brief General Eigenvalue Problem AX = BXE using ScaLAPACK PZHEGVX.
408!> \param amatrix ...
409!> \param bmatrix ...
410!> \param eigenvectors ...
411!> \param eigenvalues ...
412! **************************************************************************************************
413 SUBROUTINE cp_cfm_geeig_scalapack(amatrix, bmatrix, eigenvectors, eigenvalues)
414
415 TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
416 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
417
418 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_geeig_scalapack'
419
420#if defined(__parallel)
421 REAL(kind=dp), PARAMETER :: orfac = -1.0_dp, &
422 vl = 0.0_dp, &
423 vu = 0.0_dp
424
425 COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
426 COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a, b, z
427 INTEGER :: handle, info, liwork, lwork, lrwork, &
428 m, n, nb, neig, npcol, nprow, nz
429 INTEGER, DIMENSION(9) :: desca, descb, descz
430 INTEGER, DIMENSION(:), ALLOCATABLE :: iclustr, ifail, iwork
431 REAL(kind=dp) :: abstol
432 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: gap, rwork, w
433
434 INTEGER :: mq0, nn, np0, npe
435 INTEGER, EXTERNAL :: iceil, numroc
436 REAL(kind=dp), EXTERNAL :: dlamch
437#if defined (__HAS_IEEE_EXCEPTIONS)
438 LOGICAL, DIMENSION(5) :: halt
439#endif
440#else
441 INTEGER :: handle
442#endif
443
444 CALL timeset(routinen, handle)
445
446#if defined(__parallel)
447 n = amatrix%matrix_struct%nrow_global
448 neig = min(SIZE(eigenvalues), n)
449
450 IF (neig == 0) THEN
451 CALL timestop(handle)
452 RETURN
453 END IF
454
455 IF (amatrix%matrix_struct%nrow_block /= amatrix%matrix_struct%ncol_block) THEN
456 cpabort("ERROR in "//routinen//": Invalid blocksize (no square blocks) found")
457 END IF
458
459 a => amatrix%local_data
460 b => bmatrix%local_data
461 z => eigenvectors%local_data
462 desca(:) = amatrix%matrix_struct%descriptor(:)
463 descb(:) = bmatrix%matrix_struct%descriptor(:)
464 descz(:) = eigenvectors%matrix_struct%descriptor(:)
465
466 nprow = amatrix%matrix_struct%context%num_pe(1)
467 npcol = amatrix%matrix_struct%context%num_pe(2)
468 npe = nprow*npcol
469 nb = amatrix%matrix_struct%nrow_block
470 nn = max(n, nb, 2)
471 np0 = numroc(nn, nb, 0, 0, nprow)
472 mq0 = max(numroc(nn, nb, 0, 0, npcol), nb)
473
474 lwork = n + (np0 + mq0 + nb)*nb
475 lrwork = 4*n + max(5*nn, np0*mq0) + iceil(neig, npe)*nn + max(0, neig - 1)*n
476 liwork = 6*max(n, npe + 1, 4)
477
478 ALLOCATE (gap(npe))
479 gap = 0.0_dp
480 ALLOCATE (iclustr(2*npe))
481 iclustr = 0
482 ALLOCATE (ifail(n))
483 ifail = 0
484 ALLOCATE (iwork(liwork))
485 ALLOCATE (rwork(lrwork))
486 ALLOCATE (w(n))
487 ALLOCATE (work(lwork))
488
489 abstol = 2.0_dp*dlamch("S")
490
491#if defined (__HAS_IEEE_EXCEPTIONS)
492 CALL ieee_get_halting_mode(ieee_all, halt)
493 CALL ieee_set_halting_mode(ieee_all, .false.)
494#endif
495 CALL pzhegvx(1, "V", "I", "U", n, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb, &
496 vl, vu, 1, neig, abstol, m, nz, w(1), orfac, z(1, 1), 1, 1, descz, &
497 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, ifail(1), &
498 iclustr(1), gap(1), info)
499#if defined (__HAS_IEEE_EXCEPTIONS)
500 CALL ieee_set_halting_mode(ieee_all, halt)
501#endif
502
503 IF (info /= 0 .OR. m < neig .OR. nz < neig) THEN
504 cpabort("ERROR in PZHEGVX (ScaLAPACK), info="//trim(cp_to_string(info)))
505 END IF
506
507 eigenvalues(:) = 0.0_dp
508 eigenvalues(1:neig) = w(1:neig)
509
510 DEALLOCATE (gap, iclustr, ifail, iwork, rwork, w, work)
511#else
512 mark_used(amatrix)
513 mark_used(bmatrix)
514 mark_used(eigenvectors)
515 mark_used(eigenvalues)
516 cpabort("ERROR in "//routinen//": PZHEGVX requested without ScaLAPACK support")
517#endif
518
519 CALL timestop(handle)
520
521 END SUBROUTINE cp_cfm_geeig_scalapack
522
523! **************************************************************************************************
524!> \brief General Eigenvalue Problem AX = BXE
525!> Use canonical orthogonalization
526!> \param amatrix ...
527!> \param bmatrix ...
528!> \param eigenvectors ...
529!> \param eigenvalues ...
530!> \param work ...
531!> \param epseig ...
532! **************************************************************************************************
533 SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
534
535 TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
536 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
537 TYPE(cp_cfm_type), INTENT(IN) :: work
538 REAL(kind=dp), INTENT(IN) :: epseig
539
540 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_geeig_canon'
541
542 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
543 INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
544 nmo, nx
545 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
546
547 CALL timeset(routinen, handle)
548
549 ! Test sizes
550 CALL cp_cfm_get_info(amatrix, nrow_global=nao)
551 nmo = SIZE(eigenvalues)
552 ALLOCATE (evals(nao), cevals(nao))
553
554 ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
555 CALL cp_cfm_scale(-z_one, bmatrix)
556 CALL cp_cfm_heevd(bmatrix, work, evals)
557 evals(:) = -evals(:)
558 nc = nao
559 DO i = 1, nao
560 IF (evals(i) < epseig) THEN
561 nc = i - 1
562 EXIT
563 END IF
564 END DO
565 cpassert(nc /= 0)
566
567 IF (nc /= nao) THEN
568 IF (nc < nmo) THEN
569 ! Copy NULL space definition to last vectors of eigenvectors (if needed)
570 ncol = nmo - nc
571 CALL cp_cfm_to_cfm(work, eigenvectors, ncol, nc + 1, nc + 1)
572 END IF
573 ! Set NULL space in eigenvector matrix of S to zero
574 DO icol = nc + 1, nao
575 DO irow = 1, nao
576 CALL cp_cfm_set_element(work, irow, icol, z_zero)
577 END DO
578 END DO
579 ! Set small eigenvalues to a dummy save value
580 evals(nc + 1:nao) = 1.0_dp
581 END IF
582 ! Calculate U*s**(-1/2)
583 cevals(:) = cmplx(1.0_dp/sqrt(evals(:)), 0.0_dp, kind=dp)
584 CALL cp_cfm_column_scale(work, cevals)
585 ! Reduce to get U^(-C) * H * U^(-1)
586 CALL cp_cfm_gemm("C", "N", nao, nao, nao, z_one, work, amatrix, z_zero, bmatrix)
587 CALL cp_cfm_gemm("N", "N", nao, nao, nao, z_one, bmatrix, work, z_zero, amatrix)
588 IF (nc /= nao) THEN
589 ! set diagonal values to save large value
590 DO icol = nc + 1, nao
591 CALL cp_cfm_set_element(amatrix, icol, icol, cmplx(10000.0_dp, 0.0_dp, kind=dp))
592 END DO
593 END IF
594 ! Diagonalize
595 CALL cp_cfm_heevd(amatrix, bmatrix, evals)
596 eigenvalues(1:nmo) = evals(1:nmo)
597 nx = min(nc, nmo)
598 ! Restore vectors C = U^(-1) * C*
599 CALL cp_cfm_gemm("N", "N", nao, nx, nc, z_one, work, bmatrix, z_zero, eigenvectors)
600
601 DEALLOCATE (evals)
602
603 CALL timestop(handle)
604
605 END SUBROUTINE cp_cfm_geeig_canon
606
607END MODULE cp_cfm_diag
cp_cfm_basic_linalg::cp_cfm_scale
Definition cp_cfm_basic_linalg.F:60
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:56
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
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:563
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:980
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
Definition cp_cfm_basic_linalg.F:633
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:1062
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
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_cfm_diag.F:309
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:74
cp_cfm_diag::cp_cfm_geeig_canon
subroutine, public cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
General Eigenvalue Problem AX = BXE Use canonical orthogonalization.
Definition cp_cfm_diag.F:534
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_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:182
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:302
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:125
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:630
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:106
cp_fm_diag::fm_diag_type_dlaf
integer, parameter, public fm_diag_type_dlaf
Definition cp_fm_diag.F:106
cp_fm_diag::fm_diag_type_scalapack
integer, parameter, public fm_diag_type_scalapack
Definition cp_fm_diag.F:106
cp_fm_diag::direct_generalized_diagonalization
logical, save, public direct_generalized_diagonalization
Definition cp_fm_diag.F:100
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:314
cp_fm_diag::diag_type
integer, save, public diag_type
Definition cp_fm_diag.F:88
cp_fm_diag::cusolver_n_min
integer, parameter, public cusolver_n_min
Definition cp_fm_diag.F:94
cp_fm_diag::diag_check_requested
logical function, public diag_check_requested()
Return whether diagonalization checks should be performed.
Definition cp_fm_diag.F:299
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
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:46
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:143
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
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_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:69