71#if defined (__parallel)
74#if defined (__HAS_IEEE_EXCEPTIONS)
75 USE ieee_exceptions,
ONLY: ieee_get_halting_mode, &
76 ieee_set_halting_mode, &
79#include "../base/base_uses.f90"
85 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'cp_fm_diag'
94 INTEGER,
SAVE :: elpa_neigvec_min = 0
98 INTEGER,
SAVE,
PUBLIC :: dlaf_neigvec_min = 0
102 REAL(kind=
dp),
SAVE :: eps_check_diag = -1.0_dp
109#if defined(__CUSOLVERMP)
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
158 LOGICAL,
SAVE :: initialized = .false.
160 fallback_applied = .false.
162 IF (diag_lib ==
"ScaLAPACK")
THEN
164 ELSE IF (diag_lib ==
"ELPA")
THEN
171 fallback_applied = .true.
173 ELSE IF (diag_lib ==
"cuSOLVER")
THEN
175 ELSE IF (diag_lib ==
"DLAF")
THEN
179 cpabort(
"ERROR in diag_init: CP2K was not compiled with DLA-Future support")
182 cpabort(
"ERROR in diag_init: Initialization of unknown diagonalization library requested")
198 dlaf_neigvec_min = dlaf_neigvec_min_input
200 mark_used(dlaf_neigvec_min_input)
203 elpa_neigvec_min = elpa_neigvec_min_input
204 eps_check_diag = eps_check_diag_input
237 TYPE(
cp_fm_type),
INTENT(IN) :: matrix, eigenvectors
238 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: eigenvalues
239 INTEGER,
INTENT(OUT),
OPTIONAL :: info
241 CHARACTER(LEN=*),
PARAMETER :: routinen =
'choose_eigv_solver'
246 IF (
PRESENT(info)) info = 0
249 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
252 IF (matrix%matrix_struct%nrow_global < elpa_neigvec_min)
THEN
254 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
260 IF (matrix%matrix_struct%nrow_global < 64)
THEN
262 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
269 IF (matrix%matrix_struct%nrow_global < dlaf_neigvec_min)
THEN
271 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
278 cpabort(
"ERROR in "//routinen//
": Invalid diagonalization type requested")
281 CALL check_diag(matrix, eigenvectors, nvec=
SIZE(eigenvalues))
291 SUBROUTINE check_diag(matrix, eigenvectors, nvec)
293 TYPE(
cp_fm_type),
INTENT(IN) :: matrix, eigenvectors
294 INTEGER,
INTENT(IN) :: nvec
296 CHARACTER(LEN=*),
PARAMETER :: routinen =
'check_diag'
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)
304 INTEGER :: il, jl, ipcol, iprow, &
305 mypcol, myprow, npcol, nprow
306 INTEGER,
DIMENSION(9) :: desca
309 CALL timeset(routinen, handle)
313#if defined(__CHECK_DIAG)
314 check_eigenvectors = .true.
315 IF (eps_check_diag >= 0.0_dp)
THEN
316 eps_warning = eps_check_diag
319 IF (eps_check_diag >= 0.0_dp)
THEN
320 check_eigenvectors = .true.
321 eps_warning = eps_check_diag
323 check_eigenvectors = .false.
326 eps_abort = 10.0_dp*eps_warning
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
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
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
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
370 IF (eps > eps_warning)
THEN
372 diag_type_name =
"SYEVD"
374 diag_type_name =
"ELPA"
376 diag_type_name =
"CUSOLVER"
378 diag_type_name =
"DLAF"
380 cpabort(
"Unknown diag_type")
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")
389 cpwarn(
"Check of matrix diagonalization failed in routine "//routinen)
394 CALL timestop(handle)
396 END SUBROUTINE check_diag
413 TYPE(
cp_fm_type),
INTENT(IN) :: matrix, eigenvectors
414 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: eigenvalues
415 INTEGER,
INTENT(OUT),
OPTIONAL :: info
417 CHARACTER(LEN=*),
PARAMETER :: routinen =
'cp_fm_syevd'
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
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
431 CALL timeset(routinen, handle)
435 n = matrix%matrix_struct%nrow_global
438#if defined(__parallel)
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)
450 CALL cp_fm_syevd_base(matrix_new, eigenvectors_new, eig)
459 m => matrix%local_data
468 CALL dsyevd(
'V',
'U', n, m(1, 1), nl, eig(1), work(1), lwork, iwork(1), liwork, myinfo)
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))
475 cpabort(trim(message))
482 ALLOCATE (work(lwork))
487 ALLOCATE (iwork(liwork))
490 CALL dsyevd(
'V',
'U', n, m(1, 1), nl, eig(1), work(1), lwork, iwork(1), liwork, myinfo)
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))
497 cpabort(trim(message))
507 IF (
PRESENT(info)) info = myinfo
509 nmo =
SIZE(eigenvalues, 1)
511 eigenvalues(1:n) = eig(1:n)
513 eigenvalues(1:nmo) = eig(1:nmo)
518 CALL check_diag(matrix, eigenvectors, n)
520 CALL timestop(handle)
531 SUBROUTINE cp_fm_syevd_base(matrix, eigenvectors, eigenvalues, info)
533 TYPE(
cp_fm_type),
INTENT(IN) :: matrix, eigenvectors
534 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: eigenvalues
535 INTEGER,
INTENT(OUT),
OPTIONAL :: info
537 CHARACTER(LEN=*),
PARAMETER :: routinen =
'cp_fm_syevd_base'
539 CHARACTER(LEN=2*default_string_length) :: message
540 INTEGER :: handle, myinfo
541#if defined(__parallel)
543 INTEGER :: liwork, lwork, &
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
554 CALL timeset(routinen, handle)
558#if defined(__parallel)
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(:)
567 v => eigenvectors%local_data
568 descv(:) = eigenvectors%matrix_struct%descriptor(:)
570 liwork = 7*n + 8*context%num_pe(2) + 2
571 ALLOCATE (iwork(liwork))
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)
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))
585 cpabort(trim(message))
595 lwork = nint(work(1) + 100000)
598 ALLOCATE (work(lwork))
602#if defined (__HAS_IEEE_EXCEPTIONS)
603 CALL ieee_get_halting_mode(ieee_all, halt)
604 CALL ieee_set_halting_mode(ieee_all, .false.)
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)
610#if defined (__HAS_IEEE_EXCEPTIONS)
611 CALL ieee_set_halting_mode(ieee_all, halt)
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))
618 cpabort(trim(message))
622 IF (
PRESENT(info)) info = myinfo
628 mark_used(eigenvectors)
629 mark_used(eigenvalues)
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))
637 CALL timestop(handle)
639 END SUBROUTINE cp_fm_syevd_base
655 SUBROUTINE cp_fm_syevx(matrix, eigenvectors, eigenvalues, neig, work_syevx)
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
665 CHARACTER(LEN=*),
PARAMETER :: routinen =
"cp_fm_syevx"
667#if defined(__parallel)
668 REAL(kind=
dp),
PARAMETER :: orfac = -1.0_dp
670 REAL(kind=
dp),
PARAMETER :: vl = 0.0_dp, &
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, &
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
687 REAL(kind=
dp),
EXTERNAL :: dlamch
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
700 INTEGER,
EXTERNAL :: ilaenv
704 n = matrix%matrix_struct%nrow_global
706 IF (
PRESENT(neig)) neig_local = neig
707 IF (neig_local == 0)
RETURN
709 CALL timeset(routinen, handle)
711 needs_evecs =
PRESENT(eigenvectors)
714 ionode = logger%para_env%is_source()
715 n = matrix%matrix_struct%nrow_global
718 work_syevx_local = 1.0_dp
719 IF (
PRESENT(work_syevx)) work_syevx_local = work_syevx
722 IF (needs_evecs)
THEN
729 abstol = 2.0_dp*dlamch(
"S")
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)
738 eigenvalues(:) = 0.0_dp
739#if defined(__parallel)
741 IF (matrix%matrix_struct%nrow_block /= matrix%matrix_struct%ncol_block)
THEN
742 cpabort(
"ERROR in "//routinen//
": Invalid blocksize (no square blocks) found")
745 a => matrix%local_data
746 desca(:) = matrix%matrix_struct%descriptor(:)
748 IF (needs_evecs)
THEN
749 z => eigenvectors%local_data
750 descz(:) = eigenvectors%matrix_struct%descriptor(:)
753 z => matrix%local_data
760 nb = matrix%matrix_struct%nrow_block
762 np0 = numroc(nn, nb, 0, 0, nprow)
763 nq0 = max(numroc(nn, nb, 0, 0, npcol), nb)
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))
769 lwork = 5*n + max(5*nn, nb*(np0 + 1))
771 liwork = 6*max(n, npe + 1, 4)
775 ALLOCATE (iclustr(2*npe))
779 ALLOCATE (iwork(liwork))
780 ALLOCATE (work(lwork))
784#if defined (__HAS_IEEE_EXCEPTIONS)
785 CALL ieee_get_halting_mode(ieee_all, halt)
786 CALL ieee_set_halting_mode(ieee_all, .false.)
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)
792#if defined (__HAS_IEEE_EXCEPTIONS)
793 CALL ieee_set_halting_mode(ieee_all, halt)
801 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T12,1X,I10))") &
804 "liwork = ", liwork, &
807 WRITE (unit=output_unit, fmt=
"(/,T3,A,(T12,6(1X,I10)))") &
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)))") &
815 cpabort(
"ERROR in PDSYEVX (ScaLAPACK)")
828 a => matrix%local_data
829 IF (needs_evecs)
THEN
830 z => eigenvectors%local_data
833 z => matrix%local_data
838 nb = max(ilaenv(1,
"DSYTRD",
"U", n, -1, -1, -1), &
839 ilaenv(1,
"DORMTR",
"U", n, -1, -1, -1))
841 lwork = max((nb + 3)*n, 8*n) + n
846 ALLOCATE (iwork(liwork))
847 ALLOCATE (work(lwork))
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)
859 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T12,1X,I10))") &
862 WRITE (unit=output_unit, fmt=
"(/,T3,A,(T12,6(1X,I10)))") &
865 cpabort(
"Error in DSYEVX (ScaLAPACK)")
875 eigenvalues(1:neig_local) = w(1:neig_local)
878 IF (needs_evecs)
CALL check_diag(matrix, eigenvectors, neig_local)
880 CALL timestop(handle)
893 SUBROUTINE cp_fm_svd(matrix_a, matrix_eigvl, matrix_eigvr_t, eigval, info)
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
901 CHARACTER(LEN=*),
PARAMETER :: routinen =
'cp_fm_svd'
903 CHARACTER(LEN=2*default_string_length) :: message
904 INTEGER :: handle, n, m, myinfo, lwork
905 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: a
907 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: work
909#if defined(__parallel)
910 INTEGER,
DIMENSION(9) :: desca, descu, descvt
912 CALL timeset(routinen, handle)
915 matrix_struct=matrix_a%matrix_struct, &
918 a => matrix_lu%local_data
919 m = matrix_lu%matrix_struct%nrow_global
920 n = matrix_lu%matrix_struct%ncol_global
929#if defined(__parallel)
931 desca(:) = matrix_lu%matrix_struct%descriptor(:)
932 descu(:) = matrix_eigvl%matrix_struct%descriptor(:)
933 descvt(:) = matrix_eigvr_t%matrix_struct%descriptor(:)
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)
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))
944 cpabort(trim(message))
950 ALLOCATE (work(lwork))
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)
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))
960 cpabort(trim(message))
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)
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))
973 cpabort(trim(message))
980 ALLOCATE (work(lwork))
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)
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))
991 cpabort(trim(message))
1000 IF (
PRESENT(info)) info = myinfo
1002 CALL timestop(handle)
1015 SUBROUTINE cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
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), &
1031 CHARACTER(LEN=*),
PARAMETER :: routinen =
'cp_fm_power'
1033 INTEGER :: handle, icol_global, &
1035 ncol_global, npcol, nprow, &
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
1043#if defined(__parallel)
1044 INTEGER :: icol_local, ipcol, iprow, irow_global, &
1048 CALL timeset(routinen, handle)
1050 my_verbose = .false.
1051 IF (
PRESENT(verbose)) my_verbose = verbose
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)
1061 nrow_global = matrix%matrix_struct%nrow_global
1062 ncol_global = matrix%matrix_struct%ncol_global
1064 ALLOCATE (eigenvalues(ncol_global))
1066 eigenvalues(:) = 0.0_dp
1072 IF (
PRESENT(eigvals))
THEN
1073 eigvals(1) = eigenvalues(1)
1074 eigvals(2) = eigenvalues(ncol_global)
1077#if defined(__parallel)
1078 eigenvectors => work%local_data
1082 DO icol_global = 1, ncol_global
1084 IF (eigenvalues(icol_global) < threshold)
THEN
1086 n_dependent = n_dependent + 1
1088 ipcol = work%matrix_struct%g2p_col(icol_global)
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
1103 f = eigenvalues(icol_global)**p
1105 ipcol = work%matrix_struct%g2p_col(icol_global)
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)
1125 eigenvectors => work%local_data
1129 DO icol_global = 1, ncol_global
1131 IF (eigenvalues(icol_global) < threshold)
THEN
1133 n_dependent = n_dependent + 1
1134 eigenvectors(1:nrow_global, icol_global) = 0.0_dp
1138 f = eigenvalues(icol_global)**p
1139 eigenvectors(1:nrow_global, icol_global) = &
1140 f*eigenvectors(1:nrow_global, icol_global)
1147 CALL cp_fm_syrk(
"U",
"N", ncol_global, 1.0_dp, work, 1, 1, 0.0_dp, matrix)
1151 IF (matrix%matrix_struct%para_env%is_source() .AND. my_verbose)
THEN
1152 condition_number = abs(eigenvalues(ncol_global)/eigenvalues(1))
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
1159 "WARNING: matrix has a negative eigenvalue, tighten EPS_DEFAULT"
1161 IF (condition_number > 1.0e12_dp)
THEN
1163 "WARNING: high condition number => possibly ill-conditioned matrix"
1167 DEALLOCATE (eigenvalues)
1169 CALL timestop(handle)
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
1198 CHARACTER(len=*),
PARAMETER :: routinen =
'cp_fm_block_jacobi'
1201 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: a, ev
1203 REAL(kind=
dp) :: tan_theta, tau, c, s
1205 REAL(kind=
dp),
DIMENSION(:),
ALLOCATABLE :: c_ip
1207#if defined(__parallel)
1210 INTEGER :: myprow, mypcol, nprow, npcol, block_dim_row, block_dim_col, &
1211 info, ev_row_block_size, iam, mynumrows, mype, npe, &
1213 REAL(kind=
dp),
DIMENSION(:, :),
ALLOCATABLE :: a_loc, ev_loc
1214 INTEGER,
DIMENSION(9) :: desca, descz, desc_a_block, &
1218 INTEGER,
EXTERNAL :: numroc
1223 CALL timeset(routinen, handle)
1225#if defined(__parallel)
1226 context => matrix%matrix_struct%context
1227 allgrp = matrix%matrix_struct%para_env
1229 myprow = context%mepos(1)
1230 mypcol = context%mepos(2)
1231 nprow = context%num_pe(1)
1232 npcol = context%num_pe(2)
1234 n = matrix%matrix_struct%nrow_global
1236 a => matrix%local_data
1237 desca(:) = matrix%matrix_struct%descriptor(:)
1238 ev => eigenvectors%local_data
1239 descz(:) = eigenvectors%matrix_struct%descriptor(:)
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))
1249 mype = matrix%matrix_struct%para_env%mepos
1250 npe = matrix%matrix_struct%para_env%num_pe
1252 CALL ictxt_loc%gridinit(matrix%matrix_struct%para_env,
'Row-major', nprow*npcol, 1)
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)
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())
1260 CALL allgrp%bcast(a_loc, 0)
1267 ev_row_block_size = n/(nprow*npcol)
1268 mynumrows = numroc(n, ev_row_block_size, iam, 0, nprow*npcol)
1270 ALLOCATE (ev_loc(mynumrows, n), c_ip(mynumrows))
1272 CALL descinit(desc_ev_loc, n, n, ev_row_block_size, n, 0, 0, ictxt_loc%get_handle(), &
1275 CALL pdgemr2d(n, n, ev, 1, 1, descz, ev_loc, 1, 1, desc_ev_loc, context%get_handle())
1281 DO q = start_sec_block, n
1283 DO p = 1, (start_sec_block - 1)
1285 IF (abs(a_loc(p, q_loc)) > thresh)
THEN
1287 tau = (eigval(q) - eigval(p))/(2.0_dp*a_loc(p, q_loc))
1289 tan_theta = sign(1.0_dp, tau)/(abs(tau) + sqrt(1.0_dp + tau*tau))
1292 c = 1.0_dp/sqrt(1.0_dp + tan_theta*tan_theta)
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)
1312 CALL pdgemr2d(n, n, ev_loc, 1, 1, desc_ev_loc, ev, 1, 1, descz, context%get_handle())
1315 DEALLOCATE (a_loc, ev_loc, c_ip)
1317 CALL ictxt_loc%gridexit()
1321 n = matrix%matrix_struct%nrow_global
1325 a => matrix%local_data
1326 ev => eigenvectors%local_data
1333 DO q = start_sec_block, n
1334 DO p = 1, (start_sec_block - 1)
1336 IF (abs(a(p, q)) > thresh)
THEN
1338 tau = (eigval(q) - eigval(p))/(2.0_dp*a(p, q))
1340 tan_theta = sign(1.0_dp, tau)/(abs(tau) + sqrt(1.0_dp + tau*tau))
1343 c = 1.0_dp/sqrt(1.0_dp + tan_theta*tan_theta)
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)
1368 CALL timestop(handle)
1381 SUBROUTINE cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
1383 TYPE(
cp_fm_type),
INTENT(IN) :: amatrix, bmatrix, eigenvectors
1384 REAL(kind=
dp),
DIMENSION(:) :: eigenvalues
1387 CHARACTER(len=*),
PARAMETER :: routinen =
'cp_fm_geeig'
1389 INTEGER :: handle, nao, nmo
1391 CALL timeset(routinen, handle)
1394 nmo =
SIZE(eigenvalues)
1414 eigenvalues=eigenvalues)
1420 CALL timestop(handle)
1436 TYPE(
cp_fm_type),
INTENT(IN) :: amatrix, bmatrix, eigenvectors
1437 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: eigenvalues
1439 REAL(kind=
dp),
INTENT(IN) :: epseig
1441 CHARACTER(len=*),
PARAMETER :: routinen =
'cp_fm_geeig_canon'
1443 INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
1445 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: evals
1447 CALL timeset(routinen, handle)
1451 nmo =
SIZE(eigenvalues)
1452 ALLOCATE (evals(nao))
1457 evals(:) = -evals(:)
1460 IF (evals(i) < epseig)
THEN
1471 CALL cp_fm_to_fm(work, eigenvectors, ncol, nc + 1, nc + 1)
1474 DO icol = nc + 1, nao
1480 evals(nc + 1:nao) = 1.0_dp
1483 evals(:) = 1.0_dp/sqrt(evals(:))
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)
1490 DO icol = nc + 1, nao
1495 CALL choose_eigv_solver(matrix=amatrix, eigenvectors=bmatrix, eigenvalues=eigenvalues)
1498 CALL cp_fm_gemm(
"N",
"N", nao, nx, nc, 1.0_dp, work, bmatrix, 0.0_dp, eigenvectors)
1502 CALL timestop(handle)
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,...
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...
subroutine, public cp_fm_block_jacobi(matrix, eigenvectors, eigval, thresh, start_sec_block)
...
real(kind=dp), parameter, public eps_check_diag_default
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
integer, parameter, public fm_diag_type_cusolver
subroutine, public cp_fm_svd(matrix_a, matrix_eigvl, matrix_eigvr_t, eigval, info)
decomposes a quadratic matrix into its singular value decomposition
integer, parameter, public fm_diag_type_dlaf
integer, parameter, public fm_diag_type_scalapack
subroutine, public cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
Computes all eigenvalues and vectors of a real symmetric matrix significantly faster than syevx,...
subroutine, public diag_finalize()
Finalize the diagonalization library.
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.
integer, parameter, public fm_diag_type_default
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...
integer, save, public diag_type
integer, parameter, public fm_diag_type_elpa
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....
subroutine, public cp_fm_diag_dlaf(matrix, eigenvectors, eigenvalues)
...
subroutine, public cp_fm_diag_gen_dlaf(a_matrix, b_matrix, eigenvectors, eigenvalues)
...
subroutine, public set_elpa_print(flag)
Sets a flag that determines if additional information about the ELPA diagonalization should be printe...
subroutine, public set_elpa_kernel(requested_kernel)
Sets the active ELPA kernel.
subroutine, public cp_fm_diag_elpa(matrix, eigenvectors, eigenvalues)
Driver routine to diagonalize a FM matrix with the ELPA library.
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 = ....
subroutine, public initialize_elpa_library()
Initialize the ELPA library.
subroutine, public finalize_elpa_library()
Finalize the ELPA library.
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
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.
integer, parameter, public dp
integer, parameter, public default_string_length
Machine interface based on Fortran 2003 and POSIX.
integer, parameter, public default_output_unit
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
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
type of a logger, at the moment it contains just a print level starting at which level it should be l...