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