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