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