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