Skip to content
Permalink
Browse files

QS: Workarounds to allow for explicitly shaped force_thread array as …

…needed for OMP REDUCTION
  • Loading branch information
oschuett committed Jan 11, 2020
1 parent fe33886 commit 597fd47498f20c64e174034c8826a69f9b7252d7
Showing with 194 additions and 72 deletions.
  1. +8 −9 src/core_ae.F
  2. +19 −21 src/core_ppl.F
  3. +10 −11 src/core_ppnl.F
  4. +53 −9 src/qs_integrate_potential_product.F
  5. +51 −8 src/qs_kinetic.F
  6. +53 −14 src/qs_overlap.F
@@ -106,18 +106,19 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
REAL(KIND=dp) :: alpha_c, core_charge, core_radius, dab, &
dac, dbc, f0, rab2, rac2, rbc2, zeta_c
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ff
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: force_thread, habd, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: habd, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab, pab, verf, vnuc
REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, rab, rac, rbc
REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc, pv_thread
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, rpgfa, rpgfb, sphi_a, &
sphi_b, zeta, zetb
TYPE(all_potential_type), POINTER :: all_potential
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: ap_iterator
TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(all_potential_type), POINTER :: all_potential
REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, rpgfa, rpgfb, sphi_a, &
sphi_b, zeta, zetb
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
TYPE(sgp_potential_type), POINTER :: sgp_potential

!$ INTEGER(kind=omp_lock_kind), &
@@ -146,7 +147,6 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
END DO
END IF
END IF
ALLOCATE (force_thread(3, natom))
force_thread = 0.0_dp

IF (use_virial) THEN
@@ -470,7 +470,6 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
!$OMP END DO
DEALLOCATE (atom_of_kind, kind_of)
END IF
DEALLOCATE (force_thread)

IF (calculate_forces .AND. use_virial) THEN
virial%pv_virial = virial%pv_virial + pv_thread
@@ -115,23 +115,24 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
INTEGER, DIMENSION(nexp_max) :: nct_ppl
LOGICAL :: dokp, ecp_local, found, lpotextended
REAL(KIND=dp) :: alpha, dab, dac, dbc, f0, ppl_radius
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: force_thread, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ppl_fwork, ppl_work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: hab, pab
REAL(KIND=dp), DIMENSION(1:10) :: aloc, bloc
REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, rab, rac, rbc
REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc, pv_thread
REAL(KIND=dp), DIMENSION(4, nexp_max) :: cval_ppl
REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, h_block, p_block, rpgfa, &
rpgfb, sphi_a, sphi_b, zeta, zetb
REAL(KIND=dp), DIMENSION(nexp_max) :: alpha_ppl
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: ap_iterator
TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gth_potential_type), POINTER :: gth_potential
REAL(KIND=dp), DIMENSION(nexp_max) :: alpha_ppl
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, h_block, p_block, rpgfa, &
rpgfb, sphi_a, sphi_b, zeta, zetb
REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(4, nexp_max) :: cval_ppl
REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
TYPE(sgp_potential_type), POINTER :: sgp_potential

!$ INTEGER(kind=omp_lock_kind), &
@@ -160,7 +161,6 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
END DO
END IF
END IF
ALLOCATE (force_thread(3, natom))
force_thread = 0.0_dp

maxder = ncoset(nder)
@@ -528,7 +528,6 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
!$OMP END DO
DEALLOCATE (atom_of_kind, kind_of)
END IF
DEALLOCATE (force_thread)

IF (calculate_forces .AND. use_virial) THEN
virial%pv_virial = virial%pv_virial + pv_thread
@@ -582,18 +581,19 @@ SUBROUTINE build_core_ppl_ri(lri_ppl_coef, force, virial, calculate_forces, use_
LOGICAL :: ecp_local, lpotextended
REAL(KIND=dp) :: alpha, dac, ppl_radius
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: va, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dva, dvas, force_thread
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dva, dvas
REAL(KIND=dp), DIMENSION(1:10) :: aloc, bloc
REAL(KIND=dp), DIMENSION(3) :: force_a, rac
REAL(KIND=dp), DIMENSION(3, 3) :: pv_thread
REAL(KIND=dp), DIMENSION(4, nexp_max) :: cval_ppl
TYPE(gto_basis_set_type), POINTER :: basis_set
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gth_potential_type), POINTER :: gth_potential
REAL(KIND=dp), DIMENSION(nexp_max) :: alpha_ppl
REAL(KIND=dp), DIMENSION(:, :), POINTER :: bcon, cval_lpot, rpgfa, sphi_a, zeta
REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
set_radius_a
REAL(KIND=dp), DIMENSION(:, :), POINTER :: bcon, cval_lpot, rpgfa, sphi_a, zeta
REAL(KIND=dp), DIMENSION(nexp_max) :: alpha_ppl
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gto_basis_set_type), POINTER :: basis_set
REAL(KIND=dp), DIMENSION(4, nexp_max) :: cval_ppl
REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
TYPE(sgp_potential_type), POINTER :: sgp_potential

!$ INTEGER(kind=omp_lock_kind), &
@@ -610,7 +610,6 @@ SUBROUTINE build_core_ppl_ri(lri_ppl_coef, force, virial, calculate_forces, use_
nkind = SIZE(atomic_kind_set)
natom = SIZE(particle_set)

ALLOCATE (force_thread(3, natom))
force_thread = 0.0_dp
pv_thread = 0.0_dp
ALLOCATE (atom_of_kind(natom), kind_of(natom))
@@ -792,7 +791,6 @@ SUBROUTINE build_core_ppl_ri(lri_ppl_coef, force, virial, calculate_forces, use_
force(ikind)%gth_ppl(3, atom_a) = force(ikind)%gth_ppl(3, atom_a) + force_thread(3, iatom)
END DO
END IF
DEALLOCATE (force_thread)
DEALLOCATE (atom_of_kind, kind_of)

IF (calculate_forces .AND. use_virial) THEN
@@ -112,21 +112,22 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
LOGICAL, DIMENSION(0:9) :: is_nonlocal
REAL(KIND=dp) :: dac, f0, ppnl_radius
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: force_thread, sab, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work
REAL(KIND=dp), DIMENSION(1) :: rprjc, zetc
REAL(KIND=dp), DIMENSION(3) :: fa, fb, rab, rac, rbc
REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc, pv_thread
REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
TYPE(clist_type), POINTER :: clist
TYPE(alist_type), POINTER :: alist_ac, alist_bc
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint, c_nl
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, h_block, h_nl, p_block, rpgfa, &
sphi_a, vprj_ppnl, zeta
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint, c_nl
TYPE(alist_type), POINTER :: alist_ac, alist_bc
TYPE(clist_type), POINTER :: clist
TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a
REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential
TYPE(sgp_potential_type), POINTER :: sgp_potential
@@ -422,7 +423,6 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
CALL sap_sort(sap_int)
! *** All integrals needed have been calculated and stored in sap_int
! *** We now calculate the Hamiltonian matrix elements
ALLOCATE (force_thread(3, natom))
force_thread = 0.0_dp
pv_thread = 0.0_dp

@@ -590,7 +590,6 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
!$OMP END DO
DEALLOCATE (atom_of_kind, kind_of)
END IF
DEALLOCATE (force_thread)

IF (calculate_forces .AND. use_virial) THEN
virial%pv_virial = virial%pv_virial + pv_thread
@@ -144,14 +144,62 @@ SUBROUTINE integrate_v_rspace(v_rspace, hmat, hmat_kp, pmat, pmat_kp, &
TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external

CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace', &
routineP = moduleN//':'//routineN
INTEGER :: natom

CALL get_qs_env(qs_env, natom=natom)
CALL integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env, &
calculate_forces, force_adm, ispin, compute_tau, gapw, &
basis_type, pw_env_external, task_list_external, natom)

END SUBROUTINE integrate_v_rspace

! **************************************************************************************************
!> \brief Implementation of integrate_v_rspace with the additional natom argument passed in to
!> allow for explicitly shaped force_thread array which is needed for OMP REDUCTION.
!> \param v_rspace ...
!> \param hmat ...
!> \param hmat_kp ...
!> \param pmat ...
!> \param pmat_kp ...
!> \param qs_env ...
!> \param calculate_forces ...
!> \param force_adm ...
!> \param ispin ...
!> \param compute_tau ...
!> \param gapw ...
!> \param basis_type ...
!> \param pw_env_external ...
!> \param task_list_external ...
!> \param natom ...
! **************************************************************************************************
SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env, &
calculate_forces, force_adm, ispin, compute_tau, gapw, &
basis_type, pw_env_external, task_list_external, natom)

TYPE(pw_p_type) :: v_rspace
TYPE(dbcsr_p_type), INTENT(INOUT), OPTIONAL :: hmat
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: hmat_kp
TYPE(dbcsr_p_type), INTENT(IN), OPTIONAL :: pmat
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: pmat_kp
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL, INTENT(IN) :: calculate_forces
LOGICAL, INTENT(IN), OPTIONAL :: force_adm
INTEGER, INTENT(IN), OPTIONAL :: ispin
LOGICAL, INTENT(IN), OPTIONAL :: compute_tau, gapw
CHARACTER(len=*), INTENT(IN), OPTIONAL :: basis_type
TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external
INTEGER, INTENT(IN) :: natom

CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace_low'

CHARACTER(len=default_string_length) :: my_basis_type
INTEGER :: atom_a, bcol, brow, handle, i, iatom, igrid_level, ikind, ikind_old, ilevel, img, &
ipair, ipgf, ipgf_new, iset, iset_new, iset_old, itask, ithread, jatom, jkind, jkind_old, &
jpgf, jpgf_new, jset, jset_new, jset_old, maxco, maxpgf, maxset, maxsgf_set, na1, na2, &
natom, nb1, nb2, ncoa, ncob, nimages, nkind, nseta, nsetb, nthread, offs_dv, sgfa, sgfb
nb1, nb2, ncoa, ncob, nimages, nkind, nseta, nsetb, nthread, offs_dv, sgfa, sgfb
INTEGER(KIND=int_8), DIMENSION(:), POINTER :: atom_pair_recv, atom_pair_send
INTEGER(kind=int_8), DIMENSION(:, :), POINTER :: tasks
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
@@ -164,9 +212,9 @@ SUBROUTINE integrate_v_rspace(v_rspace, hmat, hmat_kp, pmat, pmat_kp, &
p_duplicated, pab_required, scatter, use_subpatch, use_virial
REAL(KIND=dp) :: admm_scal_fac, dab, eps_gvg_rspace, &
rab2, scalef, zetp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: force_thread
REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra, rab, rab_inv, rb
REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b, pv_thread
REAL(KIND=dp), DIMENSION(3, natom) :: force_thread
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: dist_ab, h_block, hab, p_block, pab, &
rpgfa, rpgfb, sphi_a, sphi_b, work, &
@@ -301,10 +349,8 @@ SUBROUTINE integrate_v_rspace(v_rspace, hmat, hmat_kp, pmat, pmat_kp, &
CPASSERT(do_kp)
END IF
nkind = SIZE(qs_kind_set)
natom = SIZE(particle_set)
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
pv_thread = 0.0_dp
ALLOCATE (force_thread(3, natom))
force_thread = 0.0_dp

map_consistent = dft_control%qs_control%map_consistent
@@ -393,7 +439,6 @@ SUBROUTINE integrate_v_rspace(v_rspace, hmat, hmat_kp, pmat, pmat_kp, &
NULLIFY (hdabt, hadbt, hdab, hadb)

! get maximum numbers
natom = SIZE(particle_set)
maxset = 0
maxpgf = 0
DO ikind = 1, nkind
@@ -761,7 +806,6 @@ SUBROUTINE integrate_v_rspace(v_rspace, hmat, hmat_kp, pmat, pmat_kp, &
!$OMP END DO
DEALLOCATE (atom_of_kind, kind_of)
END IF
DEALLOCATE (force_thread)
IF (use_virial) THEN
virial%pv_virial = virial%pv_virial + pv_thread
END IF
@@ -818,6 +862,6 @@ SUBROUTINE integrate_v_rspace(v_rspace, hmat, hmat_kp, pmat, pmat_kp, &

CALL timestop(handle)

END SUBROUTINE integrate_v_rspace
END SUBROUTINE integrate_v_rspace_low

END MODULE qs_integrate_potential_product

0 comments on commit 597fd47

Please sign in to comment.
You can’t perform that action at this time.