74#include "./base/base_uses.f90"
80 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_kpoint_operators'
101 INTEGER,
PARAMETER,
PRIVATE :: kpoint_part_re = 1, &
106 INTEGER,
PARAMETER,
PRIVATE :: kpoint_slot_idle = 0, &
107 kpoint_slot_in_flight = 1, &
108 kpoint_slot_local = 2
112 TYPE kpoint_operator_rs_route_type
115 INTEGER :: nimage = 0, nvalue = 0, nrow = 0, ncol = 0
116 INTEGER(KIND=int_8) :: storage = 0
117 INTEGER,
ALLOCATABLE :: row_size(:), col_size(:), entry(:, :), first(:), tile(:, :, :)
118 INTEGER,
ALLOCATABLE :: send_entry(:), send_first(:), scount(:, :), rcount(:, :), &
119 sdisp(:, :), rdisp(:, :)
120 REAL(kind=
dp),
ALLOCATABLE :: send_values(:)
121 END TYPE kpoint_operator_rs_route_type
125 TYPE kpoint_operator_block_type
126 REAL(kind=
dp),
POINTER ::
data(:, :) => null()
127 END TYPE kpoint_operator_block_type
129 TYPE kpoint_operator_rs_type
130 INTEGER :: matrix_row = 0
134 REAL(kind=
dp),
CONTIGUOUS,
POINTER :: values(:) => null()
135 TYPE(kpoint_operator_block_type),
ALLOCATABLE :: blocks(:)
136 LOGICAL,
ALLOCATABLE :: nonzero(:)
137 END TYPE kpoint_operator_rs_type
141 TYPE kpoint_operator_rs_batch_type
143 COMPLEX(KIND=dp),
ALLOCATABLE :: data(:, :, :)
144 END TYPE kpoint_operator_rs_batch_type
148 INTEGER(KIND=int_8),
PARAMETER :: rs_batch_bytes = 32_int_8*1024_int_8**2
170 POINTER :: work_sab => null()
180 TYPE(kpoint_operator_rs_type),
ALLOCATABLE :: rs(:, :)
181 TYPE(kpoint_operator_rs_route_type),
ALLOCATABLE :: rs_route
182 TYPE(kpoint_operator_rs_batch_type),
ALLOCATABLE :: rs_batch(:, :)
189 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: status
207 INTEGER,
INTENT(IN) :: matrix_row, ik
209 LOGICAL,
INTENT(IN) :: use_grid
211 CHARACTER(LEN=*),
PARAMETER :: routineN =
'kpoint_operator_densify'
214 REAL(KIND=
dp),
DIMENSION(:, :),
POINTER :: xkp
216 CALL timeset(routinen, handle)
219 xkp => ctx%kpoints%xkp
220 IF (ctx%kpoints%use_real_wfn)
THEN
222 cpassert(
PRESENT(grid))
226 CALL rskp_transform(rmatrix=ctx%rmatrix, rsmat=rsmat, ispin=matrix_row, &
227 xkp=xkp(:, ik), cell_to_index=ctx%kpoints%cell_to_index, &
228 sab_nl=ctx%kpoints%sab_nl)
230 ELSE IF (use_grid)
THEN
231 cpassert(
PRESENT(grid))
236 CALL rskp_transform(rmatrix=ctx%rmatrix, cmatrix=ctx%cmatrix, rsmat=rsmat, ispin=matrix_row, &
237 xkp=xkp(:, ik), cell_to_index=ctx%kpoints%cell_to_index, &
238 sab_nl=ctx%kpoints%sab_nl)
240 CALL copy_dbcsr_to_fm(ctx%rmatrix, ctx%fmwork(kpoint_part_re), ctx%copy_plan(kpoint_part_re))
241 IF (.NOT. ctx%kpoints%use_real_wfn)
THEN
242 CALL copy_dbcsr_to_fm(ctx%cmatrix, ctx%fmwork(kpoint_part_im), ctx%copy_plan(kpoint_part_im))
245 CALL timestop(handle)
262 SUBROUTINE kpoint_operator_launch(kpoints, ik, fmwork, fmlocal, info, target_re, target_im)
265 INTEGER,
INTENT(IN) :: ik
266 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fmwork
269 TYPE(
cp_fm_type),
INTENT(IN),
OPTIONAL :: target_re, target_im
271 CHARACTER(LEN=*),
PARAMETER :: routineN =
'kpoint_operator_launch'
277 CALL timeset(routinen, handle)
278 cpassert(.NOT.
PRESENT(target_im) .OR.
PRESENT(target_re))
279 IF (.NOT. kpoints%use_real_wfn)
THEN
280 cpassert(
PRESENT(target_re) .EQV.
PRESENT(target_im))
283 para_env => kpoints%blacs_env_all%para_env
284 IF (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
THEN
285 IF (
PRESENT(target_re))
THEN
287 IF (.NOT. kpoints%use_real_wfn)
THEN
292 IF (.NOT. kpoints%use_real_wfn)
THEN
298 IF (.NOT. kpoints%use_real_wfn)
THEN
303 CALL timestop(handle)
305 END SUBROUTINE kpoint_operator_launch
321 SUBROUTINE kpoint_operator_receive(kpoints, fmlocal, ik, info, cmat, fm_re, fm_im, mat_re, mat_im)
324 INTEGER,
INTENT(IN) :: ik
328 TYPE(
dbcsr_type),
INTENT(INOUT),
OPTIONAL :: mat_re, mat_im
332 IF (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
THEN
333 IF (
PRESENT(cmat))
THEN
336 IF (.NOT. kpoints%use_real_wfn)
THEN
340 ELSE IF (
PRESENT(fm_re))
THEN
342 IF (kpoints%use_real_wfn)
THEN
343 IF (
PRESENT(fm_im)) fm_im%local_data(:, :) = 0.0_dp
350 IF (kpoints%use_real_wfn)
THEN
351 fmlocal%local_data(:, :) = 0.0_dp
358 DO part = 1, merge(1, 2, kpoints%use_real_wfn)
361 END SUBROUTINE kpoint_operator_receive
382 INTEGER,
INTENT(IN) :: nspin
384 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_operator_context_create'
388 CALL timeset(routinen, handle)
392 cpassert(
ASSOCIATED(kpoints))
393 ctx%kpoints => kpoints
395 cpassert(
ASSOCIATED(ctx%ao_ao_fmstruct))
398 IF (.NOT.
ASSOCIATED(ctx%ao_ao_fmstruct%para_env, kpoints%blacs_env_all%para_env))
THEN
399 CALL cp_abort(__location__,
"Operator context requires the global parallel environment")
402 CALL timestop(handle)
418 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_operator_context_release'
420 INTEGER :: handle, ispin, iwork, slot
422 CALL timeset(routinen, handle)
426 IF (
ALLOCATED(ctx%status))
THEN
427 cpassert(all(ctx%status == kpoint_slot_idle))
432 NULLIFY (ctx%rmatrix, ctx%cmatrix, ctx%work_sab)
433 IF (
ALLOCATED(ctx%copy_plan))
DEALLOCATE (ctx%copy_plan)
434 IF (
ALLOCATED(ctx%rs))
THEN
435 DO slot = 1, kpoint_num_slots
436 DO ispin = 1, ctx%nspin
437 CALL kpoint_operator_rs_clear(ctx%rs(ispin, slot))
440 DEALLOCATE (ctx%rs, ctx%rs_batch)
442 IF (
ALLOCATED(ctx%rs_route))
THEN
444 DEALLOCATE (ctx%rs_route)
446 IF (
ALLOCATED(ctx%fmwork))
THEN
447 DO iwork = 1,
SIZE(ctx%fmwork)
450 DEALLOCATE (ctx%fmwork)
452 IF (
ALLOCATED(ctx%fmimag))
THEN
454 DEALLOCATE (ctx%fmimag)
456 IF (
ALLOCATED(ctx%fmlocal))
THEN
458 DEALLOCATE (ctx%fmlocal)
460 IF (
ALLOCATED(ctx%info))
DEALLOCATE (ctx%info)
461 IF (
ALLOCATED(ctx%status))
DEALLOCATE (ctx%status)
462 NULLIFY (ctx%kpoints, ctx%ao_ao_fmstruct, ctx%group_fm_pools)
464 CALL timestop(handle)
484 SUBROUTINE kpoint_operator_ensure(ctx, rsmat)
489 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_operator_ensure'
493 CALL timeset(routinen, handle)
498 IF (
ASSOCIATED(ctx%work_sab))
THEN
499 cpassert(
ASSOCIATED(ctx%work_sab, ctx%kpoints%sab_nl))
501 IF (.NOT.
ASSOCIATED(ctx%rmatrix))
THEN
502 ALLOCATE (ctx%rmatrix)
503 CALL dbcsr_create(ctx%rmatrix, template=rsmat(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
505 IF (.NOT. ctx%kpoints%use_real_wfn)
THEN
506 ALLOCATE (ctx%cmatrix)
507 CALL dbcsr_create(ctx%cmatrix, template=rsmat(1, 1)%matrix, matrix_type=dbcsr_type_antisymmetric)
510 ctx%work_sab => ctx%kpoints%sab_nl
513 IF (.NOT.
ALLOCATED(ctx%fmwork))
THEN
514 ALLOCATE (ctx%fmwork(merge(1, 2, ctx%kpoints%use_real_wfn)))
515 ALLOCATE (ctx%copy_plan(
SIZE(ctx%fmwork)))
516 CALL cp_fm_create(ctx%fmwork(kpoint_part_re), ctx%ao_ao_fmstruct, &
517 name=
"KPOINT-OPERATOR-WORK-RE")
518 IF (.NOT. ctx%kpoints%use_real_wfn)
THEN
519 CALL cp_fm_create(ctx%fmwork(kpoint_part_im), ctx%ao_ao_fmstruct, &
520 name=
"KPOINT-OPERATOR-WORK-IM")
524 CALL kpoint_operator_ensure_group(ctx)
526 CALL timestop(handle)
528 END SUBROUTINE kpoint_operator_ensure
534 SUBROUTINE kpoint_operator_ensure_group(ctx)
538 IF (.NOT.
ALLOCATED(ctx%fmlocal))
THEN
539 CALL mpools_get(ctx%kpoints%mpools, ao_ao_fm_pools=ctx%group_fm_pools)
540 ALLOCATE (ctx%fmlocal)
544 END SUBROUTINE kpoint_operator_ensure_group
565 INTEGER,
INTENT(IN) :: slot
566 LOGICAL,
INTENT(OUT),
OPTIONAL :: prepared
567 INTEGER,
INTENT(IN),
OPTIONAL :: matrix_row
568 TYPE(
cp_fm_type),
INTENT(IN),
OPTIONAL :: target
570 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_operator_prepare'
571 INTEGER(KIND=int_8),
PARAMETER :: max_bytes = 64_int_8*1024_int_8**2
573 INTEGER :: compatible, handle, ispin, jslot, nrow, &
575 INTEGER(KIND=int_8) :: retained, storage
580 CALL timeset(routinen, handle)
581 IF (
PRESENT(prepared)) prepared = .false.
582 cpassert(
ASSOCIATED(ctx%kpoints))
584 cpassert(
ASSOCIATED(ctx%kpoints%sab_nl))
585 IF (
ASSOCIATED(ctx%work_sab))
THEN
586 cpassert(
ASSOCIATED(ctx%work_sab, ctx%kpoints%sab_nl))
588 IF (
ALLOCATED(ctx%status))
THEN
589 cpassert(all(ctx%status == kpoint_slot_idle))
591 IF (.NOT.
ALLOCATED(ctx%rs))
THEN
592 ALLOCATE (ctx%rs(ctx%nspin, kpoint_num_slots), ctx%rs_batch(ctx%nspin, kpoint_num_slots))
594 DO ispin = 1, ctx%nspin
595 CALL kpoint_operator_rs_clear(ctx%rs(ispin, slot))
596 IF (
ALLOCATED(ctx%rs_batch(ispin, slot)%data))
DEALLOCATE (ctx%rs_batch(ispin, slot)%data)
597 ctx%rs_batch(ispin, slot)%first = 0
600 IF (
PRESENT(matrix_row))
THEN
601 cpassert(matrix_row >= 1 .AND. matrix_row <=
SIZE(rsmat, 1))
604 cpassert(
SIZE(rsmat, 1) == ctx%nspin)
606 para_env => ctx%kpoints%blacs_env_all%para_env
607 CALL kpoint_operator_ensure_group(ctx)
608 fms => ctx%fmlocal%matrix_struct
609 IF (
PRESENT(
TARGET)) fms =>
target%matrix_struct
610 cpassert(
ASSOCIATED(fms))
613 compatible = merge(1, 0,
ASSOCIATED(fms%para_env, ctx%kpoints%para_env_kp))
614 CALL para_env%min(compatible)
615 IF (compatible == 0)
THEN
616 CALL timestop(handle)
619 IF (.NOT.
ALLOCATED(ctx%rs_route))
THEN
620 ALLOCATE (ctx%rs_route)
621 CALL kpoint_operator_rs_plan_create(ctx%kpoints, rsmat, fms, max_bytes, ctx%rs_route, ready)
622 IF (.NOT. ready)
THEN
623 DEALLOCATE (ctx%rs_route)
624 CALL timestop(handle)
627 ctx%work_sab => ctx%kpoints%sab_nl
632 IF (
PRESENT(matrix_row)) row = matrix_row
633 IF (.NOT. kpoint_operator_rs_compatible(ctx%rs_route, rsmat, row)) compatible = 0
635 CALL para_env%min(compatible)
639 DO ispin = 1, ctx%nspin
640 IF (.NOT.
ALLOCATED(ctx%rs(ispin, jslot)%blocks)) cycle
641 retained = retained + 128_int_8*
SIZE(ctx%rs(ispin, jslot)%blocks, kind=
int_8)
642 IF (
ASSOCIATED(ctx%rs(ispin, jslot)%values))
THEN
643 retained = retained + 8_int_8*
SIZE(ctx%rs(ispin, jslot)%values, kind=
int_8)
647 storage = ctx%rs_route%storage + retained + int(nrow,
int_8)* &
648 (merge(0_int_8, 8_int_8*ctx%rs_route%nvalue, para_env%num_pe == 1) + &
649 128_int_8*
SIZE(ctx%rs_route%entry, 2))
650 CALL para_env%max(storage)
651 ready = compatible == 1 .AND. storage <= max_bytes
655 IF (
PRESENT(matrix_row)) row = matrix_row
656 CALL kpoint_operator_rs_capture(rsmat, row, ctx%rs_route, ctx%kpoints, ctx%rs(ispin, slot))
659 IF (
PRESENT(prepared)) prepared = ready
660 CALL timestop(handle)
668 SUBROUTINE kpoint_operator_rs_clear(snapshot)
670 TYPE(kpoint_operator_rs_type),
INTENT(INOUT) :: snapshot
672 IF (
ALLOCATED(snapshot%blocks))
DEALLOCATE (snapshot%blocks, snapshot%nonzero)
673 IF (
ASSOCIATED(snapshot%values))
DEALLOCATE (snapshot%values)
674 NULLIFY (snapshot%source)
675 snapshot%matrix_row = 0
677 END SUBROUTINE kpoint_operator_rs_clear
688 FUNCTION kpoint_operator_rs_match(ctx, source, ispin, slot, matrix_row)
RESULT(irs)
692 INTEGER,
INTENT(IN) :: ispin, slot, matrix_row
698 IF (.NOT.
ALLOCATED(ctx%rs))
RETURN
700 IF (.NOT.
ASSOCIATED(ctx%rs(i, slot)%source)) i = 1
701 IF (
ASSOCIATED(ctx%rs(i, slot)%source, source) .AND. ctx%rs(i, slot)%matrix_row == matrix_row) irs = i
703 END FUNCTION kpoint_operator_rs_match
717 SUBROUTINE kpoint_operator_rs_plan_create(kpoints, source, matrix_struct, max_bytes, route, ready)
722 INTEGER(KIND=int_8),
INTENT(IN) :: max_bytes
723 TYPE(kpoint_operator_rs_route_type),
INTENT(OUT) :: route
724 LOGICAL,
INTENT(OUT) :: ready
726 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_operator_rs_plan_create'
728 INTEGER :: compatible, dest, handle, i, iblock, &
729 nblock, ncol, nrow, p, pos, rank
730 INTEGER(KIND=int_8) :: nr(2), ns(2), storage
731 INTEGER(KIND=int_8),
ALLOCATABLE :: recv_size(:, :), send_size(:, :)
732 INTEGER,
ALLOCATABLE :: identities(:, :), owner(:), &
733 position(:, :), send_metadata(:, :)
734 INTEGER,
POINTER :: col_size(:), row_size(:)
735 LOGICAL,
ALLOCATABLE :: needed(:)
738 CALL timeset(routinen, handle)
740 world => kpoints%blacs_env_all%para_env
741 inter => kpoints%para_env_inter_kp
742 CALL dbcsr_get_info(source(1, 1)%matrix, row_blk_size=row_size, col_blk_size=col_size)
744 SIZE(source, 2), block_template=source(1, 1)%matrix)
745 compatible = merge(1, 0, route%traversal%symmetric)
746 IF (sum(row_size) /= matrix_struct%nrow_global .OR. &
747 sum(col_size) /= matrix_struct%ncol_global) compatible = 0
748 IF (
SIZE(row_size) ==
SIZE(col_size))
THEN
749 IF (any(row_size /= col_size)) compatible = 0
753 CALL world%min(compatible)
754 IF (compatible == 0)
THEN
755 CALL timestop(handle)
758 route%nimage =
SIZE(source, 2)
759 ALLOCATE (route%row_size(
SIZE(row_size)), route%col_size(
SIZE(col_size)))
760 route%row_size(:) = row_size
761 route%col_size(:) = col_size
762 ALLOCATE (identities(2, 0:world%num_pe - 1), owner(0:kpoints%para_env_kp%num_pe - 1), &
763 needed(0:kpoints%para_env_kp%num_pe - 1))
764 CALL world%allgather([kpoints%para_env_kp%mepos, inter%mepos], identities)
766 DO rank = 0, world%num_pe - 1
767 IF (identities(2, rank) == 0) owner(identities(1, rank)) = rank
769 cpassert(all(owner >= 0))
770 ALLOCATE (send_size(2, 0:world%num_pe - 1), recv_size(2, 0:world%num_pe - 1))
771 send_size(:, :) = 0_int_8
772 DO i = 1, route%traversal%nentry
773 nrow = row_size(route%traversal%row(i))
774 ncol = col_size(route%traversal%col(i))
775 CALL kpoint_operator_rs_owners(route, matrix_struct, i, needed)
776 DO p = 0,
SIZE(owner) - 1
777 IF (.NOT. needed(p)) cycle
779 send_size(1, dest) = send_size(1, dest) + 1_int_8
780 send_size(2, dest) = send_size(2, dest) + int(nrow,
int_8)*ncol
783 CALL world%alltoall(send_size, recv_size, 2)
786 ns = sum(send_size, dim=2)
787 nr = sum(recv_size, dim=2)
788 CALL inter%bcast(nr, 0)
791 storage = 40_int_8*route%traversal%nentry + &
792 4_int_8*(
SIZE(row_size) +
SIZE(col_size) + matrix_struct%nrow_global + matrix_struct%ncol_global) + &
793 48_int_8*ns(1) + merge(0_int_8, 8_int_8*ns(2), world%num_pe == 1) + 68_int_8*nr(1) + &
794 128_int_8*world%num_pe + 8_int_8*
SIZE(owner) + 8_int_8
795 route%storage = storage
796 storage = storage + merge(0_int_8, 8_int_8*nr(2), world%num_pe == 1) + 128_int_8*nr(1)
798 IF (maxval(ns) >= int(huge(0),
int_8) .OR. maxval(nr) >= int(huge(0),
int_8))
THEN
799 storage = max(storage, max_bytes + 1_int_8)
801 CALL world%max(storage)
802 IF (storage > max_bytes)
THEN
803 CALL timestop(handle)
808 route%nvalue = int(nr(2))
809 ALLOCATE (route%scount(2, 0:world%num_pe - 1), route%rcount(2, 0:world%num_pe - 1), &
810 route%sdisp(2, 0:world%num_pe - 1), route%rdisp(2, 0:world%num_pe - 1))
811 route%scount(:, :) = int(send_size)
812 route%rcount(:, :) = int(recv_size)
813 IF (
SIZE(owner) == 1)
CALL inter%bcast(route%rcount, 0)
814 route%sdisp(:, 0) = 0
815 route%rdisp(:, 0) = 0
816 DO rank = 1, world%num_pe - 1
817 route%sdisp(:, rank) = route%sdisp(:, rank - 1) + route%scount(:, rank - 1)
818 route%rdisp(:, rank) = route%rdisp(:, rank - 1) + route%rcount(:, rank - 1)
820 ALLOCATE (route%send_entry(int(ns(1))), route%send_first(int(ns(1)) + 1), &
821 route%send_values(merge(0, int(ns(2)), world%num_pe == 1)), send_metadata(8, int(ns(1))), &
822 position(2, 0:world%num_pe - 1), route%entry(8, nblock), &
823 route%first(nblock + 1))
824 position(:, :) = route%sdisp
825 DO i = 1, route%traversal%nentry
826 nrow = row_size(route%traversal%row(i))
827 ncol = col_size(route%traversal%col(i))
828 CALL kpoint_operator_rs_owners(route, matrix_struct, i, needed)
829 DO p = 0,
SIZE(owner) - 1
830 IF (.NOT. needed(p)) cycle
832 iblock = position(1, dest) + 1
833 pos = position(2, dest) + 1
834 route%send_entry(iblock) = i
835 route%send_first(iblock) = pos
836 send_metadata(1:4, iblock) = [route%traversal%row_offset(i), route%traversal%col_offset(i), nrow, ncol]
837 send_metadata(5:7, iblock) = nint(route%traversal%symmetry_sign(i))*route%traversal%cell(:, i)
838 send_metadata(8, iblock) = merge(1, 0, route%traversal%row(i) /= route%traversal%col(i))
839 position(:, dest) = position(:, dest) + [1, nrow*ncol]
842 route%send_first(int(ns(1)) + 1) = int(ns(2)) + 1
843 IF (
SIZE(owner) == 1)
THEN
844 CALL world%allgatherv(send_metadata, route%entry, 8*route%rcount(1, :), 8*route%rdisp(1, :))
846 CALL world%alltoall(send_metadata, 8*route%scount(1, :), 8*route%sdisp(1, :), &
847 route%entry, 8*route%rcount(1, :), 8*route%rdisp(1, :))
848 CALL inter%bcast(route%entry, 0)
852 route%first(i + 1) = route%first(i) + route%entry(3, i)*route%entry(4, i)
854 cpassert(route%first(nblock + 1) == route%nvalue + 1)
855 route%matrix_struct => matrix_struct
857 CALL kpoint_operator_rs_tile_map(route)
859 CALL timestop(handle)
861 END SUBROUTINE kpoint_operator_rs_plan_create
868 SUBROUTINE kpoint_operator_rs_tile_map(route)
869 TYPE(kpoint_operator_rs_route_type),
INTENT(INOUT) :: route
871 INTEGER :: c0, i, nc, nr, r0, rows(0:route%matrix_struct%nrow_global), &
872 cols(0:route%matrix_struct%ncol_global)
876 DO i = 1,
SIZE(route%matrix_struct%row_indices)
877 r0 = route%matrix_struct%row_indices(i)
878 IF (r0 >= 1 .AND. r0 <= ubound(rows, 1)) rows(r0) = 1
880 DO i = 1,
SIZE(route%matrix_struct%col_indices)
881 c0 = route%matrix_struct%col_indices(i)
882 IF (c0 >= 1 .AND. c0 <= ubound(cols, 1)) cols(c0) = 1
884 DO i = 1, ubound(rows, 1)
885 rows(i) = rows(i) + rows(i - 1)
887 DO i = 1, ubound(cols, 1)
888 cols(i) = cols(i) + cols(i - 1)
890 route%nrow = rows(ubound(rows, 1))
891 route%ncol = cols(ubound(cols, 1))
892 ALLOCATE (route%tile(4, 2,
SIZE(route%entry, 2)))
893 DO i = 1,
SIZE(route%entry, 2)
894 r0 = route%entry(1, i)
895 c0 = route%entry(2, i)
896 nr = route%entry(3, i)
897 nc = route%entry(4, i)
898 route%tile(:, 1, i) = [rows(r0 - 1) + 1, rows(r0 + nr - 1), cols(c0 - 1) + 1, cols(c0 + nc - 1)]
899 route%tile(:, 2, i) = [rows(c0 - 1) + 1, rows(c0 + nc - 1), cols(r0 - 1) + 1, cols(r0 + nr - 1)]
901 END SUBROUTINE kpoint_operator_rs_tile_map
910 SUBROUTINE kpoint_operator_rs_owners(route, fms, i, needed)
912 TYPE(kpoint_operator_rs_route_type),
INTENT(IN) :: route
914 INTEGER,
INTENT(IN) :: i
915 LOGICAL,
INTENT(OUT) :: needed(0:)
917 INTEGER :: c0, j, nc, nr, pc, pr, r0
918 LOGICAL :: rows(2, 0:fms%context%num_pe(1) - 1), cols(2, 0:fms%context%num_pe(2) - 1)
920 r0 = route%traversal%row_offset(i)
921 c0 = route%traversal%col_offset(i)
922 nr = route%row_size(route%traversal%row(i))
923 nc = route%col_size(route%traversal%col(i))
926 DO j = r0, r0 + nr - 1
927 rows(1, fms%g2p_row(j)) = .true.
928 cols(2, fms%g2p_col(j)) = .true.
930 DO j = c0, c0 + nc - 1
931 rows(2, fms%g2p_row(j)) = .true.
932 cols(1, fms%g2p_col(j)) = .true.
934 IF (route%traversal%row(i) == route%traversal%col(i)) rows(2, :) = .false.
936 DO pc = 0, fms%context%num_pe(2) - 1
937 DO pr = 0, fms%context%num_pe(1) - 1
938 needed(fms%context%blacs2mpi(pr, pc)) = &
939 (rows(1, pr) .AND. cols(1, pc)) .OR. (rows(2, pr) .AND. cols(2, pc))
943 END SUBROUTINE kpoint_operator_rs_owners
952 FUNCTION kpoint_operator_rs_compatible(route, source, matrix_row)
RESULT(compatible)
954 TYPE(kpoint_operator_rs_route_type),
INTENT(IN) :: route
956 INTEGER,
INTENT(IN) :: matrix_row
957 LOGICAL :: compatible
960 INTEGER,
POINTER :: col_size(:), row_size(:)
963 IF (
SIZE(source, 2) /= route%nimage)
RETURN
964 IF (matrix_row < 1 .OR. matrix_row >
SIZE(source, 1))
RETURN
965 DO image = 1,
SIZE(source, 2)
966 CALL dbcsr_get_info(source(matrix_row, image)%matrix, row_blk_size=row_size, col_blk_size=col_size)
967 IF (
SIZE(row_size) /=
SIZE(route%row_size) .OR.
SIZE(col_size) /=
SIZE(route%col_size))
RETURN
968 IF (any(row_size /= route%row_size) .OR. any(col_size /= route%col_size))
RETURN
972 END FUNCTION kpoint_operator_rs_compatible
982 SUBROUTINE kpoint_operator_rs_capture(source, matrix_row, route, kpoints, snapshot)
985 INTEGER,
INTENT(IN) :: matrix_row
986 TYPE(kpoint_operator_rs_route_type),
INTENT(INOUT) :: route
988 TYPE(kpoint_operator_rs_type),
INTENT(OUT) :: snapshot
990 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_operator_rs_capture'
992 INTEGER :: col, handle, i, iblock, ncol, nrow, pos, &
995 REAL(kind=
dp),
POINTER :: block(:, :)
998 CALL timeset(routinen, handle)
999 world => kpoints%blacs_env_all%para_env
1000 inter => kpoints%para_env_inter_kp
1001 ALLOCATE (snapshot%blocks(
SIZE(route%entry, 2)), snapshot%nonzero(
SIZE(route%entry, 2)))
1002 snapshot%source => source
1003 snapshot%matrix_row = matrix_row
1004 IF (world%num_pe == 1)
THEN
1006 DO iblock = 1,
SIZE(route%send_entry)
1007 i = route%send_entry(iblock)
1009 route%traversal%row(i), route%traversal%col(i), block, found)
1010 snapshot%nonzero(iblock) = found
1012 cpassert(
SIZE(block, 1) == route%entry(3, iblock))
1013 cpassert(
SIZE(block, 2) == route%entry(4, iblock))
1014 snapshot%blocks(iblock)%data => block
1017 CALL timestop(handle)
1020 ALLOCATE (snapshot%values(route%nvalue))
1021 route%send_values(:) = 0.0_dp
1024 DO iblock = 1,
SIZE(route%send_entry)
1025 i = route%send_entry(iblock)
1027 route%traversal%row(i), route%traversal%col(i), block, found)
1028 IF (.NOT. found) cycle
1029 nrow = route%row_size(route%traversal%row(i))
1030 ncol = route%col_size(route%traversal%col(i))
1031 cpassert(
SIZE(block, 1) == nrow .AND.
SIZE(block, 2) == ncol)
1032 pos = route%send_first(iblock)
1035 route%send_values(pos) = block(row, col)
1041 IF (kpoints%para_env_kp%num_pe == 1)
THEN
1042 CALL world%allgatherv(route%send_values, snapshot%values, route%rcount(2, :), route%rdisp(2, :))
1044 CALL world%alltoall(route%send_values, route%scount(2, :), route%sdisp(2, :), &
1045 snapshot%values, route%rcount(2, :), route%rdisp(2, :))
1046 CALL inter%bcast(snapshot%values, 0)
1048 DO iblock = 1,
SIZE(route%entry, 2)
1049 snapshot%nonzero(iblock) = any(snapshot%values(route%first(iblock):route%first(iblock + 1) - 1) /= 0.0_dp)
1050 snapshot%blocks(iblock)%data(1:route%entry(3, iblock), 1:route%entry(4, iblock)) => &
1051 snapshot%values(route%first(iblock):route%first(iblock + 1) - 1)
1053 snapshot%source => source
1054 snapshot%matrix_row = matrix_row
1055 CALL timestop(handle)
1057 END SUBROUTINE kpoint_operator_rs_capture
1072 SUBROUTINE kpoint_operator_rs_build(ctx, irs, slot, ik, TARGET, fm_re, fm_im)
1075 INTEGER,
INTENT(IN) :: irs, slot, ik
1076 TYPE(
cp_cfm_type),
INTENT(INOUT),
OPTIONAL :: target
1077 TYPE(
cp_fm_type),
INTENT(INOUT),
OPTIONAL :: fm_re, fm_im
1079 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_operator_rs_build'
1081 INTEGER :: first, handle, k, last, ncol, nrow, width
1082 INTEGER(KIND=int_8) :: bytes_per_k
1083 LOGICAL :: layout_ok
1084 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: xkp
1086 CALL timeset(routinen, handle)
1087 IF (
PRESENT(
TARGET))
THEN
1090 cpassert(
PRESENT(fm_re))
1092 IF (
PRESENT(fm_im))
THEN
1098 xkp => ctx%kpoints%xkp
1099 first = ctx%kpoints%kp_range(1)
1100 last = ctx%kpoints%kp_range(2)
1101 cpassert(ik >= first .AND. ik <= last)
1102 nrow = ctx%rs_route%nrow
1103 ncol = ctx%rs_route%ncol
1104 bytes_per_k = 16_int_8*(int(nrow,
int_8)*ncol +
SIZE(ctx%rs_route%entry, 2, kind=
int_8))
1106 int(rs_batch_bytes/(int(ctx%nspin*kpoint_num_slots,
int_8)*max(1_int_8, bytes_per_k))))
1107 IF (width < 2 .OR. nrow == 0 .OR. ncol == 0)
THEN
1109 CALL kpoint_operator_rs_build_batch(ctx%rs(irs, slot), ctx%rs_route, xkp(:, ik:ik), &
1110 TARGET=
TARGET, fm_re=fm_re, fm_im=fm_im)
1112 associate(work => ctx%rs_batch(irs, slot))
1113 IF (.NOT.
ALLOCATED(work%data))
ALLOCATE (work%data(width, nrow, ncol))
1114 IF (work%first == 0 .OR. ik < work%first .OR. ik >= work%first + width)
THEN
1115 work%first = first + ((ik - first)/width)*width
1116 CALL kpoint_operator_rs_build_batch(ctx%rs(irs, slot), ctx%rs_route, &
1117 xkp(:, work%first:min(work%first + width - 1, last)), batch=work%data)
1119 k = ik - work%first + 1
1121 IF (
PRESENT(
TARGET))
THEN
1122 target%local_data(1:nrow, 1:ncol) = work%data(k, :, :)
1123 target%local_data(nrow + 1:, :) =
z_zero
1124 target%local_data(:, ncol + 1:) =
z_zero
1126 cpassert(
PRESENT(fm_re))
1127 fm_re%local_data(1:nrow, 1:ncol) = real(work%data(k, :, :),
dp)
1128 IF (
PRESENT(fm_im)) fm_im%local_data(1:nrow, 1:ncol) = aimag(work%data(k, :, :))
1129 fm_re%local_data(nrow + 1:, :) = 0.0_dp
1130 fm_re%local_data(:, ncol + 1:) = 0.0_dp
1131 IF (
PRESENT(fm_im)) fm_im%local_data(nrow + 1:, :) = 0.0_dp
1132 IF (
PRESENT(fm_im)) fm_im%local_data(:, ncol + 1:) = 0.0_dp
1136 IF (ctx%kpoints%use_real_wfn)
THEN
1137 IF (
PRESENT(
TARGET))
target%local_data(:, :) = cmplx(real(
target%local_data,
dp), 0.0_dp,
dp)
1138 IF (
PRESENT(fm_im)) fm_im%local_data(:, :) = 0.0_dp
1140 CALL timestop(handle)
1142 END SUBROUTINE kpoint_operator_rs_build
1160 SUBROUTINE kpoint_operator_rs_build_batch(snapshot, route, xkp, batch, TARGET, fm_re, fm_im, phase_work, parallel_columns)
1162 TYPE(kpoint_operator_rs_type),
INTENT(IN) :: snapshot
1163 TYPE(kpoint_operator_rs_route_type),
INTENT(IN) :: route
1164 REAL(kind=dp),
INTENT(IN) :: xkp(:, :)
1165 COMPLEX(KIND=dp),
CONTIGUOUS,
INTENT(INOUT), &
1166 OPTIONAL,
TARGET :: batch(:, :, :)
1167 TYPE(cp_cfm_type),
INTENT(INOUT),
OPTIONAL :: target
1168 TYPE(cp_fm_type),
INTENT(INOUT),
OPTIONAL :: fm_re, fm_im
1169 COMPLEX(KIND=dp),
INTENT(INOUT),
OPTIONAL,
TARGET :: phase_work(:, :)
1170 LOGICAL,
INTENT(IN),
OPTIONAL :: parallel_columns
1172 COMPLEX(KIND=dp),
ALLOCATABLE,
TARGET :: owned_phase(:, :)
1173 COMPLEX(KIND=dp),
CONTIGUOUS,
POINTER :: output(:, :, :)
1174 COMPLEX(KIND=dp),
POINTER :: phase(:, :)
1175 INTEGER :: first, iblock, k, last, nbatch, ncol, &
1178 REAL(kind=dp) :: arg
1180 nbatch =
SIZE(xkp, 2)
1181 cpassert(nbatch >= 1 .AND. nbatch <=
rs_batch_size .AND.
SIZE(xkp, 1) == 3)
1182 cpassert(count([
PRESENT(batch),
PRESENT(
TARGET),
PRESENT(fm_re)]) == 1)
1184 IF (
PRESENT(batch))
THEN
1185 cpassert(
SIZE(batch, 1) >= nbatch)
1188 cpassert(nbatch == 1)
1189 IF (
PRESENT(
TARGET))
THEN
1190 output(1:1, 1:
SIZE(
target%local_data, 1), 1:
SIZE(
target%local_data, 2)) =>
target%local_data
1192 cpassert(
PRESENT(fm_re))
1193 fm_re%local_data(:, :) = 0.0_dp
1194 IF (
PRESENT(fm_im)) fm_im%local_data(:, :) = 0.0_dp
1197 IF (
ASSOCIATED(output)) output(:, :, :) = z_zero
1198 IF (route%nrow == 0 .OR. route%ncol == 0)
RETURN
1199 IF (
PRESENT(phase_work))
THEN
1200 cpassert(
SIZE(phase_work, 1) >= nbatch)
1201 cpassert(
SIZE(phase_work, 2) ==
SIZE(route%entry, 2))
1202 phase => phase_work(1:nbatch, :)
1204 ALLOCATE (owned_phase(nbatch,
SIZE(route%entry, 2)))
1205 phase => owned_phase
1207 DO iblock = 1,
SIZE(route%entry, 2)
1208 IF (.NOT. snapshot%nonzero(iblock)) cycle
1210 arg = twopi*dot_product(real(route%entry(5:7, iblock), dp), xkp(:, k))
1211 phase(k, iblock) = cmplx(cos(arg), sin(arg), kind=dp)
1216 IF (
PRESENT(parallel_columns)) parallel = parallel_columns
1219 IF (nthreads > 1)
THEN
1226 first = ncol*thread/nthreads + 1
1227 last = ncol*(thread + 1)/nthreads
1228 CALL kpoint_operator_rs_accumulate(snapshot, route, phase, first, last, output, fm_re, fm_im)
1232 CALL kpoint_operator_rs_accumulate(snapshot, route, phase, 1, ncol, output, fm_re, fm_im)
1235 END SUBROUTINE kpoint_operator_rs_build_batch
1248 SUBROUTINE kpoint_operator_rs_accumulate(snapshot, route, phase, first, last, output, fm_re, fm_im)
1249 TYPE(kpoint_operator_rs_type),
INTENT(IN) :: snapshot
1250 TYPE(kpoint_operator_rs_route_type),
INTENT(IN) :: route
1251 COMPLEX(KIND=dp),
INTENT(IN) :: phase(:, :)
1252 INTEGER,
INTENT(IN) :: first, last
1253 COMPLEX(KIND=dp),
CONTIGUOUS,
INTENT(IN),
POINTER :: output(:, :, :)
1254 TYPE(cp_fm_type),
INTENT(INOUT),
OPTIONAL :: fm_re, fm_im
1257 COMPLEX(KIND=dp),
CONTIGUOUS,
POINTER :: full(:, :)
1258 INTEGER :: c0, col, hi, iblock, j, k, lo, nc, nr, &
1259 part, r0, row, source_col, source_row
1260 REAL(kind=dp) :: element
1262 IF (first > last)
RETURN
1264 IF (
ASSOCIATED(output))
THEN
1265 IF (
SIZE(output, 1) == 1 .AND. &
1266 route%nrow == route%matrix_struct%nrow_global .AND. &
1267 route%ncol == route%matrix_struct%ncol_global)
THEN
1268 full(1:
SIZE(output, 2), 1:
SIZE(output, 3)) => output
1271 DO iblock = 1,
SIZE(route%entry, 2)
1272 IF (.NOT. snapshot%nonzero(iblock)) cycle
1273 r0 = route%entry(1, iblock)
1274 c0 = route%entry(2, iblock)
1275 nr = route%entry(3, iblock)
1276 nc = route%entry(4, iblock)
1277 IF (
ASSOCIATED(full))
THEN
1280 hi = min(last, c0 + nc - 1)
1282 full(r0:r0 + nr - 1, lo:hi) = full(r0:r0 + nr - 1, lo:hi) + &
1283 phase(1, iblock)*snapshot%blocks(iblock)%data(:, lo - c0 + 1:hi - c0 + 1)
1285 IF (route%entry(8, iblock) /= 0)
THEN
1287 hi = min(last, r0 + nr - 1)
1289 full(c0:c0 + nc - 1, lo:hi) = full(c0:c0 + nc - 1, lo:hi) + &
1290 conjg(phase(1, iblock))* &
1291 transpose(snapshot%blocks(iblock)%data(lo - r0 + 1:hi - r0 + 1, :))
1297 IF (part == 2 .AND. route%entry(8, iblock) == 0) cycle
1298 lo = max(first, route%tile(3, part, iblock))
1299 hi = min(last, route%tile(4, part, iblock))
1301 factors(1:
SIZE(phase, 1)) = phase(:, iblock)
1302 IF (part == 2) factors(1:
SIZE(phase, 1)) = conjg(factors(1:
SIZE(phase, 1)))
1304 col = route%matrix_struct%col_indices(j)
1305 DO row = route%tile(1, part, iblock), route%tile(2, part, iblock)
1306 source_row = route%matrix_struct%row_indices(row)
1309 source_col = source_row
1312 element = snapshot%blocks(iblock)%data(source_row - r0 + 1, source_col - c0 + 1)
1313 IF (
ASSOCIATED(output))
THEN
1315 DO k = 1,
SIZE(phase, 1)
1316 output(k, row, j) = output(k, row, j) + factors(k)*element
1319 value = factors(1)*element
1320 fm_re%local_data(row, j) = fm_re%local_data(row, j) + real(
value, dp)
1321 IF (
PRESENT(fm_im)) fm_im%local_data(row, j) = fm_im%local_data(row, j) + aimag(
value)
1327 END SUBROUTINE kpoint_operator_rs_accumulate
1336 INTEGER,
INTENT(OUT) :: nentry
1338 cpassert(
ALLOCATED(ctx%rs_route))
1339 nentry =
SIZE(ctx%rs_route%entry, 2)
1357 INTEGER,
INTENT(IN) :: ispin, slot, ik
1358 TYPE(cp_cfm_type),
INTENT(INOUT),
OPTIONAL :: target
1359 COMPLEX(KIND=dp),
INTENT(INOUT),
TARGET :: phase(:, :)
1360 TYPE(cp_fm_type),
INTENT(INOUT),
OPTIONAL :: fm_re, fm_im
1362 REAL(kind=dp),
POINTER :: xkp(:, :)
1365 cpassert(ik >= ctx%kpoints%kp_range(1))
1366 cpassert(ik <= ctx%kpoints%kp_range(2))
1368 IF (.NOT.
ASSOCIATED(ctx%rs(row, slot)%source)) row = 1
1369 cpassert(
ASSOCIATED(ctx%rs(row, slot)%source))
1370 xkp => ctx%kpoints%xkp
1371 CALL kpoint_operator_rs_build_batch(ctx%rs(row, slot), ctx%rs_route, xkp(:, ik:ik), &
1372 TARGET=
TARGET, fm_re=fm_re, fm_im=fm_im, &
1373 phase_work=phase, parallel_columns=.false.)
1374 IF (ctx%kpoints%use_real_wfn)
THEN
1375 IF (
PRESENT(
TARGET))
target%local_data(:, :) = cmplx(real(
target%local_data, dp), 0.0_dp, dp)
1376 IF (
PRESENT(fm_im)) fm_im%local_data(:, :) = 0.0_dp
1407 INTEGER,
INTENT(IN) :: ikp, ispin
1408 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: rsmat
1409 INTEGER,
INTENT(IN) :: slot
1410 TYPE(rskp_grid_type),
INTENT(IN),
OPTIONAL :: grid
1411 LOGICAL,
INTENT(IN) :: use_grid
1412 INTEGER,
INTENT(IN) :: matrix_row
1414 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_operator_start'
1416 INTEGER :: handle, igroup, ik, irs, kpmax
1417 INTEGER,
DIMENSION(:, :),
POINTER :: kp_dist
1419 CALL timeset(routinen, handle)
1421 cpassert(
ASSOCIATED(ctx%kpoints))
1422 cpassert(
ASSOCIATED(ctx%kpoints%sab_nl))
1427 cpassert(ispin >= 1 .AND. ispin <= ctx%nspin)
1428 cpassert(matrix_row >= 1 .AND. matrix_row <=
SIZE(rsmat, 1))
1429 kpmax = maxval(ctx%kpoints%kp_dist(2, :) - ctx%kpoints%kp_dist(1, :) + 1)
1430 cpassert(ikp >= 1 .AND. ikp <= kpmax)
1434 IF (.NOT.
ALLOCATED(ctx%status))
THEN
1435 ALLOCATE (ctx%status(kpmax, ctx%nspin, kpoint_num_slots))
1436 ctx%status(:, :, :) = kpoint_slot_idle
1438 cpassert(ctx%status(ikp, ispin, slot) == kpoint_slot_idle)
1439 IF (.NOT. use_grid)
THEN
1440 irs = kpoint_operator_rs_match(ctx, rsmat, ispin, slot, matrix_row)
1442 cpassert(
ASSOCIATED(ctx%work_sab, ctx%kpoints%sab_nl))
1443 ctx%status(ikp, ispin, slot) = kpoint_slot_local
1444 CALL timestop(handle)
1448 CALL kpoint_operator_ensure(ctx, rsmat)
1453 IF (.NOT.
ALLOCATED(ctx%info))
THEN
1456 ALLOCATE (ctx%info(
SIZE(ctx%kpoints%xkp, 2), ctx%nspin, &
1457 kpoint_num_slots, merge(1, 2, ctx%kpoints%use_real_wfn)))
1460 kp_dist => ctx%kpoints%kp_dist
1461 DO igroup = 1, ctx%kpoints%nkp_groups
1462 ik = kp_dist(1, igroup) + ikp - 1
1463 IF (ik > kp_dist(2, igroup)) cycle
1465 CALL kpoint_operator_launch(ctx%kpoints, ik, ctx%fmwork, ctx%fmlocal, &
1466 ctx%info(ik, ispin, slot, :))
1468 ctx%status(ikp, ispin, slot) = kpoint_slot_in_flight
1470 CALL timestop(handle)
1493 INTEGER,
INTENT(IN) :: ikp, ispin, slot
1494 TYPE(cp_cfm_type),
OPTIONAL :: cmat
1495 TYPE(cp_fm_type),
OPTIONAL :: fm_re
1497 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_operator_finish'
1499 INTEGER :: handle, igroup, ik, irs, kpmax
1500 INTEGER,
DIMENSION(2) :: kp_range
1501 INTEGER,
DIMENSION(:, :),
POINTER :: kp_dist
1502 LOGICAL :: layout_ok
1504 CALL timeset(routinen, handle)
1506 cpassert(
ASSOCIATED(ctx%kpoints))
1507 cpassert(
ALLOCATED(ctx%status))
1508 kpmax =
SIZE(ctx%status, 1)
1509 cpassert(ikp >= 1 .AND. ikp <= kpmax)
1510 cpassert(ispin >= 1 .AND. ispin <= ctx%nspin)
1512 cpassert(ctx%status(ikp, ispin, slot) /= kpoint_slot_idle)
1516 cpassert(
PRESENT(cmat) .NEQV.
PRESENT(fm_re))
1517 cpassert(.NOT.
PRESENT(fm_re) .OR. ctx%kpoints%use_real_wfn)
1518 kp_range = ctx%kpoints%kp_range
1519 IF (ctx%status(ikp, ispin, slot) == kpoint_slot_local)
THEN
1521 IF (.NOT.
ASSOCIATED(ctx%rs(irs, slot)%source)) irs = 1
1522 ik = kp_range(1) + ikp - 1
1523 IF (ik <= kp_range(2))
THEN
1524 CALL kpoint_operator_rs_build(ctx, irs, slot, ik, cmat, fm_re=fm_re)
1526 ctx%status(ikp, ispin, slot) = kpoint_slot_idle
1527 CALL timestop(handle)
1530 cpassert(ctx%status(ikp, ispin, slot) == kpoint_slot_in_flight)
1531 IF (
PRESENT(cmat))
THEN
1532 layout_ok = cp_fm_struct_equivalent(cmat%matrix_struct, ctx%fmlocal%matrix_struct)
1534 layout_ok = cp_fm_struct_equivalent(fm_re%matrix_struct, ctx%fmlocal%matrix_struct)
1537 kp_dist => ctx%kpoints%kp_dist
1538 DO igroup = 1, ctx%kpoints%nkp_groups
1539 ik = kp_dist(1, igroup) + ikp - 1
1540 IF (ik > kp_dist(2, igroup)) cycle
1541 CALL kpoint_operator_receive(ctx%kpoints, ctx%fmlocal, ik, ctx%info(ik, ispin, slot, :), &
1542 cmat=cmat, fm_re=fm_re)
1544 ctx%status(ikp, ispin, slot) = kpoint_slot_idle
1546 CALL timestop(handle)
1588 mat_re, mat_im, matrix_row, slot, grid)
1591 INTEGER,
INTENT(IN) :: ikp, ispin
1592 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: rsmat
1593 TYPE(cp_cfm_type),
OPTIONAL :: cmat
1594 TYPE(cp_fm_type),
OPTIONAL :: fm_re, fm_im
1595 TYPE(dbcsr_type),
INTENT(INOUT),
OPTIONAL :: mat_re, mat_im
1596 INTEGER,
INTENT(IN) :: matrix_row
1597 INTEGER,
INTENT(IN),
OPTIONAL :: slot
1598 TYPE(rskp_grid_type),
INTENT(IN),
OPTIONAL :: grid
1600 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_operator_get'
1602 INTEGER :: handle, igroup, ik, irs, kpmax
1603 INTEGER,
DIMENSION(:, :),
POINTER :: kp_dist
1604 LOGICAL :: layout_ok
1605 TYPE(copy_info_type) :: info(2)
1607 CALL timeset(routinen, handle)
1609 cpassert(count([
PRESENT(cmat),
PRESENT(fm_re),
PRESENT(mat_re)]) == 1)
1610 IF (
PRESENT(fm_re) .OR.
PRESENT(fm_im))
THEN
1611 cpassert(
PRESENT(fm_re))
1612 cpassert(
PRESENT(fm_im) .OR. ctx%kpoints%use_real_wfn)
1614 IF (
PRESENT(mat_re) .OR.
PRESENT(mat_im))
THEN
1615 cpassert(
PRESENT(mat_re) .AND.
PRESENT(mat_im))
1618 cpassert(
ASSOCIATED(ctx%kpoints))
1619 cpassert(
ASSOCIATED(ctx%kpoints%sab_nl))
1620 cpassert(ispin >= 1 .AND. ispin <= ctx%nspin)
1621 cpassert(matrix_row >= 1 .AND. matrix_row <=
SIZE(rsmat, 1))
1622 kpmax = maxval(ctx%kpoints%kp_dist(2, :) - ctx%kpoints%kp_dist(1, :) + 1)
1623 cpassert(ikp >= 1 .AND. ikp <= kpmax)
1624 IF (
PRESENT(slot) .AND. .NOT.
PRESENT(grid))
THEN
1626 irs = kpoint_operator_rs_match(ctx, rsmat, ispin, slot, matrix_row)
1628 ik = ctx%kpoints%kp_range(1) + ikp - 1
1629 IF (ik <= ctx%kpoints%kp_range(2))
THEN
1630 IF (
PRESENT(mat_re))
THEN
1631 IF (.NOT.
ALLOCATED(ctx%fmimag))
THEN
1632 ALLOCATE (ctx%fmimag)
1633 CALL cp_fm_create(ctx%fmimag, ctx%fmlocal%matrix_struct)
1635 CALL kpoint_operator_rs_build(ctx, irs, slot, ik, fm_re=ctx%fmlocal, fm_im=ctx%fmimag)
1636 CALL copy_fm_to_dbcsr_bc(ctx%fmlocal, mat_re)
1637 CALL copy_fm_to_dbcsr_bc(ctx%fmimag, mat_im)
1639 CALL kpoint_operator_rs_build(ctx, irs, slot, ik,
TARGET=cmat, fm_re=fm_re, fm_im=fm_im)
1642 CALL timestop(handle)
1646 CALL kpoint_operator_ensure(ctx, rsmat)
1647 IF (
PRESENT(cmat))
THEN
1650 layout_ok = cp_fm_struct_equivalent(cmat%matrix_struct, ctx%fmlocal%matrix_struct)
1654 kp_dist => ctx%kpoints%kp_dist
1655 DO igroup = 1, ctx%kpoints%nkp_groups
1656 ik = kp_dist(1, igroup) + ikp - 1
1657 IF (ik > kp_dist(2, igroup)) cycle
1661 CALL kpoint_operator_launch(ctx%kpoints, ik, ctx%fmwork, ctx%fmlocal, info, &
1662 target_re=fm_re, target_im=fm_im)
1663 CALL kpoint_operator_receive(ctx%kpoints, ctx%fmlocal, ik, info, cmat=cmat, fm_re=fm_re, fm_im=fm_im, &
1664 mat_re=mat_re, mat_im=mat_im)
1667 CALL timestop(handle)
1685 matrix_re, matrix_im)
1686 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_rs
1687 TYPE(kpoint_type),
POINTER :: kpoints
1688 TYPE(kpoint_env_type),
POINTER :: kp
1689 INTEGER,
INTENT(IN) :: ispin
1690 TYPE(cp_fm_type),
INTENT(IN) :: cache_re, cache_im
1691 TYPE(dbcsr_type),
POINTER :: matrix_re, matrix_im
1693 INTEGER :: nkp_groups
1694 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1695 REAL(dp),
DIMENSION(:, :),
POINTER :: xkp
1696 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
1699 NULLIFY (cell_to_index, sab_nl, xkp)
1700 cpassert(
ASSOCIATED(matrix_rs))
1701 cpassert(ispin >= 1 .AND. ispin <=
SIZE(matrix_rs, 1))
1702 cpassert(
ASSOCIATED(kpoints))
1703 cpassert(
ASSOCIATED(kp))
1705 CALL dbcsr_init_p(matrix_re)
1706 CALL dbcsr_init_p(matrix_im)
1707 nkp_groups = kpoints%nkp_groups
1709 cell_to_index => kpoints%cell_to_index
1710 sab_nl => kpoints%sab_nl
1711 IF (nkp_groups == 1)
THEN
1712 cpassert(
ASSOCIATED(xkp))
1713 cpassert(
ASSOCIATED(cell_to_index))
1714 cpassert(
ASSOCIATED(sab_nl))
1715 CALL dbcsr_create(matrix_re, template=matrix_rs(1, 1)%matrix, &
1716 matrix_type=dbcsr_type_symmetric)
1717 CALL dbcsr_create(matrix_im, template=matrix_rs(1, 1)%matrix, &
1718 matrix_type=dbcsr_type_antisymmetric)
1719 CALL cp_dbcsr_alloc_block_from_nbl(matrix_re, sab_nl)
1720 CALL cp_dbcsr_alloc_block_from_nbl(matrix_im, sab_nl)
1721 CALL dbcsr_set(matrix_re, 0.0_dp)
1722 CALL dbcsr_set(matrix_im, 0.0_dp)
1723 CALL rskp_transform(rmatrix=matrix_re, cmatrix=matrix_im, rsmat=matrix_rs, ispin=ispin, &
1724 xkp=xkp(1:3, kp%nkpoint), cell_to_index=cell_to_index, sab_nl=sab_nl)
1726 CALL copy_fm_to_dbcsr_bc(cache_re, matrix_re)
1727 CALL copy_fm_to_dbcsr_bc(cache_im, matrix_im)
1743 TYPE(kpoint_type),
POINTER :: kpoints
1744 TYPE(cp_fm_type),
INTENT(IN) :: ao_ao_fm
1745 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: rsmat
1746 INTEGER,
INTENT(IN) :: slot
1747 TYPE(cp_cfm_type),
DIMENSION(:),
INTENT(INOUT) :: cmat
1748 INTEGER,
INTENT(IN) :: matrix_row
1750 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_operator_get_all'
1752 COMPLEX(KIND=dp),
ALLOCATABLE :: phase(:, :)
1753 INTEGER :: handle, ik, ikp, kplocal, kpmax, nentry, &
1755 LOGICAL :: layout_ok, prepared
1758 CALL timeset(routinen, handle)
1759 kplocal =
SIZE(kpoints%kp_env)
1760 cpassert(
SIZE(cmat) == kplocal .AND. kplocal > 0)
1763 IF (prepared .AND. kpoints%para_env_kp%num_pe == 1)
THEN
1765 layout_ok = cp_fm_struct_equivalent(cmat(ikp)%matrix_struct, ctx%rs_route%matrix_struct)
1773 ALLOCATE (phase(1, nentry))
1776 ik = kpoints%kp_range(1) + ikp - 1
1783 kpmax = maxval(kpoints%kp_dist(2, :) - kpoints%kp_dist(1, :) + 1)
1786 matrix_row=matrix_row, slot=slot)
1790 CALL timestop(handle)
1804 TYPE(kpoint_type),
POINTER :: kpoints
1805 TYPE(cp_fm_type),
INTENT(IN) :: ao_ao_fm
1806 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, matrix_s
1807 TYPE(dbcsr_p_type),
DIMENSION(:, :),
OPTIONAL, &
1810 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_operator_store'
1811 INTEGER :: handle, ic, ikp, ispin, ncomp, nspin
1812 LOGICAL :: kinetic, reuse
1813 TYPE(cp_fm_struct_type),
POINTER :: matrix_struct
1814 TYPE(kpoint_env_type),
POINTER :: kp
1817 CALL timeset(routinen, handle)
1818 nspin =
SIZE(matrix_ks, 1)
1819 ncomp = merge(1, 2, kpoints%use_real_wfn)
1820 kinetic =
PRESENT(matrix_t)
1821 IF (kinetic) kinetic =
ASSOCIATED(matrix_t)
1823 CALL kpoint_operator_ensure_group(ctx)
1824 matrix_struct => ctx%fmlocal%matrix_struct
1825 DO ikp = 1,
SIZE(kpoints%kp_env)
1826 kp => kpoints%kp_env(ikp)%kpoint_env
1828 IF (
ASSOCIATED(kp%ot_hmat))
THEN
1829 IF (all(shape(kp%ot_hmat) == [ncomp, nspin]))
THEN
1830 reuse = cp_fm_struct_equivalent(kp%ot_hmat(1, 1)%matrix_struct, matrix_struct)
1833 IF (.NOT. reuse)
THEN
1834 CALL cp_fm_release(kp%ot_hmat)
1835 ALLOCATE (kp%ot_hmat(ncomp, nspin))
1838 CALL cp_fm_create(kp%ot_hmat(ic, ispin), matrix_struct)
1842 CALL allocate_spin_free(kp%ot_smat)
1843 IF (kinetic)
CALL allocate_spin_free(kp%ot_tmat)
1849 CALL timestop(handle)
1857 SUBROUTINE allocate_spin_free(cache)
1858 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: cache
1864 IF (
ASSOCIATED(cache))
THEN
1865 IF (
SIZE(cache) == ncomp) reuse = cp_fm_struct_equivalent(cache(1)%matrix_struct, matrix_struct)
1868 CALL cp_fm_release(cache)
1869 ALLOCATE (cache(ncomp))
1871 CALL cp_fm_create(cache(ic), matrix_struct)
1873 END SUBROUTINE allocate_spin_free
1881 SUBROUTINE store_slot(rsmat, slot, nrows)
1882 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: rsmat
1883 INTEGER,
INTENT(IN) :: slot, nrows
1885 COMPLEX(KIND=dp),
ALLOCATABLE :: phase(:, :)
1886 INTEGER :: ikp, ispin, kplocal, kpmax, nentry, &
1888 LOGICAL :: prepared, ready, use_grid
1889 TYPE(rskp_grid_type),
ALLOCATABLE :: grids(:)
1893 use_grid = kpoints%lattice_fft == lattice_fft_on .OR. &
1894 (kpoints%lattice_fft == lattice_fft_auto .AND.
SIZE(kpoints%xkp, 2) >= 27)
1897 IF (.NOT. use_grid .OR. nworkers > 1 .OR. kpoints%blacs_env_all%para_env%num_pe > 1)
THEN
1904 kplocal =
SIZE(kpoints%kp_env)
1905 IF (prepared .AND. kpoints%para_env_kp%num_pe == 1)
THEN
1907 nworkers = min(nworkers, kplocal*nrows)
1910 ALLOCATE (phase(1, nentry))
1914 CALL store_round(rsmat, slot, ikp, ispin, phase)
1921 kpmax = maxval(kpoints%kp_dist(2, :) - kpoints%kp_dist(1, :) + 1)
1923 use_grid = use_grid .AND. .NOT. prepared
1925 CALL kpoint_operator_ensure(ctx, rsmat)
1926 ALLOCATE (grids(nrows))
1928 CALL rskp_transform_grid_prepare(grids(ispin), ctx%rmatrix, rsmat, ispin, &
1929 kpoints%xkp, kpoints%nkp_grid, kpoints%cell_to_index, &
1930 kpoints%sab_nl, ready, &
1931 max_storage_bytes=512_int_8*1024_int_8**2/int(nspin + 1, int_8))
1932 use_grid = use_grid .AND. ready
1938 CALL store_round(rsmat, slot, ikp, ispin, grid=grids(ispin))
1940 CALL store_round(rsmat, slot, ikp, ispin)
1944 IF (
ALLOCATED(grids))
THEN
1946 CALL rskp_transform_grid_release(grids(ispin))
1951 END SUBROUTINE store_slot
1962 SUBROUTINE store_round(rsmat, slot, ikp, ispin, phase, grid)
1963 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: rsmat
1964 INTEGER,
INTENT(IN) :: slot, ikp, ispin
1965 COMPLEX(KIND=dp),
INTENT(INOUT),
OPTIONAL :: phase(:, :)
1966 TYPE(rskp_grid_type),
INTENT(IN),
OPTIONAL :: grid
1969 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: cache
1970 TYPE(cp_fm_type),
POINTER :: imag
1971 TYPE(kpoint_env_type),
POINTER :: kp
1974 kp => kpoints%kp_env(min(ikp,
SIZE(kpoints%kp_env)))%kpoint_env
1977 cache => kp%ot_hmat(:, ispin)
1984 IF (
SIZE(cache) == 2) imag => cache(2)
1985 IF (
PRESENT(phase))
THEN
1986 ik = kpoints%kp_range(1) + ikp - 1
1990 matrix_row=ispin, slot=slot, grid=grid)
1992 END SUBROUTINE store_round
Basic linear algebra operations for complex full matrices.
subroutine, public cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b). where b is a real matrix (adapted from cp_cf...
Represents a complex full matrix distributed on many processors.
subroutine, public dbcsr_get_readonly_block_p(matrix, row, col, block, found, row_size, col_size)
Like dbcsr_get_block_p() but with matrix being INTENT(IN). When invoking this routine,...
subroutine, public dbcsr_deallocate_matrix(matrix)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_set(matrix, alpha)
...
Routines that link DBCSR and CP2K concepts together.
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_fm_to_dbcsr_bc(fm, bc_mat)
Copy a BLACS matrix to a dbcsr matrix with a special block-cyclic distribution, which requires no com...
pool for for elements that are retained and released
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
subroutine, public fm_pool_give_back_fm(pool, element)
returns the element to the pool
represent the structure of a full matrix
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_retain(fmstruct)
retains a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_start_copy_general(source, destination, para_env, info)
Initiates the copy operation: get distribution data, post MPI isend and irecvs.
subroutine, public cp_fm_cleanup_copy_general(info)
Completes the copy operation: wait for comms clean up MPI state.
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_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
subroutine, public cp_fm_finish_copy_general(destination, info)
Completes the copy operation: wait for comms, unpack, clean up MPI state.
Defines the basic variable types.
integer, parameter, public int_8
integer, parameter, public dp
Routines needed for kpoint calculation.
subroutine, public rskp_transform(rmatrix, cmatrix, rsmat, ispin, xkp, cell_to_index, sab_nl, is_complex, rs_sign)
Transformation of real space matrices to a kpoint.
subroutine, public rskp_transform_grid_extract(grid, ikp, rmatrix, cmatrix)
Extract one reciprocal-grid matrix from a prepared local DBCSR block cache.
subroutine, public rskp_transform_grid_release(grid)
Release a reciprocal-grid transformation cache.
subroutine, public kp_transform_plan_create(plan, sab_nl, cell_to_index, nimg, block_template, group_entries)
Build the immutable neighbor-list traversal shared by R-to-K and K-to-R transforms.
subroutine, public rskp_transform_grid_prepare(grid, rmatrix, rsmat, ispin, xkp, nkp_grid, cell_to_index, sab_nl, used_fft, is_complex, rs_sign, max_storage_bytes)
Prepare a batched real-cell to complete reciprocal-grid transform for local DBCSR blocks....
Types and basic routines needed for a kpoint calculation.
Definition of mathematical constants and functions.
complex(kind=dp), parameter, public z_one
complex(kind=dp), parameter, public gaussi
real(kind=dp), parameter, public twopi
complex(kind=dp), parameter, public z_zero
Interface to the message passing library MPI.
Assembly of real and complex k-point operators from real-space DBCSR matrices. Outputs use the distri...
subroutine, public kpoint_operator_context_create(ctx, kpoints, ao_ao_fm, nspin)
Create an empty operator context. The context owns no buffers after this call; they are created on fi...
subroutine, public kpoint_operator_prepare(ctx, rsmat, slot, prepared, matrix_row, target)
Prepare an immutable operator slot for group-local assembly. Refresh after the source values change,...
subroutine, public kpoint_operator_finish(ctx, ikp, ispin, slot, cmat, fm_re)
Finish one started instance: build its prepared group-local operator, or merge received parts into cm...
subroutine, public kpoint_operator_work_size(ctx, nentry)
Size per-worker phase scratch for direct assembly into a solver matrix.
subroutine, public kpoint_operator_get(ctx, ikp, ispin, rsmat, cmat, fm_re, fm_im, mat_re, mat_im, matrix_row, slot, grid)
Assemble one operator for one global round and spin in one call: run the single-group transfer of the...
integer, parameter, public kpoint_spin_free
Row of rsmat that carries the spin-free operator image. S and T carry no spin copies: drivers launch ...
subroutine, public kpoint_operator_evaluate(ctx, ispin, slot, ik, phase, target, fm_re, fm_im)
Evaluate an immutable prepared operator with exclusively caller-owned output and scratch....
integer, parameter, public kpoint_slot_s
integer, parameter rs_batch_size
subroutine kpoint_operator_densify(ctx, rsmat, matrix_row, ik, grid, use_grid)
Transform and densify one operator using the context-owned workspaces and routing plans.
subroutine, public kpoint_operator_store(kpoints, ao_ao_fm, matrix_ks, matrix_s, matrix_t)
Prepare owned split-FM operators for OT without entering an SCF eigensolver. H (or the initialization...
integer, parameter, public kpoint_slot_t
integer, parameter, public kpoint_slot_ks
subroutine, public kpoint_operator_get_local(matrix_rs, kpoints, kp, ispin, cache_re, cache_im, matrix_re, matrix_im)
Read an operator inside a group-local OT loop, after collective cache preparation....
subroutine, public kpoint_operator_get_all(kpoints, ao_ao_fm, rsmat, slot, cmat, matrix_row)
Fill caller-owned matrices for one operator row, including independent history snapshots....
subroutine, public kpoint_operator_start(ctx, ikp, ispin, rsmat, slot, grid, use_grid, matrix_row)
Start one operator for one global round and spin. The service walks all groups with a k point in this...
subroutine, public kpoint_operator_context_release(ctx)
Release the context together with every buffer it created. Call once, at the exit of the routine that...
wrapper for the pools of matrixes
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
Define the neighbor list data types and the corresponding functionality.
Represent a complex full matrix.
represent a pool of elements with the same structure
keeps the information about the structure of a full matrix
Stores the state of a copy between cp_fm_start_copy_general and cp_fm_finish_copy_general.
Keeps information about a specific k-point.
Contains information about kpoints.
stores all the informations relevant to an mpi environment
Caller-owned assembly context for one driver call. Created and released in the same driver routine an...