(git:f89ef83)
Loading...
Searching...
No Matches
cp_fm_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 some of the diagonalization schemes available for
10!> cp_fm_type. cp_fm_power also moved here as it is very related
11!> \note
12!> first version : most routines imported
13!> \par History
14!> - unused Jacobi routines removed, cosmetics (05.04.06,MK)
15!> \author Joost VandeVondele (2003-08)
16! **************************************************************************************************
17MODULE cp_fm_diag
18 USE cp_blacs_types, ONLY: cp_blacs_type
19 USE cp_blacs_env, ONLY: cp_blacs_env_type
20 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale, &
21 cp_fm_gemm, &
22 cp_fm_scale, &
23 cp_fm_syrk, &
24 cp_fm_triangular_invert, &
25 cp_fm_triangular_multiply, &
26 cp_fm_uplo_to_full
27 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
28 USE cp_fm_diag_utils, ONLY: cp_fm_redistribute_end, &
29 cp_fm_redistribute_start
30 USE cp_fm_elpa, ONLY: cp_fm_diag_elpa, &
31 finalize_elpa_library, &
32 initialize_elpa_library, &
33 set_elpa_kernel
34 USE cp_cfm_elpa, ONLY: set_elpa_c_kernel
35 USE cp_fm_cusolver_api, ONLY: cp_fm_diag_cusolver, &
36 cp_fm_general_cusolver
37#if defined(__DLAF)
38 USE cp_fm_dlaf_api, ONLY: cp_fm_diag_dlaf, cp_fm_diag_gen_dlaf
39 USE cp_dlaf_utils_api, ONLY: cp_dlaf_initialize, cp_dlaf_finalize
40#endif
41 USE cp_fm_types, ONLY: cp_fm_get_info, &
42 cp_fm_set_element, &
43 cp_fm_to_fm, &
44 cp_fm_type, &
45 cp_fm_create, &
46 cp_fm_get_info, &
47 cp_fm_release, &
48 cp_fm_set_all, &
49 cp_fm_to_fm, &
50 cp_fm_to_fm_submat, &
51 cp_fm_type
52 USE cp_fm_struct, ONLY: cp_fm_struct_equivalent, &
53 cp_fm_struct_create, &
54 cp_fm_struct_release, &
55 cp_fm_struct_type
56 USE cp_log_handling, ONLY: cp_logger_get_default_unit_nr, &
57 cp_get_default_logger, &
58 cp_logger_get_default_io_unit, &
59 cp_logger_type, &
60 cp_to_string
61 USE cp_log_handling, ONLY: cp_get_default_logger, &
62 cp_logger_get_default_unit_nr, &
63 cp_logger_get_unit_nr, &
64 cp_logger_type
65 USE kinds, ONLY: default_string_length, &
66 dp
67 USE machine, ONLY: default_output_unit, &
68 m_memory
69 USE parallel_gemm_api, ONLY: parallel_gemm
70#if defined (__parallel)
71 USE message_passing, ONLY: mp_comm_type
72#endif
73#if defined (__HAS_IEEE_EXCEPTIONS)
74 USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
75 ieee_set_halting_mode, &
76 ieee_all
77#endif
78#include "../base/base_uses.f90"
79
80 IMPLICIT NONE
81 PRIVATE
82
83 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_diag'
84
85 REAL(kind=dp), PARAMETER, PUBLIC :: eps_check_diag_default = 5.0e-14_dp
86
87 ! Placeholder eigenvalue for linearly dependent modes removed by canonical orthogonalization.
88 ! Pushes unphysical states above the real spectrum so downstream routines can ignore them.
89 REAL(kind=dp), PARAMETER, PUBLIC :: set_removed_eigval_to = 10000.0_dp
90
91 ! The following saved variables are diagonalization global
92 ! Stores the default library for diagonalization
93 INTEGER, SAVE, PUBLIC :: diag_type = 0
94 ! Whether the diagonalization library was requested explicitly (e.g. via PREFERRED_DIAG_LIBRARY);
95 ! complex matrices use ELPA only in that case.
96 LOGICAL, SAVE, PUBLIC :: diag_lib_explicit = .false.
97 ! Minimum number of eigenvectors for the use of the ELPA eigensolver.
98 ! The ScaLAPACK eigensolver is used as fallback for all smaller cases.
99 INTEGER, SAVE, PUBLIC :: elpa_neigvec_min = 0
100 ! Minimum matrix size for the use of the cuSOLVERMp eigensolver.
101 ! Smaller matrices use the ScaLAPACK fallback to avoid GPU launch overheads.
102 INTEGER, PARAMETER, PUBLIC :: cusolver_n_min = 64
103 ! Minimum number of eigenvectors for the use of the DLAF eigensolver.
104 ! The ScaLAPACK eigensolver is used as fallback for all smaller cases.
105 INTEGER, SAVE, PUBLIC :: dlaf_neigvec_min = 0
106 LOGICAL, SAVE, PUBLIC :: direct_generalized_diagonalization = .false.
107 ! Threshold value for the orthonormality check of the eigenvectors obtained
108 ! after a diagonalization. A negative value disables the check.
109 REAL(kind=dp), SAVE :: eps_check_diag = -1.0_dp
110
111 ! Constants for the diag_type above
112 INTEGER, PARAMETER, PUBLIC :: fm_diag_type_scalapack = 101, &
113 fm_diag_type_elpa = 102, &
114 fm_diag_type_cusolver = 103, &
115 fm_diag_type_dlaf = 104
116#if defined(__CUSOLVERMP)
117 INTEGER, PARAMETER, PUBLIC :: fm_diag_type_default = fm_diag_type_cusolver
118#elif defined(__ELPA)
119 INTEGER, PARAMETER, PUBLIC :: fm_diag_type_default = fm_diag_type_elpa
120#else
121 INTEGER, PARAMETER, PUBLIC :: fm_diag_type_default = fm_diag_type_scalapack
122#endif
123
124 ! Public subroutines
125 PUBLIC :: choose_eigv_solver, &
126 cp_fm_block_jacobi, &
127 cp_fm_power, &
128 cp_fm_syevd, &
129 cp_fm_syevx, &
130 cp_fm_svd, &
131 cp_fm_geeig, &
132 cp_fm_geeig_canon, &
133 diag_check_requested, &
134 diag_check_warning_threshold, &
135 diag_init, &
136 diag_finalize
137
138CONTAINS
139
140! **************************************************************************************************
141!> \brief Setup the diagonalization library to be used
142!> \param diag_lib diag_library flag from GLOBAL section in input
143!> \param fallback_applied .TRUE. if support for the requested library was not compiled-in and fallback
144!> to ScaLAPACK was applied, .FALSE. otherwise.
145!> \param elpa_kernel integer that determines which ELPA kernel to use for diagonalization
146!> \param elpa_c_kernel ...
147!> \param elpa_neigvec_min_input ...
148!> \param elpa_qr logical that determines if ELPA should try to use QR to accelerate the
149!> diagonalization procedure of suitably sized matrices
150!> \param elpa_print logical that determines if information about the ELPA diagonalization should
151!> be printed
152!> \param elpa_one_stage logical that enables the one-stage solver
153!> \param dlaf_neigvec_min_input ...
154!> \param eps_check_diag_input ...
155!> \param direct_generalized_diagonalization_input ...
156!> \param diag_lib_explicit_input ...
157!> \par History
158!> - Add support for DLA-Future (05.09.2023, RMeli)
159!> \author MI 11.2013
160! **************************************************************************************************
161 SUBROUTINE diag_init(diag_lib, fallback_applied, elpa_kernel, elpa_c_kernel, elpa_neigvec_min_input, &
162 elpa_qr, elpa_print, elpa_one_stage, dlaf_neigvec_min_input, eps_check_diag_input, &
163 direct_generalized_diagonalization_input, diag_lib_explicit_input)
164 CHARACTER(LEN=*), INTENT(IN) :: diag_lib
165 LOGICAL, INTENT(OUT) :: fallback_applied
166 INTEGER, INTENT(IN) :: elpa_kernel
167 INTEGER, INTENT(IN), OPTIONAL :: elpa_c_kernel
168 INTEGER, INTENT(IN) :: elpa_neigvec_min_input
169 LOGICAL, INTENT(IN) :: elpa_qr, elpa_print, elpa_one_stage
170 INTEGER, INTENT(IN) :: dlaf_neigvec_min_input
171 REAL(kind=dp), INTENT(IN) :: eps_check_diag_input
172 LOGICAL, INTENT(IN), OPTIONAL :: direct_generalized_diagonalization_input, &
173 diag_lib_explicit_input
174
175 LOGICAL, SAVE :: initialized = .false.
176
177 fallback_applied = .false.
178
179 IF (diag_lib == "ScaLAPACK") THEN
180 diag_type = fm_diag_type_scalapack
181 ELSE IF (diag_lib == "ELPA") THEN
182#if defined (__ELPA)
183 ! ELPA is requested and available
184 diag_type = fm_diag_type_elpa
185#else
186 ! ELPA library requested but not linked, switch back to SL
187 diag_type = fm_diag_type_scalapack
188 fallback_applied = .true.
189#endif
190 ELSE IF (diag_lib == "cuSOLVER") THEN
191 diag_type = fm_diag_type_cusolver
192 ELSE IF (diag_lib == "DLAF") THEN
193#if defined (__DLAF)
194 diag_type = fm_diag_type_dlaf
195#else
196 cpabort("ERROR in diag_init: CP2K was not compiled with DLA-Future support")
197#endif
198 ELSE
199 cpabort("ERROR in diag_init: Initialization of unknown diagonalization library requested")
200 END IF
201
202 ! Complex matrices use ELPA only when the diagonalization library was
203 ! requested explicitly (e.g. via PREFERRED_DIAG_LIBRARY ELPA).
204 diag_lib_explicit = .true.
205 IF (PRESENT(diag_lib_explicit_input)) diag_lib_explicit = diag_lib_explicit_input
206
207 ! Initialization of requested diagonalization library
208 IF (.NOT. initialized .AND. diag_type == fm_diag_type_elpa) THEN
209 CALL initialize_elpa_library(one_stage=elpa_one_stage, qr=elpa_qr, should_print=elpa_print)
210 CALL set_elpa_kernel(elpa_kernel)
211 IF (PRESENT(elpa_c_kernel)) CALL set_elpa_c_kernel(elpa_c_kernel)
212 initialized = .true.
213 END IF
214#if defined(__DLAF)
215 IF (.NOT. initialized .AND. diag_type == fm_diag_type_dlaf) THEN
216 CALL cp_dlaf_initialize()
217 initialized = .true.
218 END IF
219 dlaf_neigvec_min = dlaf_neigvec_min_input
220#else
221 mark_used(dlaf_neigvec_min_input)
222#endif
223
224 elpa_neigvec_min = elpa_neigvec_min_input
225 eps_check_diag = eps_check_diag_input
226 IF (PRESENT(direct_generalized_diagonalization_input)) THEN
227 direct_generalized_diagonalization = direct_generalized_diagonalization_input
228 ELSE
229 direct_generalized_diagonalization = .false.
230 END IF
231
232 END SUBROUTINE diag_init
233
234! **************************************************************************************************
235!> \brief Finalize the diagonalization library
236! **************************************************************************************************
237 SUBROUTINE diag_finalize()
238#if defined (__ELPA)
239 IF (diag_type == fm_diag_type_elpa) &
240 CALL finalize_elpa_library()
241#endif
242#if defined (__DLAF)
243 IF (diag_type == fm_diag_type_dlaf) &
244 CALL cp_dlaf_finalize()
245#endif
246 END SUBROUTINE diag_finalize
247
248! **************************************************************************************************
249!> \brief Choose the Eigensolver depending on which library is available
250!> ELPA seems to be unstable for small systems
251!> \param matrix ...
252!> \param eigenvectors ...
253!> \param eigenvalues ...
254!> \param info ...
255!> \par info If present returns error code and prevents program stops.
256!> Works currently only for cp_fm_syevd with ScaLAPACK.
257!> Other solvers will end the program regardless of PRESENT(info).
258!> \par History
259!> - Do not use ELPA for small matrices and use instead ScaLAPACK as fallback (10.05.2021, MK)
260! **************************************************************************************************
261 SUBROUTINE choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
262
263 TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
264 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
265 INTEGER, INTENT(OUT), OPTIONAL :: info
266
267 CHARACTER(LEN=*), PARAMETER :: routinen = 'choose_eigv_solver'
268
269 ! Sample peak memory
270 CALL m_memory()
271
272 IF (PRESENT(info)) info = 0 ! Default for solvers that do not return an info.
273
274 IF (diag_type == fm_diag_type_scalapack) THEN
275 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
276
277 ELSE IF (diag_type == fm_diag_type_elpa) THEN
278 IF (matrix%matrix_struct%nrow_global < elpa_neigvec_min) THEN
279 ! We don't trust ELPA with very small matrices.
280 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
281 ELSE
282 CALL cp_fm_diag_elpa(matrix, eigenvectors, eigenvalues)
283 END IF
284
285 ELSE IF (diag_type == fm_diag_type_cusolver) THEN
286 IF (matrix%matrix_struct%nrow_global < cusolver_n_min) THEN
287 ! We don't trust cuSolver with very small matrices.
288 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
289 ELSE
290 CALL cp_fm_diag_cusolver(matrix, eigenvectors, eigenvalues)
291 END IF
292
293#if defined(__DLAF)
294 ELSE IF (diag_type == fm_diag_type_dlaf) THEN
295 IF (matrix%matrix_struct%nrow_global < dlaf_neigvec_min) THEN
296 ! Use ScaLAPACK for small matrices
297 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
298 ELSE
299 CALL cp_fm_diag_dlaf(matrix, eigenvectors, eigenvalues)
300 END IF
301#endif
302
303 ELSE
304 cpabort("ERROR in "//routinen//": Invalid diagonalization type requested")
305 END IF
306
307 CALL check_diag(matrix, eigenvectors, nvec=SIZE(eigenvalues))
308
309 END SUBROUTINE choose_eigv_solver
310
311! **************************************************************************************************
312!> \brief Return whether diagonalization checks should be performed.
313!> \return ...
314! **************************************************************************************************
315 FUNCTION diag_check_requested() RESULT(check_requested)
316 LOGICAL :: check_requested
317
318#if defined(__CHECK_DIAG)
319 check_requested = .true.
320#else
321 check_requested = eps_check_diag >= 0.0_dp
322#endif
323
324 END FUNCTION diag_check_requested
325
326! **************************************************************************************************
327!> \brief Return the warning threshold for diagonalization checks.
328!> \return ...
329! **************************************************************************************************
330 FUNCTION diag_check_warning_threshold() RESULT(eps_warning)
331 REAL(kind=dp) :: eps_warning
332
333 eps_warning = eps_check_diag_default
334 IF (eps_check_diag >= 0.0_dp) THEN
335 eps_warning = eps_check_diag
336 END IF
337
338 END FUNCTION diag_check_warning_threshold
339
340! **************************************************************************************************
341!> \brief Check result of diagonalization, i.e. the orthonormality of the eigenvectors
342!> \param matrix Work matrix
343!> \param eigenvectors Eigenvectors to be checked
344!> \param nvec ...
345! **************************************************************************************************
346 SUBROUTINE check_diag(matrix, eigenvectors, nvec)
347
348 TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
349 INTEGER, INTENT(IN) :: nvec
350
351 CHARACTER(LEN=*), PARAMETER :: routinen = 'check_diag'
352
353 CHARACTER(LEN=default_string_length) :: diag_type_name
354 REAL(kind=dp) :: eps, eps_abort, eps_warning, gold, test
355 INTEGER :: handle, i, j, ncol, nrow, output_unit
356 LOGICAL :: check_eigenvectors
357#if defined(__parallel)
358 TYPE(cp_blacs_env_type), POINTER :: context
359 INTEGER :: il, jl, ipcol, iprow, &
360 mypcol, myprow, npcol, nprow
361 INTEGER, DIMENSION(9) :: desca
362#endif
363
364 CALL timeset(routinen, handle)
365
366 output_unit = default_output_unit
367 check_eigenvectors = diag_check_requested()
368 eps_warning = diag_check_warning_threshold()
369 eps_abort = 10.0_dp*eps_warning
370
371 gold = 0.0_dp
372 test = 0.0_dp
373 eps = 0.0_dp
374
375 IF (check_eigenvectors) THEN
376#if defined(__parallel)
377 nrow = eigenvectors%matrix_struct%nrow_global
378 ncol = min(eigenvectors%matrix_struct%ncol_global, nvec)
379 CALL cp_fm_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, eigenvectors, eigenvectors, 0.0_dp, matrix)
380 context => matrix%matrix_struct%context
381 myprow = context%mepos(1)
382 mypcol = context%mepos(2)
383 nprow = context%num_pe(1)
384 npcol = context%num_pe(2)
385 desca(:) = matrix%matrix_struct%descriptor(:)
386 outer: DO j = 1, ncol
387 DO i = 1, ncol
388 CALL infog2l(i, j, desca, nprow, npcol, myprow, mypcol, il, jl, iprow, ipcol)
389 IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
390 gold = merge(0.0_dp, 1.0_dp, i /= j)
391 test = matrix%local_data(il, jl)
392 eps = abs(test - gold)
393 IF (eps > eps_warning) EXIT outer
394 END IF
395 END DO
396 END DO outer
397#else
398 nrow = SIZE(eigenvectors%local_data, 1)
399 ncol = min(SIZE(eigenvectors%local_data, 2), nvec)
400 CALL dgemm("T", "N", ncol, ncol, nrow, 1.0_dp, &
401 eigenvectors%local_data(1, 1), nrow, &
402 eigenvectors%local_data(1, 1), nrow, &
403 0.0_dp, matrix%local_data(1, 1), nrow)
404 outer: DO j = 1, ncol
405 DO i = 1, ncol
406 gold = merge(0.0_dp, 1.0_dp, i /= j)
407 test = matrix%local_data(i, j)
408 eps = abs(test - gold)
409 IF (eps > eps_warning) EXIT outer
410 END DO
411 END DO outer
412#endif
413 IF (eps > eps_warning) THEN
414 IF (diag_type == fm_diag_type_scalapack) THEN
415 diag_type_name = "SYEVD"
416 ELSE IF (diag_type == fm_diag_type_elpa) THEN
417 diag_type_name = "ELPA"
418 ELSE IF (diag_type == fm_diag_type_cusolver) THEN
419 diag_type_name = "CUSOLVER"
420 ELSE IF (diag_type == fm_diag_type_dlaf) THEN
421 diag_type_name = "DLAF"
422 ELSE
423 cpabort("Unknown diag_type")
424 END IF
425 WRITE (unit=output_unit, fmt="(/,T2,A,/,T2,A,I0,A,I0,A,F0.15,/,T2,A,F0.0,A,ES10.3)") &
426 "The eigenvectors returned by "//trim(diag_type_name)//" are not orthonormal", &
427 "Matrix element (", i, ", ", j, ") = ", test, &
428 "The deviation from the expected value ", gold, " is", eps
429 IF (eps > eps_abort) THEN
430 cpabort("ERROR in "//routinen//": Check of matrix diagonalization failed")
431 ELSE
432 cpwarn("Check of matrix diagonalization failed in routine "//routinen)
433 END IF
434 END IF
435 END IF
436
437 CALL timestop(handle)
438
439 END SUBROUTINE check_diag
440
441! **************************************************************************************************
442!> \brief Check C^T*S*C = I for a generalized eigenvalue problem.
443!> \param overlap original overlap matrix S; used as work matrix and overwritten
444!> \param eigenvectors eigenvectors C to be checked
445!> \param scratch work matrix
446!> \param nvec ...
447! **************************************************************************************************
448 SUBROUTINE check_generalized_diag(overlap, eigenvectors, scratch, nvec)
449
450 TYPE(cp_fm_type), INTENT(IN) :: eigenvectors
451 TYPE(cp_fm_type), INTENT(INOUT) :: overlap, scratch
452 INTEGER, INTENT(IN) :: nvec
453
454 CHARACTER(LEN=*), PARAMETER :: routinen = 'check_generalized_diag'
455
456 CHARACTER(LEN=default_string_length) :: diag_type_name
457 REAL(kind=dp) :: eps, eps_abort, eps_warning, gold, test
458 INTEGER :: handle, i, j, ncol, nrow, output_unit
459#if defined(__parallel)
460 TYPE(cp_blacs_env_type), POINTER :: context
461 INTEGER :: il, jl, ipcol, iprow, &
462 mypcol, myprow, npcol, nprow
463 INTEGER, DIMENSION(9) :: desca
464#endif
465
466 CALL timeset(routinen, handle)
467
468 IF (.NOT. diag_check_requested()) THEN
469 CALL timestop(handle)
470 RETURN
471 END IF
472
473 output_unit = default_output_unit
474 eps_warning = diag_check_warning_threshold()
475 eps_abort = 10.0_dp*eps_warning
476
477 nrow = eigenvectors%matrix_struct%nrow_global
478 ncol = min(eigenvectors%matrix_struct%ncol_global, nvec)
479
480 CALL parallel_gemm("N", "N", nrow, ncol, nrow, 1.0_dp, overlap, eigenvectors, 0.0_dp, scratch)
481 CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, eigenvectors, scratch, 0.0_dp, overlap)
482
483 gold = 0.0_dp
484 test = 0.0_dp
485 eps = 0.0_dp
486
487#if defined(__parallel)
488 context => overlap%matrix_struct%context
489 myprow = context%mepos(1)
490 mypcol = context%mepos(2)
491 nprow = context%num_pe(1)
492 npcol = context%num_pe(2)
493 desca(:) = overlap%matrix_struct%descriptor(:)
494 outer: DO j = 1, ncol
495 DO i = 1, ncol
496 CALL infog2l(i, j, desca, nprow, npcol, myprow, mypcol, il, jl, iprow, ipcol)
497 IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
498 gold = merge(0.0_dp, 1.0_dp, i /= j)
499 test = overlap%local_data(il, jl)
500 eps = abs(test - gold)
501 IF (eps > eps_warning) EXIT outer
502 END IF
503 END DO
504 END DO outer
505#else
506 outer: DO j = 1, ncol
507 DO i = 1, ncol
508 gold = merge(0.0_dp, 1.0_dp, i /= j)
509 test = overlap%local_data(i, j)
510 eps = abs(test - gold)
511 IF (eps > eps_warning) EXIT outer
512 END DO
513 END DO outer
514#endif
515
516 IF (eps > eps_warning) THEN
517 IF (diag_type == fm_diag_type_scalapack) THEN
518 diag_type_name = "SYGVX"
519 ELSE IF (diag_type == fm_diag_type_elpa) THEN
520 diag_type_name = "ELPA"
521 ELSE IF (diag_type == fm_diag_type_cusolver) THEN
522 diag_type_name = "CUSOLVER"
523 ELSE IF (diag_type == fm_diag_type_dlaf) THEN
524 diag_type_name = "DLAF"
525 ELSE
526 cpabort("Unknown diag_type")
527 END IF
528 WRITE (unit=output_unit, fmt="(/,T2,A,/,T2,A,I0,A,I0,A,F0.15,/,T2,A,F0.0,A,ES10.3)") &
529 "The generalized eigenvectors returned by "//trim(diag_type_name)//" are not S-orthonormal", &
530 "Matrix element (", i, ", ", j, ") = ", test, &
531 "The deviation from the expected value ", gold, " is", eps
532 IF (eps > eps_abort) THEN
533 cpabort("ERROR in "//routinen//": Check of generalized matrix diagonalization failed")
534 ELSE
535 cpwarn("Check of generalized matrix diagonalization failed in routine "//routinen)
536 END IF
537 END IF
538
539 CALL timestop(handle)
540
541 END SUBROUTINE check_generalized_diag
542
543! **************************************************************************************************
544!> \brief Issues an error messages and exits (optionally only warns).
545!> \param mesg message to be issued
546!> \param info error code (optional)
547!> \param warn only warn (optional)
548! **************************************************************************************************
549 SUBROUTINE cp_fm_error(mesg, info, warn)
550 CHARACTER(LEN=*), INTENT(IN) :: mesg
551 INTEGER, INTENT(IN), OPTIONAL :: info
552 LOGICAL, INTENT(IN), OPTIONAL :: warn
553
554 CHARACTER(LEN=2*default_string_length) :: message
555 LOGICAL :: warning
556
557 IF (PRESENT(info)) THEN
558 WRITE (message, "(A,A,I0,A)") mesg, " (INFO = ", info, ")"
559 ELSE
560 WRITE (message, "(A)") mesg
561 END IF
562
563 IF (PRESENT(warn)) THEN
564 warning = warn
565 ELSE ! abort
566 warning = .false.
567 END IF
568
569 IF (warning) THEN
570 cpwarn(trim(message))
571 ELSE
572 cpabort(trim(message))
573 END IF
574 END SUBROUTINE cp_fm_error
575
576! **************************************************************************************************
577!> \brief Computes all eigenvalues and vectors of a real symmetric matrix
578!> significantly faster than syevx, scales also much better.
579!> Needs workspace to allocate all the eigenvectors
580!> \param matrix ...
581!> \param eigenvectors ...
582!> \param eigenvalues ...
583!> \param info ...
584!> \par matrix is supposed to be in upper triangular form, and overwritten by this routine
585!> \par info If present returns error code and prevents program stops.
586!> Works currently only for scalapack.
587!> Other solvers will end the program regardless of PRESENT(info).
588! **************************************************************************************************
589 SUBROUTINE cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
590
591 TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
592 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
593 INTEGER, INTENT(OUT), OPTIONAL :: info
594
595 CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_fm_syevd'
596
597 INTEGER :: handle, myinfo, n, nmo
598 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eig
599#if defined(__parallel)
600 TYPE(cp_fm_type) :: eigenvectors_new, matrix_new
601#else
602 INTEGER :: liwork, lwork
603 INTEGER, DIMENSION(:), POINTER :: iwork
604 REAL(kind=dp), DIMENSION(:, :), POINTER :: m
605 REAL(kind=dp), DIMENSION(:), POINTER :: work
606 INTEGER, TARGET :: v(1)
607 REAL(kind=dp), TARGET :: w(1)
608#endif
609
610 CALL timeset(routinen, handle)
611
612 myinfo = 0
613
614 n = matrix%matrix_struct%nrow_global
615 ALLOCATE (eig(n))
616
617#if defined(__parallel)
618 ! Determine if the input matrix needs to be redistributed before diagonalization.
619 ! Heuristics are used to determine the optimal number of CPUs for diagonalization.
620 ! The redistributed matrix is stored in matrix_new, which is just a pointer
621 ! to the original matrix if no redistribution is required
622 CALL cp_fm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new)
623
624 ! Call scalapack on CPUs that hold the new matrix
625 IF (ASSOCIATED(matrix_new%matrix_struct)) THEN
626 IF (PRESENT(info)) THEN
627 CALL cp_fm_syevd_base(matrix_new, eigenvectors_new, eig, myinfo)
628 ELSE
629 CALL cp_fm_syevd_base(matrix_new, eigenvectors_new, eig)
630 END IF
631 END IF
632 ! Redistribute results and clean up
633 CALL cp_fm_redistribute_end(matrix, eigenvectors, eig, matrix_new, eigenvectors_new)
634#else
635 ! Retrieve the optimal work array sizes first
636 lwork = -1
637 liwork = -1
638 m => matrix%local_data
639 iwork => v
640 work => w
641
642 CALL dsyevd('V', 'U', n, m(1, 1), SIZE(m, 1), eig(1), work(1), lwork, iwork(1), liwork, myinfo)
643
644 IF (myinfo /= 0) THEN
645 CALL cp_fm_error("ERROR in DSYEVD: Work space query failed", myinfo, PRESENT(info))
646 END IF
647
648 ! Reallocate work arrays and perform diagonalisation
649 lwork = nint(work(1))
650 ALLOCATE (work(lwork))
651
652 liwork = iwork(1)
653 ALLOCATE (iwork(liwork))
654
655 CALL dsyevd('V', 'U', n, m(1, 1), SIZE(m, 1), eig(1), work(1), lwork, iwork(1), liwork, myinfo)
656
657 IF (myinfo /= 0) THEN
658 CALL cp_fm_error("ERROR in DSYEVD: Matrix diagonalization failed", myinfo, PRESENT(info))
659 END IF
660
661 CALL cp_fm_to_fm(matrix, eigenvectors)
662
663 DEALLOCATE (iwork)
664 DEALLOCATE (work)
665#endif
666
667 IF (PRESENT(info)) info = myinfo
668
669 nmo = SIZE(eigenvalues, 1)
670 IF (nmo > n) THEN
671 eigenvalues(1:n) = eig(1:n)
672 ELSE
673 eigenvalues(1:nmo) = eig(1:nmo)
674 END IF
675
676 DEALLOCATE (eig)
677
678 CALL check_diag(matrix, eigenvectors, n)
679
680 CALL timestop(handle)
681
682 END SUBROUTINE cp_fm_syevd
683
684! **************************************************************************************************
685!> \brief ...
686!> \param matrix ...
687!> \param eigenvectors ...
688!> \param eigenvalues ...
689!> \param info ...
690! **************************************************************************************************
691 SUBROUTINE cp_fm_syevd_base(matrix, eigenvectors, eigenvalues, info)
692
693 TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
694 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
695 INTEGER, INTENT(OUT), OPTIONAL :: info
696
697 CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_fm_syevd_base'
698
699 INTEGER :: handle, myinfo
700#if defined(__parallel)
701 TYPE(cp_blacs_env_type), POINTER :: context
702 INTEGER :: liwork, lwork, n
703 INTEGER, DIMENSION(9) :: descm, descv
704 INTEGER, DIMENSION(:), POINTER :: iwork
705 REAL(kind=dp), DIMENSION(:), POINTER :: work
706 REAL(kind=dp), DIMENSION(:, :), POINTER :: m, v
707 REAL(kind=dp), TARGET :: w(1)
708#if defined (__HAS_IEEE_EXCEPTIONS)
709 LOGICAL, DIMENSION(5) :: halt
710#endif
711#endif
712
713 CALL timeset(routinen, handle)
714
715 myinfo = 0
716
717#if defined(__parallel)
718
719 n = matrix%matrix_struct%nrow_global
720 m => matrix%local_data
721 context => matrix%matrix_struct%context
722 descm(:) = matrix%matrix_struct%descriptor(:)
723
724 v => eigenvectors%local_data
725 descv(:) = eigenvectors%matrix_struct%descriptor(:)
726
727 liwork = 7*n + 8*context%num_pe(2) + 2
728 ALLOCATE (iwork(liwork))
729
730 ! Work space query
731 lwork = -1
732 work => w
733
734 CALL pdsyevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
735 work(1), lwork, iwork(1), liwork, myinfo)
736
737 IF (matrix%matrix_struct%para_env%is_source() .AND. (myinfo /= 0)) THEN
738 CALL cp_fm_error("ERROR in PDSYEVD: Work space query failed", myinfo, PRESENT(info))
739 END IF
740
741 lwork = nint(work(1)) ! can be insufficient due to bug in reference ScaLAPACK
742#if !defined(__SCALAPACK_NO_WA)
743 ! Query workspace for QDORMTR as called by reference ScaLAPACK (PDSYEVD).
744 CALL pdormtr('L', 'U', 'N', n, n, m(1, 1), 1, 1, descm, m(1, 1), &
745 v(1, 1), 1, 1, descv, work(1), -1, myinfo)
746
747 IF (matrix%matrix_struct%para_env%is_source() .AND. (myinfo /= 0)) THEN
748 CALL cp_fm_error("ERROR in PDORMTR: Work space query failed", myinfo, PRESENT(info))
749 END IF
750
751 IF (lwork < (work(1) + 2*n)) THEN
752 lwork = nint(work(1)) + 2*n ! still wrong by 2*N
753 END IF
754#endif
755 ALLOCATE (work(lwork))
756
757 ! Initial/documented amount of liwork is exceeded (slightly worrisome too).
758 IF (liwork < iwork(1)) THEN
759 liwork = iwork(1)
760 DEALLOCATE (iwork)
761 ALLOCATE (iwork(liwork))
762 END IF
763
764 ! ScaLAPACK takes advantage of IEEE754 exceptions for speedup.
765 ! Therefore, we disable floating point traps temporarily.
766#if defined (__HAS_IEEE_EXCEPTIONS)
767 CALL ieee_get_halting_mode(ieee_all, halt)
768 CALL ieee_set_halting_mode(ieee_all, .false.)
769#endif
770
771 CALL pdsyevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
772 work(1), lwork, iwork(1), liwork, myinfo)
773
774#if defined (__HAS_IEEE_EXCEPTIONS)
775 CALL ieee_set_halting_mode(ieee_all, halt)
776#endif
777 IF (matrix%matrix_struct%para_env%is_source() .AND. (myinfo /= 0)) THEN
778 CALL cp_fm_error("ERROR in PDSYEVD: Matrix diagonalization failed", myinfo, PRESENT(info))
779 END IF
780
781 IF (PRESENT(info)) info = myinfo
782
783 DEALLOCATE (work)
784 DEALLOCATE (iwork)
785#else
786 mark_used(matrix)
787 mark_used(eigenvectors)
788 mark_used(eigenvalues)
789 myinfo = -1
790 IF (PRESENT(info)) info = myinfo
791 CALL cp_fm_error("ERROR in "//trim(routinen)// &
792 ": Matrix diagonalization using PDSYEVD requested without ScaLAPACK support")
793#endif
794
795 CALL timestop(handle)
796
797 END SUBROUTINE cp_fm_syevd_base
798
799! **************************************************************************************************
800!> \brief compute eigenvalues and optionally eigenvectors of a real symmetric matrix using scalapack.
801!> If eigenvectors are required this routine will replicate a full matrix on each CPU...
802!> if more than a handful of vectors are needed, use cp_fm_syevd instead
803!> \param matrix ...
804!> \param eigenvectors ...
805!> \param eigenvalues ...
806!> \param neig ...
807!> \param work_syevx ...
808!> \par matrix is supposed to be in upper triangular form, and overwritten by this routine
809!> neig is the number of vectors needed (default all)
810!> work_syevx evec calculation only, is the fraction of the working buffer allowed (1.0 use full buffer)
811!> reducing this saves time, but might cause the routine to fail
812! **************************************************************************************************
813 SUBROUTINE cp_fm_syevx(matrix, eigenvectors, eigenvalues, neig, work_syevx)
814
815 ! Diagonalise the symmetric n by n matrix using the LAPACK library.
816
817 TYPE(cp_fm_type), INTENT(IN) :: matrix
818 TYPE(cp_fm_type), OPTIONAL, INTENT(IN) :: eigenvectors
819 REAL(kind=dp), OPTIONAL, INTENT(IN) :: work_syevx
820 INTEGER, INTENT(IN), OPTIONAL :: neig
821 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
822
823 CHARACTER(LEN=*), PARAMETER :: routinen = "cp_fm_syevx"
824
825#if defined(__parallel)
826 REAL(kind=dp), PARAMETER :: orfac = -1.0_dp
827#endif
828 REAL(kind=dp), PARAMETER :: vl = 0.0_dp, &
829 vu = 0.0_dp
830
831 TYPE(cp_blacs_env_type), POINTER :: context
832 TYPE(cp_logger_type), POINTER :: logger
833 CHARACTER(LEN=1) :: job_type
834 REAL(kind=dp) :: abstol, work_syevx_local
835 INTEGER :: handle, info, liwork, lwork, &
836 m, n, nb, npcol, nprow, &
837 output_unit, neig_local
838 LOGICAL :: ionode, needs_evecs
839 INTEGER, DIMENSION(:), ALLOCATABLE :: ifail, iwork
840 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: w, work
841 REAL(kind=dp), DIMENSION(:, :), POINTER :: a, z
842
843 REAL(kind=dp), EXTERNAL :: dlamch
844
845#if defined(__parallel)
846 INTEGER :: nn, np0, npe, nq0, nz
847 INTEGER, DIMENSION(9) :: desca, descz
848 INTEGER, DIMENSION(:), ALLOCATABLE :: iclustr
849 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: gap
850 INTEGER, EXTERNAL :: iceil, numroc
851#else
852 INTEGER :: nla, nlz
853 INTEGER, EXTERNAL :: ilaenv
854#endif
855#if defined (__HAS_IEEE_EXCEPTIONS)
856 LOGICAL, DIMENSION(5) :: halt
857#endif
858
859 ! by default all
860 n = matrix%matrix_struct%nrow_global
861 neig_local = n
862 IF (PRESENT(neig)) neig_local = neig
863 IF (neig_local == 0) RETURN
864
865 CALL timeset(routinen, handle)
866
867 needs_evecs = PRESENT(eigenvectors)
868
869 logger => cp_get_default_logger()
870 ionode = logger%para_env%is_source()
871 n = matrix%matrix_struct%nrow_global
872
873 ! by default allocate all needed space
874 work_syevx_local = 1.0_dp
875 IF (PRESENT(work_syevx)) work_syevx_local = work_syevx
876
877 ! set scalapack job type
878 IF (needs_evecs) THEN
879 job_type = "V"
880 ELSE
881 job_type = "N"
882 END IF
883
884 ! target the most accurate calculation of the eigenvalues
885 abstol = 2.0_dp*dlamch("S")
886
887 context => matrix%matrix_struct%context
888 nprow = context%num_pe(1)
889 npcol = context%num_pe(2)
890
891 ALLOCATE (w(n))
892 eigenvalues(:) = 0.0_dp
893#if defined(__parallel)
894
895 IF (matrix%matrix_struct%nrow_block /= matrix%matrix_struct%ncol_block) THEN
896 cpabort("ERROR in "//routinen//": Invalid blocksize (no square blocks) found")
897 END IF
898
899 a => matrix%local_data
900 desca(:) = matrix%matrix_struct%descriptor(:)
901
902 IF (needs_evecs) THEN
903 z => eigenvectors%local_data
904 descz(:) = eigenvectors%matrix_struct%descriptor(:)
905 ELSE
906 ! z will not be referenced
907 z => matrix%local_data
908 descz = desca
909 END IF
910
911 ! Get the optimal work storage size
912
913 npe = nprow*npcol
914 nb = matrix%matrix_struct%nrow_block
915 nn = max(n, nb, 2)
916 np0 = numroc(nn, nb, 0, 0, nprow)
917 nq0 = max(numroc(nn, nb, 0, 0, npcol), nb)
918
919 IF (needs_evecs) THEN
920 lwork = 5*n + max(5*nn, np0*nq0) + iceil(neig_local, npe)*nn + 2*nb*nb + &
921 int(work_syevx_local*real((neig_local - 1)*n, dp)) !!!! allocates a full matrix on every CPU !!!!!
922 ELSE
923 lwork = 5*n + max(5*nn, nb*(np0 + 1))
924 END IF
925 liwork = 6*max(n, npe + 1, 4)
926
927 ALLOCATE (gap(npe))
928 gap = 0.0_dp
929 ALLOCATE (iclustr(2*npe))
930 iclustr = 0
931 ALLOCATE (ifail(n))
932 ifail = 0
933 ALLOCATE (iwork(liwork))
934 ALLOCATE (work(lwork))
935
936 ! ScaLAPACK takes advantage of IEEE754 exceptions for speedup.
937 ! Therefore, we disable floating point traps temporarily.
938#if defined (__HAS_IEEE_EXCEPTIONS)
939 CALL ieee_get_halting_mode(ieee_all, halt)
940 CALL ieee_set_halting_mode(ieee_all, .false.)
941#endif
942 CALL pdsyevx(job_type, "I", "U", n, a(1, 1), 1, 1, desca, vl, vu, 1, neig_local, abstol, m, nz, w(1), orfac, &
943 z(1, 1), 1, 1, descz, work(1), lwork, iwork(1), liwork, ifail(1), iclustr(1), gap, info)
944#if defined (__HAS_IEEE_EXCEPTIONS)
945 CALL ieee_set_halting_mode(ieee_all, halt)
946#endif
947
948 ! Error handling
949 IF (info /= 0) THEN
950 IF (ionode) THEN
951 output_unit = cp_logger_get_unit_nr(logger, local=.false.)
952 WRITE (unit=output_unit, fmt="(/,(T3,A,T12,1X,I10))") &
953 "info = ", info, &
954 "lwork = ", lwork, &
955 "liwork = ", liwork, &
956 "nz = ", nz
957 IF (info > 0) THEN
958 WRITE (unit=output_unit, fmt="(/,T3,A,(T12,6(1X,I10)))") &
959 "ifail = ", ifail
960 WRITE (unit=output_unit, fmt="(/,T3,A,(T12,6(1X,I10)))") &
961 "iclustr = ", iclustr
962 WRITE (unit=output_unit, fmt="(/,T3,A,(T12,6(1X,E10.3)))") &
963 "gap = ", gap
964 END IF
965 END IF
966 cpabort("ERROR in PDSYEVX (ScaLAPACK)")
967 END IF
968
969 ! Release work storage
970 DEALLOCATE (gap)
971 DEALLOCATE (iclustr)
972
973#else
974
975 a => matrix%local_data
976 IF (needs_evecs) THEN
977 z => eigenvectors%local_data
978 ELSE
979 ! z will not be referenced
980 z => matrix%local_data
981 END IF
982
983 ! Get the optimal work storage size
984
985 nb = max(ilaenv(1, "DSYTRD", "U", n, -1, -1, -1), &
986 ilaenv(1, "DORMTR", "U", n, -1, -1, -1))
987
988 lwork = max((nb + 3)*n, 8*n) + n ! sun bug fix
989 liwork = 5*n
990
991 ALLOCATE (ifail(n))
992 ifail = 0
993 ALLOCATE (iwork(liwork))
994 ALLOCATE (work(lwork))
995 info = 0
996 nla = SIZE(a, 1)
997 nlz = SIZE(z, 1)
998
999 ! LAPACK takes advantage of IEEE754 exceptions for speedup.
1000 ! Therefore, we disable floating point traps temporarily.
1001#if defined (__HAS_IEEE_EXCEPTIONS)
1002 CALL ieee_get_halting_mode(ieee_all, halt)
1003 CALL ieee_set_halting_mode(ieee_all, .false.)
1004#endif
1005 CALL dsyevx(job_type, "I", "U", n, a(1, 1), nla, vl, vu, 1, neig_local, &
1006 abstol, m, w, z(1, 1), nlz, work(1), lwork, iwork(1), ifail(1), info)
1007#if defined (__HAS_IEEE_EXCEPTIONS)
1008 CALL ieee_set_halting_mode(ieee_all, halt)
1009#endif
1010
1011 ! Error handling
1012 IF (info /= 0) THEN
1013 output_unit = cp_logger_get_unit_nr(logger, local=.false.)
1014 WRITE (unit=output_unit, fmt="(/,(T3,A,T12,1X,I10))") &
1015 "info = ", info
1016 IF (info > 0) THEN
1017 WRITE (unit=output_unit, fmt="(/,T3,A,(T12,6(1X,I10)))") &
1018 "ifail = ", ifail
1019 END IF
1020 cpabort("Error in DSYEVX (ScaLAPACK)")
1021 END IF
1022
1023#endif
1024 ! Release work storage
1025 DEALLOCATE (ifail)
1026 DEALLOCATE (iwork)
1027 DEALLOCATE (work)
1028 eigenvalues(1:neig_local) = w(1:neig_local)
1029 DEALLOCATE (w)
1030
1031 IF (needs_evecs) CALL check_diag(matrix, eigenvectors, neig_local)
1032
1033 CALL timestop(handle)
1034
1035 END SUBROUTINE cp_fm_syevx
1036
1037! **************************************************************************************************
1038!> \brief decomposes a quadratic matrix into its singular value decomposition
1039!> \param matrix_a ...
1040!> \param matrix_eigvl ...
1041!> \param matrix_eigvr_t ...
1042!> \param eigval ...
1043!> \param info ...
1044!> \author Maximilian Graml
1045! **************************************************************************************************
1046 SUBROUTINE cp_fm_svd(matrix_a, matrix_eigvl, matrix_eigvr_t, eigval, info)
1047
1048 TYPE(cp_fm_type), INTENT(IN) :: matrix_a
1049 TYPE(cp_fm_type), INTENT(INOUT) :: matrix_eigvl, matrix_eigvr_t
1050 REAL(kind=dp), DIMENSION(:), POINTER, &
1051 INTENT(INOUT) :: eigval
1052 INTEGER, INTENT(OUT), OPTIONAL :: info
1053
1054 CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_fm_svd'
1055
1056 INTEGER :: handle, n, m, myinfo, lwork
1057 REAL(kind=dp), DIMENSION(:, :), POINTER :: a
1058 TYPE(cp_fm_type) :: matrix_lu
1059 REAL(kind=dp), DIMENSION(:), POINTER :: work
1060 REAL(kind=dp), TARGET :: w(1)
1061#if defined(__parallel)
1062 INTEGER, DIMENSION(9) :: desca, descu, descvt
1063#endif
1064
1065 CALL timeset(routinen, handle)
1066
1067 CALL cp_fm_create(matrix=matrix_lu, &
1068 matrix_struct=matrix_a%matrix_struct, &
1069 name="A_lu"//trim(adjustl(cp_to_string(1)))//"MATRIX")
1070 CALL cp_fm_to_fm(matrix_a, matrix_lu)
1071 a => matrix_lu%local_data
1072 m = matrix_lu%matrix_struct%nrow_global
1073 n = matrix_lu%matrix_struct%ncol_global
1074 ! Assert that incoming matrix is quadratic
1075 cpassert(m == n)
1076
1077 ! Prepare for workspace queries
1078 myinfo = 0
1079 lwork = -1
1080 work => w
1081#if defined(__parallel)
1082 ! To do: That might need a redistribution check as in cp_fm_syevd
1083 desca(:) = matrix_lu%matrix_struct%descriptor(:)
1084 descu(:) = matrix_eigvl%matrix_struct%descriptor(:)
1085 descvt(:) = matrix_eigvr_t%matrix_struct%descriptor(:)
1086
1087 ! Workspace query
1088 CALL pdgesvd('V', 'V', m, m, matrix_lu%local_data, 1, 1, desca, eigval, matrix_eigvl%local_data, &
1089 1, 1, descu, matrix_eigvr_t%local_data, 1, 1, descvt, work, lwork, myinfo)
1090
1091 IF (matrix_lu%matrix_struct%para_env%is_source() .AND. (myinfo /= 0)) THEN
1092 CALL cp_fm_error("ERROR in PDGESVD: Work space query failed", myinfo, PRESENT(info))
1093 END IF
1094
1095 lwork = nint(work(1))
1096 ALLOCATE (work(lwork))
1097
1098 CALL pdgesvd('V', 'V', m, m, matrix_lu%local_data, 1, 1, desca, eigval, matrix_eigvl%local_data, &
1099 1, 1, descu, matrix_eigvr_t%local_data, 1, 1, descvt, work, lwork, myinfo)
1100
1101 IF (matrix_lu%matrix_struct%para_env%is_source() .AND. (myinfo /= 0)) THEN
1102 CALL cp_fm_error("ERROR in PDGESVD: Matrix diagonalization failed", myinfo, PRESENT(info))
1103 END IF
1104#else
1105 ! Workspace query
1106 CALL dgesvd('S', 'S', m, m, matrix_lu%local_data, m, eigval, matrix_eigvl%local_data, &
1107 m, matrix_eigvr_t%local_data, m, work, lwork, myinfo)
1108
1109 IF (myinfo /= 0) THEN
1110 CALL cp_fm_error("ERROR in DGESVD: Work space query failed", myinfo, PRESENT(info))
1111 END IF
1112
1113 lwork = nint(work(1))
1114 ALLOCATE (work(lwork))
1115
1116 CALL dgesvd('S', 'S', m, m, matrix_lu%local_data, m, eigval, matrix_eigvl%local_data, &
1117 m, matrix_eigvr_t%local_data, m, work, lwork, myinfo)
1118
1119 IF (myinfo /= 0) THEN
1120 CALL cp_fm_error("ERROR in DGESVD: Matrix diagonalization failed", myinfo, PRESENT(info))
1121 END IF
1122
1123#endif
1124 ! Release intermediary matrices
1125 DEALLOCATE (work)
1126 CALL cp_fm_release(matrix_lu)
1127
1128 IF (PRESENT(info)) info = myinfo
1129
1130 CALL timestop(handle)
1131 END SUBROUTINE cp_fm_svd
1132
1133! **************************************************************************************************
1134!> \brief ...
1135!> \param matrix ...
1136!> \param work ...
1137!> \param exponent ...
1138!> \param threshold ...
1139!> \param n_dependent ...
1140!> \param verbose ...
1141!> \param eigvals ...
1142! **************************************************************************************************
1143 SUBROUTINE cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
1144
1145 ! Raise the real symmetric n by n matrix to the power given by
1146 ! the exponent. All eigenvectors with a corresponding eigenvalue lower
1147 ! than threshold are quenched. result in matrix
1148
1149 ! - Creation (29.03.1999, Matthias Krack)
1150 ! - Parallelised using BLACS and ScaLAPACK (06.06.2001,MK)
1151
1152 TYPE(cp_fm_type), INTENT(IN) :: matrix, work
1153 REAL(kind=dp), INTENT(IN) :: exponent, threshold
1154 INTEGER, INTENT(OUT) :: n_dependent
1155 LOGICAL, INTENT(IN), OPTIONAL :: verbose
1156 REAL(kind=dp), DIMENSION(2), INTENT(OUT), &
1157 OPTIONAL :: eigvals
1158
1159 CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_fm_power'
1160
1161 INTEGER :: handle, icol_global, &
1162 mypcol, myprow, &
1163 ncol_global, nrow_global
1164 LOGICAL :: my_verbose
1165 REAL(kind=dp) :: condition_number, f, p
1166 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: eigenvalues
1167 REAL(kind=dp), DIMENSION(:, :), POINTER :: eigenvectors
1168 TYPE(cp_blacs_env_type), POINTER :: context
1169
1170#if defined(__parallel)
1171 INTEGER :: icol_local, ipcol, iprow, irow_global, irow_local
1172#endif
1173
1174 CALL timeset(routinen, handle)
1175
1176 my_verbose = .false.
1177 IF (PRESENT(verbose)) my_verbose = verbose
1178
1179 context => matrix%matrix_struct%context
1180 myprow = context%mepos(1)
1181 mypcol = context%mepos(2)
1182 n_dependent = 0
1183 p = 0.5_dp*exponent
1184
1185 nrow_global = matrix%matrix_struct%nrow_global
1186 ncol_global = matrix%matrix_struct%ncol_global
1187
1188 ALLOCATE (eigenvalues(ncol_global))
1189 eigenvalues(:) = 0.0_dp
1190
1191 ! Compute the eigenvectors and eigenvalues
1192
1193 CALL choose_eigv_solver(matrix, work, eigenvalues)
1194
1195 IF (PRESENT(eigvals)) THEN
1196 eigvals(1) = eigenvalues(1)
1197 eigvals(2) = eigenvalues(ncol_global)
1198 END IF
1199
1200#if defined(__parallel)
1201 eigenvectors => work%local_data
1202
1203 ! Build matrix**exponent with eigenvector quenching
1204
1205 DO icol_global = 1, ncol_global
1206
1207 IF (eigenvalues(icol_global) < threshold) THEN
1208
1209 n_dependent = n_dependent + 1
1210
1211 ipcol = work%matrix_struct%g2p_col(icol_global)
1212
1213 IF (mypcol == ipcol) THEN
1214 icol_local = work%matrix_struct%g2l_col(icol_global)
1215 DO irow_global = 1, nrow_global
1216 iprow = work%matrix_struct%g2p_row(irow_global)
1217 IF (myprow == iprow) THEN
1218 irow_local = work%matrix_struct%g2l_row(irow_global)
1219 eigenvectors(irow_local, icol_local) = 0.0_dp
1220 END IF
1221 END DO
1222 END IF
1223
1224 ELSE
1225
1226 f = eigenvalues(icol_global)**p
1227
1228 ipcol = work%matrix_struct%g2p_col(icol_global)
1229
1230 IF (mypcol == ipcol) THEN
1231 icol_local = work%matrix_struct%g2l_col(icol_global)
1232 DO irow_global = 1, nrow_global
1233 iprow = work%matrix_struct%g2p_row(irow_global)
1234 IF (myprow == iprow) THEN
1235 irow_local = work%matrix_struct%g2l_row(irow_global)
1236 eigenvectors(irow_local, icol_local) = &
1237 f*eigenvectors(irow_local, icol_local)
1238 END IF
1239 END DO
1240 END IF
1241
1242 END IF
1243
1244 END DO
1245
1246#else
1247
1248 eigenvectors => work%local_data
1249
1250 ! Build matrix**exponent with eigenvector quenching
1251
1252 DO icol_global = 1, ncol_global
1253
1254 IF (eigenvalues(icol_global) < threshold) THEN
1255
1256 n_dependent = n_dependent + 1
1257 eigenvectors(1:nrow_global, icol_global) = 0.0_dp
1258
1259 ELSE
1260
1261 f = eigenvalues(icol_global)**p
1262 eigenvectors(1:nrow_global, icol_global) = &
1263 f*eigenvectors(1:nrow_global, icol_global)
1264
1265 END IF
1266
1267 END DO
1268
1269#endif
1270 CALL cp_fm_syrk("U", "N", ncol_global, 1.0_dp, work, 1, 1, 0.0_dp, matrix)
1271 CALL cp_fm_uplo_to_full(matrix, work)
1272
1273 ! Print some warnings/notes
1274 IF (matrix%matrix_struct%para_env%is_source() .AND. my_verbose) THEN
1275 condition_number = abs(eigenvalues(ncol_global)/eigenvalues(1))
1276 WRITE (unit=cp_logger_get_default_unit_nr(), fmt="(/,(T2,A,ES15.6))") &
1277 "CP_FM_POWER: smallest eigenvalue:", eigenvalues(1), &
1278 "CP_FM_POWER: largest eigenvalue: ", eigenvalues(ncol_global), &
1279 "CP_FM_POWER: condition number: ", condition_number
1280 IF (eigenvalues(1) <= 0.0_dp) THEN
1281 WRITE (unit=cp_logger_get_default_unit_nr(), fmt="(/,T2,A)") &
1282 "WARNING: matrix has a negative eigenvalue, tighten EPS_DEFAULT"
1283 END IF
1284 IF (condition_number > 1.0e12_dp) THEN
1285 WRITE (unit=cp_logger_get_default_unit_nr(), fmt="(/,T2,A)") &
1286 "WARNING: high condition number => possibly ill-conditioned matrix"
1287 END IF
1288 END IF
1289
1290 DEALLOCATE (eigenvalues)
1291
1292 CALL timestop(handle)
1293
1294 END SUBROUTINE cp_fm_power
1295
1296! **************************************************************************************************
1297!> \brief ...
1298!> \param matrix ...
1299!> \param eigenvectors ...
1300!> \param eigval ...
1301!> \param thresh ...
1302!> \param start_sec_block ...
1303! **************************************************************************************************
1304 SUBROUTINE cp_fm_block_jacobi(matrix, eigenvectors, eigval, thresh, start_sec_block)
1305
1306 ! Calculates block diagonalization of a full symmetric matrix
1307 ! It has its origin in cp_fm_syevx. This routine rotates only elements
1308 ! which are larger than a threshold values "thresh".
1309 ! start_sec_block is the start of the second block.
1310 ! IT DOES ONLY ONE SWEEP!
1311
1312 ! - Creation (07.10.2002, Martin Fengler)
1313 ! - Cosmetics (05.04.06, MK)
1314
1315 TYPE(cp_fm_type), INTENT(IN) :: eigenvectors, matrix
1316 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigval
1317 INTEGER, INTENT(IN) :: start_sec_block
1318 REAL(kind=dp), INTENT(IN) :: thresh
1319
1320 CHARACTER(len=*), PARAMETER :: routinen = 'cp_fm_block_jacobi'
1321
1322 INTEGER :: handle
1323 REAL(kind=dp), DIMENSION(:, :), POINTER :: a, ev
1324
1325 REAL(kind=dp) :: tan_theta, tau, c, s
1326 INTEGER :: q, p, n
1327 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: c_ip
1328
1329#if defined(__parallel)
1330 TYPE(cp_blacs_env_type), POINTER :: context
1331
1332 INTEGER :: nprow, npcol, block_dim_row, block_dim_col, info, &
1333 ev_row_block_size, iam, mynumrows, mype, npe, q_loc
1334 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: a_loc, ev_loc
1335 INTEGER, DIMENSION(9) :: desca, descz, &
1336 desc_a_block, &
1337 desc_ev_loc
1338 TYPE(mp_comm_type):: allgrp
1339 TYPE(cp_blacs_type) :: ictxt_loc
1340 INTEGER, EXTERNAL :: numroc
1341#endif
1342
1343 ! -------------------------------------------------------------------------
1344
1345 CALL timeset(routinen, handle)
1346
1347#if defined(__parallel)
1348 context => matrix%matrix_struct%context
1349 allgrp = matrix%matrix_struct%para_env
1350
1351 nprow = context%num_pe(1)
1352 npcol = context%num_pe(2)
1353
1354 n = matrix%matrix_struct%nrow_global
1355
1356 a => matrix%local_data
1357 desca(:) = matrix%matrix_struct%descriptor(:)
1358 ev => eigenvectors%local_data
1359 descz(:) = eigenvectors%matrix_struct%descriptor(:)
1360
1361 ! Copy the block to be rotated to the master process firstly and broadcast to all processes
1362 ! start_sec_block defines where the second block starts!
1363 ! Block will be processed together with the OO block
1364
1365 block_dim_row = start_sec_block - 1
1366 block_dim_col = n - block_dim_row
1367 ALLOCATE (a_loc(block_dim_row, block_dim_col))
1368
1369 mype = matrix%matrix_struct%para_env%mepos
1370 npe = matrix%matrix_struct%para_env%num_pe
1371 ! Get a new context
1372 CALL ictxt_loc%gridinit(matrix%matrix_struct%para_env, 'R', nprow*npcol, 1)
1373
1374 CALL descinit(desc_a_block, block_dim_row, block_dim_col, block_dim_row, &
1375 block_dim_col, 0, 0, ictxt_loc%get_handle(), block_dim_row, info)
1376
1377 CALL pdgemr2d(block_dim_row, block_dim_col, a, 1, start_sec_block, desca, &
1378 a_loc, 1, 1, desc_a_block, context%get_handle())
1379 ! Only the master (root) process received data yet
1380 CALL allgrp%bcast(a_loc, 0)
1381
1382 ! Since each process owns now the upper block, the eigenvectors can be re-sorted in such a way that
1383 ! each process has a NN*1 grid, i.e. the process owns a bunch of rows which can be modified locally
1384
1385 ! Initialize distribution of the eigenvectors
1386 iam = mype
1387 ev_row_block_size = n/(nprow*npcol)
1388 mynumrows = numroc(n, ev_row_block_size, iam, 0, nprow*npcol)
1389
1390 ALLOCATE (ev_loc(mynumrows, n), c_ip(mynumrows))
1391
1392 CALL descinit(desc_ev_loc, n, n, ev_row_block_size, n, 0, 0, ictxt_loc%get_handle(), &
1393 mynumrows, info)
1394
1395 CALL pdgemr2d(n, n, ev, 1, 1, descz, ev_loc, 1, 1, desc_ev_loc, context%get_handle())
1396
1397 ! Start block diagonalization of matrix
1398
1399 q_loc = 0
1400
1401 DO q = start_sec_block, n
1402 q_loc = q_loc + 1
1403 DO p = 1, (start_sec_block - 1)
1404
1405 IF (abs(a_loc(p, q_loc)) > thresh) THEN
1406
1407 tau = (eigval(q) - eigval(p))/(2.0_dp*a_loc(p, q_loc))
1408
1409 tan_theta = sign(1.0_dp, tau)/(abs(tau) + sqrt(1.0_dp + tau*tau))
1410
1411 ! Cos(theta)
1412 c = 1.0_dp/sqrt(1.0_dp + tan_theta*tan_theta)
1413 s = tan_theta*c
1414
1415 ! Calculate eigenvectors: Q*J
1416 ! c_ip = c*EV_loc(:,p) - s*EV_loc(:,q)
1417 ! c_iq = s*EV_loc(:,p) + c*EV_loc(:,q)
1418 ! EV(:,p) = c_ip
1419 ! EV(:,q) = c_iq
1420 CALL dcopy(mynumrows, ev_loc(1, p), 1, c_ip(1), 1)
1421 CALL dscal(mynumrows, c, ev_loc(1, p), 1)
1422 CALL daxpy(mynumrows, -s, ev_loc(1, q), 1, ev_loc(1, p), 1)
1423 CALL dscal(mynumrows, c, ev_loc(1, q), 1)
1424 CALL daxpy(mynumrows, s, c_ip(1), 1, ev_loc(1, q), 1)
1425
1426 END IF
1427
1428 END DO
1429 END DO
1430
1431 ! Copy eigenvectors back to the original distribution
1432 CALL pdgemr2d(n, n, ev_loc, 1, 1, desc_ev_loc, ev, 1, 1, descz, context%get_handle())
1433
1434 ! Release work storage
1435 DEALLOCATE (a_loc, ev_loc, c_ip)
1436
1437 CALL ictxt_loc%gridexit()
1438
1439#else
1440
1441 n = matrix%matrix_struct%nrow_global
1442
1443 ALLOCATE (c_ip(n)) ! Local eigenvalue vector
1444
1445 a => matrix%local_data ! Contains the Matrix to be worked on
1446 ev => eigenvectors%local_data ! Contains the eigenvectors up to blocksize, rest is garbage
1447
1448 ! Start matrix diagonalization
1449
1450 tan_theta = 0.0_dp
1451 tau = 0.0_dp
1452
1453 DO q = start_sec_block, n
1454 DO p = 1, (start_sec_block - 1)
1455
1456 IF (abs(a(p, q)) > thresh) THEN
1457
1458 tau = (eigval(q) - eigval(p))/(2.0_dp*a(p, q))
1459
1460 tan_theta = sign(1.0_dp, tau)/(abs(tau) + sqrt(1.0_dp + tau*tau))
1461
1462 ! Cos(theta)
1463 c = 1.0_dp/sqrt(1.0_dp + tan_theta*tan_theta)
1464 s = tan_theta*c
1465
1466 ! Calculate eigenvectors: Q*J
1467 ! c_ip = c*EV(:,p) - s*EV(:,q)
1468 ! c_iq = s*EV(:,p) + c*EV(:,q)
1469 ! EV(:,p) = c_ip
1470 ! EV(:,q) = c_iq
1471 CALL dcopy(n, ev(1, p), 1, c_ip(1), 1)
1472 CALL dscal(n, c, ev(1, p), 1)
1473 CALL daxpy(n, -s, ev(1, q), 1, ev(1, p), 1)
1474 CALL dscal(n, c, ev(1, q), 1)
1475 CALL daxpy(n, s, c_ip(1), 1, ev(1, q), 1)
1476
1477 END IF
1478
1479 END DO
1480 END DO
1481
1482 ! Release work storage
1483
1484 DEALLOCATE (c_ip)
1485
1486#endif
1487
1488 CALL timestop(handle)
1489
1490 END SUBROUTINE cp_fm_block_jacobi
1491
1492! **************************************************************************************************
1493!> \brief General Eigenvalue Problem AX = BXE.
1494!> Use cuSOLVERMp directly when requested and large enough; otherwise
1495!> reduce the problem through a Cholesky decomposition of B.
1496!> \param amatrix ...
1497!> \param bmatrix ...
1498!> \param eigenvectors ...
1499!> \param eigenvalues ...
1500!> \param work ...
1501! **************************************************************************************************
1502 SUBROUTINE cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
1503
1504 TYPE(cp_fm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
1505 REAL(kind=dp), DIMENSION(:) :: eigenvalues
1506 TYPE(cp_fm_type), INTENT(IN) :: work
1507
1508 CHARACTER(len=*), PARAMETER :: routinen = 'cp_fm_geeig'
1509
1510 INTEGER :: handle, nao, nmo
1511 LOGICAL :: check_eigenvectors
1512 TYPE(cp_fm_type) :: overlap_check, scratch_check
1513
1514 CALL timeset(routinen, handle)
1515
1516 CALL cp_fm_get_info(amatrix, nrow_global=nao)
1517 nmo = SIZE(eigenvalues)
1518 check_eigenvectors = diag_check_requested()
1519
1520 IF (diag_type == fm_diag_type_cusolver .AND. direct_generalized_diagonalization .AND. &
1521 nao >= cusolver_n_min) THEN
1522 ! Use cuSolverMP generalized eigenvalue solver without a CP2K-side
1523 ! Cholesky reduction.
1524 ! Use work as intermediate buffer since eigenvectors may be smaller (nao x nmo)
1525 IF (check_eigenvectors) THEN
1526 CALL cp_fm_create(overlap_check, bmatrix%matrix_struct)
1527 CALL cp_fm_create(scratch_check, bmatrix%matrix_struct)
1528 CALL cp_fm_to_fm(bmatrix, overlap_check)
1529 END IF
1530 CALL cp_fm_general_cusolver(amatrix, bmatrix, work, eigenvalues)
1531 IF (check_eigenvectors) THEN
1532 CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
1533 CALL cp_fm_release(scratch_check)
1534 CALL cp_fm_release(overlap_check)
1535 END IF
1536 CALL cp_fm_to_fm(work, eigenvectors, nmo)
1537#if defined(__parallel)
1538 ELSE IF (diag_type == fm_diag_type_scalapack .AND. direct_generalized_diagonalization) THEN
1539 ! Use ScaLAPACK generalized eigenvalue solver without a CP2K-side
1540 ! Cholesky reduction.
1541 IF (check_eigenvectors) THEN
1542 CALL cp_fm_create(overlap_check, bmatrix%matrix_struct)
1543 CALL cp_fm_create(scratch_check, bmatrix%matrix_struct)
1544 CALL cp_fm_to_fm(bmatrix, overlap_check)
1545 END IF
1546 CALL cp_fm_geeig_scalapack(amatrix, bmatrix, work, eigenvalues)
1547 IF (check_eigenvectors) THEN
1548 CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
1549 CALL cp_fm_release(scratch_check)
1550 CALL cp_fm_release(overlap_check)
1551 END IF
1552 CALL cp_fm_to_fm(work, eigenvectors, nmo)
1553#endif
1554#if defined(__DLAF)
1555 ELSE IF (diag_type == fm_diag_type_dlaf .AND. direct_generalized_diagonalization .AND. &
1556 nao >= dlaf_neigvec_min) THEN
1557 ! Use DLA-Future generalized eigenvalue solver for large matrices
1558 IF (check_eigenvectors) THEN
1559 CALL cp_fm_create(overlap_check, bmatrix%matrix_struct)
1560 CALL cp_fm_create(scratch_check, bmatrix%matrix_struct)
1561 CALL cp_fm_to_fm(bmatrix, overlap_check)
1562 END IF
1563 CALL cp_fm_diag_gen_dlaf(amatrix, bmatrix, work, eigenvalues)
1564 IF (check_eigenvectors) THEN
1565 CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
1566 CALL cp_fm_release(scratch_check)
1567 CALL cp_fm_release(overlap_check)
1568 END IF
1569 CALL cp_fm_to_fm(work, eigenvectors, nmo)
1570#endif
1571 ELSE
1572 ! Cholesky decompose S=U(T)U
1573 CALL cp_fm_cholesky_decompose(bmatrix)
1574 ! Invert to get U^(-1)
1575 CALL cp_fm_triangular_invert(bmatrix)
1576 ! Reduce to get U^(-T) * H * U^(-1)
1577 CALL cp_fm_triangular_multiply(bmatrix, amatrix, side="R")
1578 CALL cp_fm_triangular_multiply(bmatrix, amatrix, transpose_tr=.true.)
1579 ! Diagonalize
1580 CALL choose_eigv_solver(matrix=amatrix, eigenvectors=work, &
1581 eigenvalues=eigenvalues)
1582 ! Restore vectors C = U^(-1) * C*
1583 CALL cp_fm_triangular_multiply(bmatrix, work)
1584 CALL cp_fm_to_fm(work, eigenvectors, nmo)
1585 END IF
1586
1587 CALL timestop(handle)
1588
1589 END SUBROUTINE cp_fm_geeig
1590
1591! **************************************************************************************************
1592!> \brief General Eigenvalue Problem AX = BXE using ScaLAPACK PDSYGVX.
1593!> \param amatrix ...
1594!> \param bmatrix ...
1595!> \param eigenvectors ...
1596!> \param eigenvalues ...
1597! **************************************************************************************************
1598 SUBROUTINE cp_fm_geeig_scalapack(amatrix, bmatrix, eigenvectors, eigenvalues)
1599
1600 TYPE(cp_fm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
1601 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
1602
1603 CHARACTER(len=*), PARAMETER :: routinen = 'cp_fm_geeig_scalapack'
1604
1605#if defined(__parallel)
1606 REAL(kind=dp), PARAMETER :: orfac = -1.0_dp, &
1607 vl = 0.0_dp, &
1608 vu = 0.0_dp
1609
1610 INTEGER :: handle, info, liwork, lwork, m, n, nb, &
1611 neig, npcol, nprow, nz
1612 INTEGER, DIMENSION(9) :: desca, descb, descz
1613 INTEGER, DIMENSION(:), ALLOCATABLE :: iclustr, ifail, iwork
1614 REAL(kind=dp) :: abstol
1615 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: gap, w, work
1616 REAL(kind=dp), DIMENSION(:, :), POINTER :: a, b, z
1617
1618 INTEGER :: mq0, nn, np0, npe
1619 INTEGER, EXTERNAL :: iceil, numroc
1620 REAL(kind=dp), EXTERNAL :: dlamch
1621#if defined (__HAS_IEEE_EXCEPTIONS)
1622 LOGICAL, DIMENSION(5) :: halt
1623#endif
1624#else
1625 INTEGER :: handle
1626#endif
1627
1628 CALL timeset(routinen, handle)
1629
1630#if defined(__parallel)
1631 n = amatrix%matrix_struct%nrow_global
1632 neig = min(SIZE(eigenvalues), n)
1633
1634 IF (neig == 0) THEN
1635 CALL timestop(handle)
1636 RETURN
1637 END IF
1638
1639 IF (amatrix%matrix_struct%nrow_block /= amatrix%matrix_struct%ncol_block) THEN
1640 cpabort("ERROR in "//routinen//": Invalid blocksize (no square blocks) found")
1641 END IF
1642
1643 a => amatrix%local_data
1644 b => bmatrix%local_data
1645 z => eigenvectors%local_data
1646 desca(:) = amatrix%matrix_struct%descriptor(:)
1647 descb(:) = bmatrix%matrix_struct%descriptor(:)
1648 descz(:) = eigenvectors%matrix_struct%descriptor(:)
1649
1650 nprow = amatrix%matrix_struct%context%num_pe(1)
1651 npcol = amatrix%matrix_struct%context%num_pe(2)
1652 npe = nprow*npcol
1653 nb = amatrix%matrix_struct%nrow_block
1654 nn = max(n, nb, 2)
1655 np0 = numroc(nn, nb, 0, 0, nprow)
1656 mq0 = max(numroc(nn, nb, 0, 0, npcol), nb)
1657
1658 lwork = 5*n + max(5*nn, np0*mq0 + 2*nb*nb) + iceil(neig, npe)*nn + &
1659 max(0, neig - 1)*n
1660 liwork = 6*max(n, npe + 1, 4)
1661
1662 ALLOCATE (gap(npe))
1663 gap = 0.0_dp
1664 ALLOCATE (iclustr(2*npe))
1665 iclustr = 0
1666 ALLOCATE (ifail(n))
1667 ifail = 0
1668 ALLOCATE (iwork(liwork))
1669 ALLOCATE (w(n))
1670 ALLOCATE (work(lwork))
1671
1672 abstol = 2.0_dp*dlamch("S")
1673
1674#if defined (__HAS_IEEE_EXCEPTIONS)
1675 CALL ieee_get_halting_mode(ieee_all, halt)
1676 CALL ieee_set_halting_mode(ieee_all, .false.)
1677#endif
1678 CALL pdsygvx(1, "V", "I", "U", n, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb, &
1679 vl, vu, 1, neig, abstol, m, nz, w(1), orfac, z(1, 1), 1, 1, descz, &
1680 work(1), lwork, iwork(1), liwork, ifail(1), iclustr(1), gap(1), info)
1681#if defined (__HAS_IEEE_EXCEPTIONS)
1682 CALL ieee_set_halting_mode(ieee_all, halt)
1683#endif
1684
1685 IF (info /= 0 .OR. m < neig .OR. nz < neig) THEN
1686 cpabort("ERROR in PDSYGVX (ScaLAPACK), info="//trim(cp_to_string(info)))
1687 END IF
1688
1689 eigenvalues(:) = 0.0_dp
1690 eigenvalues(1:neig) = w(1:neig)
1691
1692 DEALLOCATE (gap, iclustr, ifail, iwork, w, work)
1693#else
1694 mark_used(amatrix)
1695 mark_used(bmatrix)
1696 mark_used(eigenvectors)
1697 mark_used(eigenvalues)
1698 cpabort("ERROR in "//routinen//": PDSYGVX requested without ScaLAPACK support")
1699#endif
1700
1701 CALL timestop(handle)
1702
1703 END SUBROUTINE cp_fm_geeig_scalapack
1704
1705! **************************************************************************************************
1706!> \brief General Eigenvalue Problem AX = BXE
1707!> Use canonical diagonalization : U*s**(-1/2)
1708!> \param amatrix ...
1709!> \param bmatrix ...
1710!> \param eigenvectors ...
1711!> \param eigenvalues ...
1712!> \param work ...
1713!> \param epseig ...
1714!> \param nmo_retained ...
1715! **************************************************************************************************
1716 SUBROUTINE cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig, &
1717 nmo_retained)
1718
1719 TYPE(cp_fm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
1720 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
1721 TYPE(cp_fm_type), INTENT(IN) :: work
1722 REAL(kind=dp), INTENT(IN) :: epseig
1723 INTEGER, INTENT(OUT), OPTIONAL :: nmo_retained
1724
1725 CHARACTER(len=*), PARAMETER :: routinen = 'cp_fm_geeig_canon'
1726
1727 INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
1728 nmo, nx
1729 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
1730
1731 CALL timeset(routinen, handle)
1732
1733 ! Test sizees
1734 CALL cp_fm_get_info(amatrix, nrow_global=nao)
1735 nmo = SIZE(eigenvalues)
1736 ALLOCATE (evals(nao))
1737
1738 ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
1739 CALL cp_fm_scale(-1.0_dp, bmatrix)
1740 CALL choose_eigv_solver(matrix=bmatrix, eigenvectors=work, eigenvalues=evals)
1741 evals(:) = -evals(:)
1742 nc = nao
1743 DO i = 1, nao
1744 IF (evals(i) < epseig) THEN
1745 nc = i - 1
1746 EXIT
1747 END IF
1748 END DO
1749 cpassert(nc /= 0)
1750
1751 IF (nc /= nao) THEN
1752 IF (nc < nmo) THEN
1753 ! Copy NULL space definition to last vectors of eigenvectors (if needed)
1754 ncol = nmo - nc
1755 CALL cp_fm_to_fm(work, eigenvectors, ncol, nc + 1, nc + 1)
1756 END IF
1757 ! Set NULL space in eigenvector matrix of S to zero
1758 DO icol = nc + 1, nao
1759 DO irow = 1, nao
1760 CALL cp_fm_set_element(work, irow, icol, 0.0_dp)
1761 END DO
1762 END DO
1763 ! Set small eigenvalues to a dummy save value
1764 evals(nc + 1:nao) = 1.0_dp
1765 END IF
1766 ! Calculate U*s**(-1/2)
1767 evals(:) = 1.0_dp/sqrt(evals(:))
1768 CALL cp_fm_column_scale(work, evals)
1769 ! Reduce to get U^(-T) * H * U^(-1)
1770 CALL cp_fm_gemm("T", "N", nao, nao, nao, 1.0_dp, work, amatrix, 0.0_dp, bmatrix)
1771 CALL cp_fm_gemm("N", "N", nao, nao, nao, 1.0_dp, bmatrix, work, 0.0_dp, amatrix)
1772 IF (nc /= nao) THEN
1773 ! set diagonal values to save large value
1774 DO icol = nc + 1, nao
1775 CALL cp_fm_set_element(amatrix, icol, icol, set_removed_eigval_to)
1776 END DO
1777 END IF
1778 ! Diagonalize
1779 CALL choose_eigv_solver(matrix=amatrix, eigenvectors=bmatrix, eigenvalues=eigenvalues)
1780 nx = min(nc, nmo)
1781 ! Restore vectors C = U^(-1) * C*
1782 CALL cp_fm_gemm("N", "N", nao, nx, nc, 1.0_dp, work, bmatrix, 0.0_dp, eigenvectors)
1783
1784 ! Number of basis modes that survived the linear-dependency filter. The remaining
1785 ! nao - nc entries of eigenvalues(:) are the placeholders set above.
1786 IF (PRESENT(nmo_retained)) nmo_retained = nc
1787
1788 DEALLOCATE (evals)
1789
1790 CALL timestop(handle)
1791
1792 END SUBROUTINE cp_fm_geeig_canon
1793
1794END MODULE cp_fm_diag
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
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_blacs_types
wrappers for the actual blacs calls. all functionality needed in the code should actually be provide ...
Definition cp_blacs_types.F:18
cp_cfm_elpa
Wrapper for ELPA (complex matrices, i.e. cp_cfm_type).
Definition cp_cfm_elpa.F:11
cp_cfm_elpa::set_elpa_c_kernel
subroutine, public set_elpa_c_kernel(requested_kernel)
Sets the active ELPA kernel for complex matrices.
Definition cp_cfm_elpa.F:191
cp_dlaf_utils_api
Definition cp_dlaf_utils_api.F:8
cp_dlaf_utils_api::cp_dlaf_finalize
subroutine, public cp_dlaf_finalize()
Finalize DLA-Future and pika runtime.
Definition cp_dlaf_utils_api.F:56
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_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_gemm
subroutine, public cp_fm_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)
computes matrix_c = beta * matrix_c + alpha * ( matrix_a ** transa ) * ( matrix_b ** transb )
Definition cp_fm_basic_linalg.F:476
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1691
cp_fm_basic_linalg::cp_fm_syrk
subroutine, public cp_fm_syrk(uplo, trans, k, alpha, matrix_a, ia, ja, beta, matrix_c)
performs a rank-k update of a symmetric matrix_c matrix_c = beta * matrix_c + alpha * matrix_a * tran...
Definition cp_fm_basic_linalg.F:666
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1585
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition cp_fm_basic_linalg.F:1503
cp_fm_basic_linalg::cp_fm_triangular_invert
subroutine, public cp_fm_triangular_invert(matrix_a, uplo_tr)
inverts a triangular matrix
Definition cp_fm_basic_linalg.F:1969
cp_fm_basic_linalg::cp_fm_triangular_multiply
subroutine, public cp_fm_triangular_multiply(triangular_matrix, matrix_b, side, transpose_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_fm_basic_linalg.F:1414
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_fm_cholesky.F:62
cp_fm_cusolver_api
Wrapper for cuSOLVERMp.
Definition cp_fm_cusolver_api.F:12
cp_fm_cusolver_api::cp_fm_general_cusolver
subroutine, public cp_fm_general_cusolver(amatrix, bmatrix, eigenvectors, eigenvalues)
Driver routine to solve generalized eigenvalue problem A*x = lambda*B*x with cuSOLVERMp.
Definition cp_fm_cusolver_api.F:113
cp_fm_cusolver_api::cp_fm_diag_cusolver
subroutine, public cp_fm_diag_cusolver(matrix, eigenvectors, eigenvalues)
Driver routine to diagonalize a FM matrix with the cuSOLVERMp library.
Definition cp_fm_cusolver_api.F:40
cp_fm_diag_utils
Auxiliary tools to redistribute cp_fm_type and cp_cfm_type matrices before and after diagonalization....
Definition cp_fm_diag_utils.F:17
cp_fm_diag_utils::cp_fm_redistribute_end
subroutine, public cp_fm_redistribute_end(matrix, eigenvectors, eig, matrix_new, eigenvectors_new)
Redistributes eigenvectors and eigenvalues back to the original communicator group.
Definition cp_fm_diag_utils.F:471
cp_fm_diag_utils::cp_fm_redistribute_start
subroutine, public cp_fm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new, caller_is_elpa, redist_info)
Determines the optimal number of CPUs for matrix diagonalization and redistributes the input matrices...
Definition cp_fm_diag_utils.F:278
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::cp_fm_block_jacobi
subroutine, public cp_fm_block_jacobi(matrix, eigenvectors, eigval, thresh, start_sec_block)
...
Definition cp_fm_diag.F:1305
cp_fm_diag::eps_check_diag_default
real(kind=dp), parameter, public eps_check_diag_default
Definition cp_fm_diag.F:85
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition cp_fm_diag.F:1144
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::cp_fm_svd
subroutine, public cp_fm_svd(matrix_a, matrix_eigvl, matrix_eigvr_t, eigval, info)
decomposes a quadratic matrix into its singular value decomposition
Definition cp_fm_diag.F:1047
cp_fm_diag::cp_fm_geeig_canon
subroutine, public cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig, nmo_retained)
General Eigenvalue Problem AX = BXE Use canonical diagonalization : U*s**(-1/2).
Definition cp_fm_diag.F:1718
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::cp_fm_syevd
subroutine, public cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
Computes all eigenvalues and vectors of a real symmetric matrix significantly faster than syevx,...
Definition cp_fm_diag.F:590
cp_fm_diag::diag_finalize
subroutine, public diag_finalize()
Finalize the diagonalization library.
Definition cp_fm_diag.F:238
cp_fm_diag::direct_generalized_diagonalization
logical, save, public direct_generalized_diagonalization
Definition cp_fm_diag.F:106
cp_fm_diag::fm_diag_type_default
integer, parameter, public fm_diag_type_default
Definition cp_fm_diag.F:121
cp_fm_diag::elpa_neigvec_min
integer, save, public elpa_neigvec_min
Definition cp_fm_diag.F:99
cp_fm_diag::diag_init
subroutine, public diag_init(diag_lib, fallback_applied, elpa_kernel, elpa_c_kernel, elpa_neigvec_min_input, elpa_qr, elpa_print, elpa_one_stage, dlaf_neigvec_min_input, eps_check_diag_input, direct_generalized_diagonalization_input, diag_lib_explicit_input)
Setup the diagonalization library to be used.
Definition cp_fm_diag.F:164
cp_fm_diag::cp_fm_geeig
subroutine, public cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE. Use cuSOLVERMp directly when requested and large enough; otherwi...
Definition cp_fm_diag.F:1503
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:262
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::cp_fm_syevx
subroutine, public cp_fm_syevx(matrix, eigenvectors, eigenvalues, neig, work_syevx)
compute eigenvalues and optionally eigenvectors of a real symmetric matrix using scalapack....
Definition cp_fm_diag.F:814
cp_fm_diag::diag_lib_explicit
logical, save, public diag_lib_explicit
Definition cp_fm_diag.F:96
cp_fm_dlaf_api
Definition cp_fm_dlaf_api.F:8
cp_fm_dlaf_api::cp_fm_diag_dlaf
subroutine, public cp_fm_diag_dlaf(matrix, eigenvectors, eigenvalues)
...
Definition cp_fm_dlaf_api.F:121
cp_fm_dlaf_api::cp_fm_diag_gen_dlaf
subroutine, public cp_fm_diag_gen_dlaf(a_matrix, b_matrix, eigenvectors, eigenvalues)
...
Definition cp_fm_dlaf_api.F:229
cp_fm_elpa
Wrapper for ELPA.
Definition cp_fm_elpa.F:12
cp_fm_elpa::set_elpa_kernel
subroutine, public set_elpa_kernel(requested_kernel)
Sets the active ELPA kernel.
Definition cp_fm_elpa.F:247
cp_fm_elpa::elpa_qr
logical, save, public elpa_qr
Definition cp_fm_elpa.F:161
cp_fm_elpa::cp_fm_diag_elpa
subroutine, public cp_fm_diag_elpa(matrix, eigenvectors, eigenvalues)
Driver routine to diagonalize a FM matrix with the ELPA library.
Definition cp_fm_elpa.F:299
cp_fm_elpa::elpa_print
logical, save, public elpa_print
Definition cp_fm_elpa.F:161
cp_fm_elpa::finalize_elpa_library
subroutine, public finalize_elpa_library()
Finalize the ELPA library.
Definition cp_fm_elpa.F:235
cp_fm_elpa::elpa_one_stage
logical, save, public elpa_one_stage
Definition cp_fm_elpa.F:168
cp_fm_elpa::initialize_elpa_library
subroutine, public initialize_elpa_library(one_stage, qr, should_print)
Initialize the ELPA library.
Definition cp_fm_elpa.F:214
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_equivalent
logical function, public cp_fm_struct_equivalent(fmstruct1, fmstruct2)
returns true if the two matrix structures are equivalent, false otherwise.
Definition cp_fm_struct.F:392
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1044
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition cp_fm_types.F:1461
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_set_element
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
Definition cp_fm_types.F:690
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:576
cp_log_handling::cp_logger_get_unit_nr
integer function, public cp_logger_get_unit_nr(logger, local)
returns the unit nr for the requested kind of log.
Definition cp_log_handling.F:506
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
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
machine::m_memory
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition machine.F:440
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
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_blacs_types::cp_blacs_type
Definition cp_blacs_types.F:30
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_comm_type
Definition message_passing.F:179