From 1df1199c145652d4a172ffedffd9b9e59ce98d27 Mon Sep 17 00:00:00 2001 From: Frederick Stein <43850145+fstein93@users.noreply.github.com> Date: Sat, 11 Nov 2023 16:30:44 +0100 Subject: [PATCH] PW: Remove module pw/fast in favor of more generic copy routines (#3105) --- src/CMakeLists.txt | 1 - src/pw/dct.F | 44 +++----- src/pw/dgs.F | 110 ++++++++++++++----- src/pw/fast.F | 193 ---------------------------------- src/pw/fft_tools.F | 6 +- src/pw/pw_gpu.F | 8 +- src/pw/pw_methods.F | 188 ++++++++++++++++++++++++++++++--- src/pw/realspace_grid_types.F | 16 +-- src/qs_collocate_density.F | 4 +- 9 files changed, 282 insertions(+), 288 deletions(-) delete mode 100644 src/pw/fast.F diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ae69571e5c..4c00283501 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1131,7 +1131,6 @@ list( pw_env_methods.F pw_env/pw_env_types.F pw_env/rs_pw_interface.F - pw/fast.F pw/fft/fft_kinds.F pw/fft/fft_lib.F pw/fft/fft_plan.F diff --git a/src/pw/dct.F b/src/pw/dct.F index ba7d8db1f0..93f37247f3 100644 --- a/src/pw/dct.F +++ b/src/pw/dct.F @@ -15,7 +15,6 @@ ! ************************************************************************************************** MODULE dct - USE fast, ONLY: copy_cr USE kinds, ONLY: dp USE message_passing, ONLY: mp_comm_type,& mp_request_null,& @@ -24,9 +23,8 @@ MODULE dct USE pw_grid_types, ONLY: pw_grid_type USE pw_grids, ONLY: pw_grid_create,& pw_grid_setup - USE pw_types, ONLY: COMPLEXDATA1D,& - REALDATA1D,& - RECIPROCALSPACE,& + USE pw_types, ONLY: REALDATA3D,& + REALSPACE,& pw_type #include "../base/base_uses.f90" @@ -708,21 +706,19 @@ CONTAINS INTEGER, INTENT(IN) :: neumann_directions INTEGER, DIMENSION(:), INTENT(IN), POINTER :: dests_shrink - INTEGER, INTENT(IN) :: srcs_shrink + INTEGER, INTENT(INOUT) :: srcs_shrink INTEGER, DIMENSION(2, 3), INTENT(IN) :: bounds_local_shftd TYPE(pw_type), INTENT(IN) :: pw_in, pw_shrinked CHARACTER(LEN=*), PARAMETER :: routineN = 'pw_shrink' - COMPLEX(dp), DIMENSION(:, :, :), POINTER :: cc3d INTEGER :: group_size, handle, i, in_space, in_use, lb1_orig, lb1_xpnd, lb2_orig, lb2_xpnd, & lb3_orig, lb3_xpnd, maxn_sendrecv, rs_mpo, send_lb1, send_lb2, send_lb3, send_ub1, & - send_ub2, send_ub3, ub1_orig, ub1_xpnd, ub2_orig, ub2_xpnd, ub3_orig, ub3_xpnd + send_ub2, send_ub3, tag, ub1_orig, ub1_xpnd, ub2_orig, ub2_xpnd, ub3_orig, ub3_xpnd INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: bounds_local_all INTEGER, DIMENSION(2, 3) :: bounds_local_xpnd REAL(dp), DIMENSION(:, :, :), POINTER :: cr3d, send_crmsg TYPE(mp_comm_type) :: rs_group - TYPE(mp_request_type) :: recv_req, send_req TYPE(pw_grid_type), POINTER :: pw_grid_orig CALL timeset(routineN, handle) @@ -734,6 +730,9 @@ CONTAINS bounds_local_xpnd = pw_in%pw_grid%bounds_local in_space = pw_in%in_space in_use = pw_in%in_use + CPASSERT(in_space == REALSPACE) + CPASSERT(in_use == REALDATA3D) + tag = 1 SELECT CASE (neumann_directions) CASE (neumannXYZ, neumannXY) @@ -745,32 +744,15 @@ CONTAINS END SELECT ! cosine transform is a real transform. The cosine transform of a 3D data must be real and 3D. - NULLIFY (cr3d, cc3d) + NULLIFY (cr3d) lb1_xpnd = bounds_local_xpnd(1, 1) ub1_xpnd = bounds_local_xpnd(2, 1) lb2_xpnd = bounds_local_xpnd(1, 2) ub2_xpnd = bounds_local_xpnd(2, 2) lb3_xpnd = bounds_local_xpnd(1, 3) ub3_xpnd = bounds_local_xpnd(2, 3) - IF ((in_use .EQ. REALDATA1D) .OR. (in_use .EQ. COMPLEXDATA1D)) THEN - IF (in_space .EQ. RECIPROCALSPACE) THEN - ALLOCATE (cr3d(lb1_xpnd:ub1_xpnd, lb2_xpnd:ub2_xpnd, lb3_xpnd:ub3_xpnd)) - ALLOCATE (cc3d(lb1_xpnd:ub1_xpnd, lb2_xpnd:ub2_xpnd, lb3_xpnd:ub3_xpnd)) - cc3d = RESHAPE(pw_in%cc, (/ub1_xpnd - lb1_xpnd + 1, ub2_xpnd - lb2_xpnd + 1, ub3_xpnd - lb3_xpnd + 1/)) - CALL copy_cr(cc3d, cr3d) - DEALLOCATE (cc3d) - ELSE - ALLOCATE (cr3d(lb1_xpnd:ub1_xpnd, lb2_xpnd:ub2_xpnd, lb3_xpnd:ub3_xpnd)) - cr3d = RESHAPE(pw_in%cr, (/ub1_xpnd - lb1_xpnd + 1, ub2_xpnd - lb2_xpnd + 1, ub3_xpnd - lb3_xpnd + 1/)) - END IF - ELSE - IF (in_space .EQ. RECIPROCALSPACE) THEN - CALL copy_cr(pw_in%cc3d, cr3d) - ELSE - ALLOCATE (cr3d(lb1_xpnd:ub1_xpnd, lb2_xpnd:ub2_xpnd, lb3_xpnd:ub3_xpnd)) - cr3d = pw_in%cr3d - END IF - END IF + ALLOCATE (cr3d(lb1_xpnd:ub1_xpnd, lb2_xpnd:ub2_xpnd, lb3_xpnd:ub3_xpnd)) + cr3d = pw_in%cr3d ! let all the nodes know about each others shifted local bounds ALLOCATE (bounds_local_all(2, 3, group_size)) @@ -788,8 +770,7 @@ CONTAINS ALLOCATE (send_crmsg(send_lb1:send_ub1, send_lb2:send_ub2, send_lb3:send_ub3)) send_crmsg = cr3d(send_lb1:send_ub1, send_lb2:send_ub2, send_lb3:send_ub3) - CALL rs_group%isend(send_crmsg, dests_shrink(i), send_req) - CALL send_req%wait() + CALL rs_group%send(send_crmsg, dests_shrink(i), tag) DEALLOCATE (send_crmsg) END IF END DO @@ -807,8 +788,7 @@ CONTAINS ELSE IF (srcs_shrink .EQ. -1) THEN ! the source is invalid ... do nothing ELSE - CALL rs_group%irecv(pw_shrinked%cr3d, srcs_shrink, recv_req) - CALL recv_req%wait() + CALL rs_group%recv(pw_shrinked%cr3d, srcs_shrink, tag) END IF DEALLOCATE (bounds_local_all) diff --git a/src/pw/dgs.F b/src/pw/dgs.F index 3250229731..169cb288d6 100644 --- a/src/pw/dgs.F +++ b/src/pw/dgs.F @@ -12,17 +12,14 @@ ! ************************************************************************************************** MODULE dgs - USE fast, ONLY: copy_cr,& - copy_cri,& - rankup,& - vc_x_vc,& - vr_x_vc USE fft_tools, ONLY: BWFFT,& FFT_RADIX_CLOSEST,& FFT_RADIX_NEXT_ODD,& fft3d,& fft_radix_operations USE kinds, ONLY: dp + USE mathconstants, ONLY: z_one,& + z_zero USE mathlib, ONLY: det_3x3,& inv_3x3 USE pw_grid_info, ONLY: pw_find_cutoff @@ -30,8 +27,12 @@ MODULE dgs pw_grid_type USE pw_grids, ONLY: pw_grid_change,& pw_grid_setup + USE pw_methods, ONLY: pw_copy,& + pw_multiply_with,& + pw_zero USE pw_types, ONLY: COMPLEXDATA3D,& - REALDATA3D,& + pw_create,& + pw_release,& pw_type USE realspace_grid_types, ONLY: realspace_grid_type USE structure_factors, ONLY: structure_factor_evaluate @@ -861,28 +862,25 @@ CONTAINS COMPLEX(KIND=dp) :: za, zb COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zs - COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: cd INTEGER :: nd(3) + TYPE(pw_type) :: cd nd = rhos1%pw_grid%npts ALLOCATE (zs(nd(1)*nd(2))) zs = 0.0_dp - ALLOCATE (cd(nd(1), nd(2), nd(3))) - cd = 0.0_dp + CALL pw_create(cd, rho0%pw_grid, COMPLEXDATA3D, rho0%in_space) + CALL pw_zero(cd) za = CMPLX(0.0_dp, 0.0_dp, KIND=dp) zb = CMPLX(charge1, 0.0_dp, KIND=dp) - CALL rankup(nd, za, cd, zb, ex1, ey1, ez1, zs) - IF (rho0%in_use == REALDATA3D) & - CALL vr_x_vc(rho0%cr3d, cd) - IF (rho0%in_use == COMPLEXDATA3D) & - CALL vc_x_vc(rho0%cc3d, cd) - CALL fft3d(BWFFT, nd, cd) - CALL copy_cr(cd, rhos1%cr3d) + CALL rankup(nd, za, cd%cc3d, zb, ex1, ey1, ez1, zs) + CALL pw_multiply_with(cd, rho0) + CALL fft3d(BWFFT, nd, cd%cc3d) + CALL pw_copy(cd, rhos1) DEALLOCATE (zs) - DEALLOCATE (cd) + CALL pw_release(cd) END SUBROUTINE dg_get_patch_1 @@ -910,31 +908,28 @@ CONTAINS COMPLEX(KIND=dp) :: za, zb COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zs - COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: cd INTEGER :: nd(3) + TYPE(pw_type) :: cd nd = rhos1%pw_grid%npts ALLOCATE (zs(nd(1)*nd(2))) zs = 0.0_dp - ALLOCATE (cd(nd(1), nd(2), nd(3))) - cd = 0.0_dp + CALL pw_create(cd, rhos1%pw_grid, COMPLEXDATA3D, rho0%in_space) + CALL pw_zero(cd) za = CMPLX(0.0_dp, 0.0_dp, KIND=dp) zb = CMPLX(charge2, 0.0_dp, KIND=dp) - CALL rankup(nd, za, cd, zb, ex2, ey2, ez2, zs) + CALL rankup(nd, za, cd%cc3d, zb, ex2, ey2, ez2, zs) za = CMPLX(0.0_dp, 1.0_dp, KIND=dp) zb = CMPLX(charge1, 0.0_dp, KIND=dp) - CALL rankup(nd, za, cd, zb, ex1, ey1, ez1, zs) - IF (rho0%in_use == REALDATA3D) & - CALL vr_x_vc(rho0%cr3d, cd) - IF (rho0%in_use == COMPLEXDATA3D) & - CALL vc_x_vc(rho0%cc3d, cd) - CALL fft3d(BWFFT, nd, cd) - CALL copy_cri(cd, rhos1%cr3d, rhos2%cr3d) + CALL rankup(nd, za, cd%cc3d, zb, ex1, ey1, ez1, zs) + CALL pw_multiply_with(cd, rho0) + CALL fft3d(BWFFT, nd, cd%cc3d) + CALL copy_cri(cd%cc3d, rhos1%cr3d, rhos2%cr3d) DEALLOCATE (zs) - DEALLOCATE (cd) + CALL pw_release(cd) END SUBROUTINE dg_get_patch_2 @@ -1139,4 +1134,61 @@ CONTAINS END SUBROUTINE dg_int_patch_folded_1d +! ************************************************************************************************** +!> \brief ... +!> \param n ... +!> \param za ... +!> \param cmat ... +!> \param zb ... +!> \param ex ... +!> \param ey ... +!> \param ez ... +!> \param scr ... +! ************************************************************************************************** + SUBROUTINE rankup(n, za, cmat, zb, ex, ey, ez, scr) +! +! cmat(i,j,k) <- za * cmat(i,j,k) + ex(i) * ey(j) * ez(k) +! + + INTEGER, DIMENSION(3), INTENT(IN) :: n + COMPLEX(KIND=dp), INTENT(IN) :: za + COMPLEX(KIND=dp), DIMENSION(:, :, :), & + INTENT(INOUT) :: cmat + COMPLEX(KIND=dp), INTENT(IN) :: zb + COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN) :: ex, ey, ez + COMPLEX(KIND=dp), DIMENSION(:), INTENT(INOUT) :: scr + + INTEGER :: n2, n3 + + n2 = n(1)*n(2) + n3 = n2*n(3) + scr(1:n2) = z_zero + CALL zgeru(n(1), n(2), zb, ex, 1, ey, 1, scr, n(1)) + CALL zscal(n3, za, cmat, 1) + CALL zgeru(n2, n(3), z_one, scr, 1, ez, 1, cmat, n2) + + END SUBROUTINE rankup + +! ************************************************************************************************** +!> \brief Copy a the real and imag. parts of a complex 3D array into two real arrays +!> \param z the complex array +!> \param r1 the real array for the real part +!> \param r2 the real array for the imaginary part +! ************************************************************************************************** + SUBROUTINE copy_cri(z, r1, r2) +! +! r1 = real ( z ) +! r2 = imag ( z ) +! + + COMPLEX(KIND=dp), INTENT(IN) :: z(:, :, :) + REAL(KIND=dp), INTENT(INOUT) :: r1(:, :, :), r2(:, :, :) + +!$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(r1,r2,z) + r1(:, :, :) = REAL(z(:, :, :), KIND=dp) + r2(:, :, :) = AIMAG(z(:, :, :)) +!$OMP END PARALLEL WORKSHARE + + END SUBROUTINE copy_cri + END MODULE dgs diff --git a/src/pw/fast.F b/src/pw/fast.F deleted file mode 100644 index f141ef835f..0000000000 --- a/src/pw/fast.F +++ /dev/null @@ -1,193 +0,0 @@ -!--------------------------------------------------------------------------------------------------! -! CP2K: A general program to perform molecular dynamics simulations ! -! Copyright 2000-2023 CP2K developers group ! -! ! -! SPDX-License-Identifier: GPL-2.0-or-later ! -!--------------------------------------------------------------------------------------------------! - -MODULE fast - - USE kinds, ONLY: dp - USE mathconstants, ONLY: z_one,& - z_zero,& - zero - - IMPLICIT NONE - - PRIVATE - - PUBLIC :: rankup, vc_x_vc, vr_x_vc, copy_cri, copy_cr, copy_rc, zero_c - - INTERFACE zero_c - MODULE PROCEDURE zero_c2, zero_c3 - END INTERFACE - -CONTAINS - -! ************************************************************************************************** -!> \brief ... -!> \param n ... -!> \param za ... -!> \param cmat ... -!> \param zb ... -!> \param ex ... -!> \param ey ... -!> \param ez ... -!> \param scr ... -! ************************************************************************************************** - SUBROUTINE rankup(n, za, cmat, zb, ex, ey, ez, scr) -! -! cmat(i,j,k) <- za * cmat(i,j,k) + ex(i) * ey(j) * ez(k) -! - - INTEGER, DIMENSION(3), INTENT(IN) :: n - COMPLEX(KIND=dp), INTENT(IN) :: za - COMPLEX(KIND=dp), DIMENSION(:, :, :), & - INTENT(INOUT) :: cmat - COMPLEX(KIND=dp), INTENT(IN) :: zb - COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN) :: ex, ey, ez - COMPLEX(KIND=dp), DIMENSION(:), INTENT(INOUT) :: scr - - INTEGER :: n2, n3 - - n2 = n(1)*n(2) - n3 = n2*n(3) - scr(1:n2) = z_zero - CALL zgeru(n(1), n(2), zb, ex, 1, ey, 1, scr, n(1)) - CALL zscal(n3, za, cmat, 1) - CALL zgeru(n2, n(3), z_one, scr, 1, ez, 1, cmat, n2) - - END SUBROUTINE rankup - -! ************************************************************************************************** -!> \brief Multiply two complex 3D arrays element-wise -!> \param cvec2 the other complex array -!> \param cvec the complex array returned as the result -! ************************************************************************************************** - SUBROUTINE vc_x_vc(cvec2, cvec) -! -! cvec(i) <- cvec(i) * cvec2(i) -! - - COMPLEX(KIND=dp), INTENT(IN) :: cvec2(:, :, :) - COMPLEX(KIND=dp), INTENT(INOUT) :: cvec(:, :, :) - -!$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(cvec,cvec2) - cvec(:, :, :) = cvec(:, :, :)*cvec2(:, :, :) -!$OMP END PARALLEL WORKSHARE - - END SUBROUTINE vc_x_vc - -! ************************************************************************************************** -!> \brief Scale a complex array element-wise by a real array -!> \param rvec the real array -!> \param cvec the complex array -! ************************************************************************************************** - SUBROUTINE vr_x_vc(rvec, cvec) -! -! cvec(i) <- cvec(i) * rvec(i) -! - - REAL(KIND=dp), INTENT(IN) :: rvec(:, :, :) - COMPLEX(KIND=dp), INTENT(INOUT) :: cvec(:, :, :) - -!$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(cvec,rvec) - cvec(:, :, :) = cvec(:, :, :)*CMPLX(rvec(:, :, :), KIND=dp) -!$OMP END PARALLEL WORKSHARE - - END SUBROUTINE vr_x_vc - -! ************************************************************************************************** -!> \brief Copy a the real and imag. parts of a complex 3D array into two real arrays -!> \param z the complex array -!> \param r1 the real array for the real part -!> \param r2 the real array for the imaginary part -! ************************************************************************************************** - SUBROUTINE copy_cri(z, r1, r2) -! -! r1 = real ( z ) -! r2 = imag ( z ) -! - - COMPLEX(KIND=dp), INTENT(IN) :: z(:, :, :) - REAL(KIND=dp), INTENT(INOUT) :: r1(:, :, :), r2(:, :, :) - -!$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(r1,r2,z) - r1(:, :, :) = REAL(z(:, :, :), KIND=dp) - r2(:, :, :) = AIMAG(z(:, :, :)) -!$OMP END PARALLEL WORKSHARE - - END SUBROUTINE copy_cri - -! ************************************************************************************************** -!> \brief Copy the real part of a complex 3D array to a real array -!> \param z the complex array -!> \param r1 the real array -! ************************************************************************************************** - SUBROUTINE copy_cr(z, r1) -! -! r1 = real ( z ) -! - - COMPLEX(KIND=dp), INTENT(IN) :: z(:, :, :) - REAL(KIND=dp), INTENT(INOUT) :: r1(:, :, :) - -!$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(r1,z) - r1(:, :, :) = REAL(z(:, :, :), KIND=dp) -!$OMP END PARALLEL WORKSHARE - - END SUBROUTINE copy_cr - -! ************************************************************************************************** -!> \brief Copy a real 3D array to complex 3D array -!> \param r1 the real array -!> \param z the complex array -! ************************************************************************************************** - SUBROUTINE copy_rc(r1, z) -! -! z = r1 -! - - REAL(KIND=dp), INTENT(IN) :: r1(:, :, :) - COMPLEX(KIND=dp), INTENT(INOUT) :: z(:, :, :) - -!$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(r1,z) - z(:, :, :) = CMPLX(r1(:, :, :), zero, KIND=dp) -!$OMP END PARALLEL WORKSHARE - - END SUBROUTINE copy_rc - -! ************************************************************************************************** -!> \brief Zero a complex 2D array (optionally with OpenMP) -!> \param z the array -! ************************************************************************************************** - SUBROUTINE zero_c2(z) -! -! z = ( 0.0_dp , 0.0_dp) -! - - COMPLEX(KIND=dp), INTENT(INOUT) :: z(:, :) - -!$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(z) - z(:, :) = z_zero -!$OMP END PARALLEL WORKSHARE - END SUBROUTINE zero_c2 - -! ************************************************************************************************** -!> \brief Zero a complex 3D array (optionally with OpenMP) -!> \param z the array -! ************************************************************************************************** - SUBROUTINE zero_c3(z) -! -! z = ( 0.0_dp , 0.0_dp) -! - - COMPLEX(KIND=dp), INTENT(INOUT) :: z(:, :, :) - -!$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(z) - z(:, :, :) = z_zero -!$OMP END PARALLEL WORKSHARE - - END SUBROUTINE zero_c3 - -END MODULE fast diff --git a/src/pw/fft_tools.F b/src/pw/fft_tools.F index 6f85203171..ab7921473d 100644 --- a/src/pw/fft_tools.F +++ b/src/pw/fft_tools.F @@ -31,7 +31,6 @@ MODULE fft_tools C_PTR,& C_SIZE_T USE cp_log_handling, ONLY: cp_logger_get_default_io_unit - USE fast, ONLY: zero_c USE fft_lib, ONLY: & fft_1dm, fft_3d, fft_alloc, fft_create_plan_1dm, fft_create_plan_3d, fft_dealloc, & fft_destroy_plan, fft_do_cleanup, fft_do_init, fft_get_lengths, fft_library @@ -39,6 +38,7 @@ MODULE fft_tools USE kinds, ONLY: dp,& dp_size,& sp + USE mathconstants, ONLY: z_zero USE message_passing, ONLY: mp_cart_type,& mp_comm_null,& mp_comm_type,& @@ -820,8 +820,8 @@ CONTAINS sbuf => fft_scratch%r1buf tbuf => fft_scratch%tbuf - CALL zero_c(sbuf) - CALL zero_c(tbuf) + sbuf = z_zero + tbuf = z_zero IF (sign == FWFFT) THEN ! cin -> gin diff --git a/src/pw/pw_gpu.F b/src/pw/pw_gpu.F index fd754fbf18..a923e46453 100644 --- a/src/pw/pw_gpu.F +++ b/src/pw/pw_gpu.F @@ -28,11 +28,11 @@ MODULE pw_gpu C_INT,& C_LOC,& C_PTR - USE fast, ONLY: zero_c USE fft_tools, ONLY: & cube_transpose_1, cube_transpose_2, fft_scratch_sizes, fft_scratch_type, get_fft_scratch, & release_fft_scratch, x_to_yz, xz_to_yz, yz_to_x, yz_to_xz USE kinds, ONLY: dp + USE mathconstants, ONLY: z_zero USE message_passing, ONLY: mp_cart_type,& mp_comm_type USE pw_grid_types, ONLY: FULLSPACE @@ -504,7 +504,7 @@ CONTAINS CALL pw_gpu_sf(pw1, pbuf, scale) ! Exchange data ( transpose of matrix ) and sort - IF (pw1%pw_grid%grid_span /= FULLSPACE) CALL zero_c(qbuf) + IF (pw1%pw_grid%grid_span /= FULLSPACE) qbuf = z_zero CALL yz_to_xz(pbuf, rs_group, r_dim, g_pos, p2p, pw1%pw_grid%para%yzp, nyzray, & bo(:, :, :, 2), qbuf, fft_scratch) @@ -515,7 +515,7 @@ CONTAINS CALL pw_gpu_f(qbuf, rbuf, -1, n(2), mx2*mz2) ! Exchange data ( transpose of matrix ) - IF (pw1%pw_grid%grid_span /= FULLSPACE) CALL zero_c(sbuf) + IF (pw1%pw_grid%grid_span /= FULLSPACE) sbuf = z_zero CALL cube_transpose_1(rbuf, bo(:, :, :, 2), bo(:, :, :, 1), sbuf, fft_scratch) @@ -541,7 +541,7 @@ CONTAINS CALL pw_gpu_sf(pw1, sbuf, scale) ! Exchange data ( transpose of matrix ) and sort - IF (pw1%pw_grid%grid_span /= FULLSPACE) CALL zero_c(tbuf) + IF (pw1%pw_grid%grid_span /= FULLSPACE) tbuf = z_zero CALL x_to_yz(sbuf, gs_group, g_pos, p2p, pw1%pw_grid%para%yzp, nyzray, & bo(:, :, :, 2), tbuf, fft_scratch) diff --git a/src/pw/pw_methods.F b/src/pw/pw_methods.F index 682674d3c6..29a201afc0 100644 --- a/src/pw/pw_methods.F +++ b/src/pw/pw_methods.F @@ -26,9 +26,6 @@ MODULE pw_methods USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,& cp_to_string - USE fast, ONLY: copy_cr,& - copy_rc,& - zero_c USE fft_tools, ONLY: BWFFT,& FWFFT,& fft3d @@ -36,6 +33,7 @@ MODULE pw_methods accurate_sum USE kinds, ONLY: dp USE machine, ONLY: m_memory + USE mathconstants, ONLY: z_zero USE pw_copy_all, ONLY: pw_copy_match USE pw_fpga, ONLY: pw_fpga_c1dr3d_3d_dp,& pw_fpga_c1dr3d_3d_sp,& @@ -71,6 +69,7 @@ MODULE pw_methods PUBLIC :: pw_integral_ab, pw_integral_a2b PUBLIC :: pw_dr2_gg, pw_integrate_function PUBLIC :: pw_set, pw_truncated + PUBLIC :: pw_copy_to_array, pw_copy_from_array CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_methods' LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. @@ -89,6 +88,16 @@ MODULE pw_methods MODULE PROCEDURE pw_set_value, pw_set_array END INTERFACE + INTERFACE pw_copy_to_array + MODULE PROCEDURE pw_copy_to_array_r + MODULE PROCEDURE pw_copy_to_array_c + END INTERFACE + + INTERFACE pw_copy_from_array + MODULE PROCEDURE pw_copy_from_array_r + MODULE PROCEDURE pw_copy_from_array_c + END INTERFACE + CONTAINS ! ************************************************************************************************** @@ -266,6 +275,16 @@ CONTAINS pw1%in_space == REALSPACE) THEN !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw1,pw2) pw2%cr3d(:, :, :) = pw1%cr3d(:, :, :) +!$OMP END PARALLEL WORKSHARE + ELSE IF (pw1%in_use == COMPLEXDATA3D .AND. pw2%in_use == REALDATA3D .AND. & + pw1%in_space == REALSPACE) THEN +!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw1,pw2) + pw2%cr3d(:, :, :) = REAL(pw1%cc3d(:, :, :), KIND=dp) +!$OMP END PARALLEL WORKSHARE + ELSE IF (pw1%in_use == REALDATA3D .AND. pw2%in_use == COMPLEXDATA3D .AND. & + pw1%in_space == REALSPACE) THEN +!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw1,pw2) + pw2%cc3d(:, :, :) = pw1%cr3d(:, :, :) !$OMP END PARALLEL WORKSHARE ELSE CPABORT("No suitable data field") @@ -319,6 +338,22 @@ CONTAINS !$OMP WORKSHARE pw2%cc3d(:, :, :) = pw1%cc3d(:, :, :) !$OMP END WORKSHARE +!$OMP END PARALLEL + ELSE IF (pw1%in_use == REALDATA3D .AND. & + pw2%in_use == COMPLEXDATA3D) THEN + ns = SIZE(pw1%cc3d) +!$OMP PARALLEL DEFAULT(NONE) SHARED(pw1, pw2) +!$OMP WORKSHARE + pw2%cc3d(:, :, :) = pw1%cr3d(:, :, :) +!$OMP END WORKSHARE +!$OMP END PARALLEL + ELSE IF (pw1%in_use == COMPLEXDATA3D .AND. & + pw2%in_use == REALDATA3D) THEN + ns = SIZE(pw1%cc3d) +!$OMP PARALLEL DEFAULT(NONE) SHARED(pw1, pw2) +!$OMP WORKSHARE + pw2%cr3d(:, :, :) = REAL(pw1%cc3d(:, :, :), KIND=dp) +!$OMP END WORKSHARE !$OMP END PARALLEL ELSE CPABORT("No suitable data field") @@ -332,6 +367,127 @@ CONTAINS END SUBROUTINE pw_copy +! ************************************************************************************************** +!> \brief ... +!> \param pw ... +!> \param array ... +! ************************************************************************************************** + SUBROUTINE pw_copy_to_array_r(pw, array) + TYPE(pw_type), INTENT(IN) :: pw + REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: array + + CHARACTER(len=*), PARAMETER :: routineN = 'pw_copy_to_array_r' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + IF (pw%in_use == REALDATA3D) THEN +!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array) + array = pw%cr3d +!$OMP END PARALLEL WORKSHARE + ELSE IF (pw%in_use == COMPLEXDATA3D) THEN +!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array) + array = REAL(pw%cc3d, KIND=dp) +!$OMP END PARALLEL WORKSHARE + ELSE + CPABORT("PW grid must be 3D!") + END IF + + CALL timestop(handle) + END SUBROUTINE pw_copy_to_array_r + +! ************************************************************************************************** +!> \brief ... +!> \param pw ... +!> \param array ... +! ************************************************************************************************** + SUBROUTINE pw_copy_to_array_c(pw, array) + TYPE(pw_type), INTENT(IN) :: pw + COMPLEX(KIND=dp), DIMENSION(:, :, :), & + INTENT(INOUT) :: array + + CHARACTER(len=*), PARAMETER :: routineN = 'pw_copy_to_array_c' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + IF (pw%in_use == REALDATA3D) THEN +!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array) + array = pw%cr3d +!$OMP END PARALLEL WORKSHARE + ELSE IF (pw%in_use == COMPLEXDATA3D) THEN +!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array) + array = pw%cc3d +!$OMP END PARALLEL WORKSHARE + ELSE + CPABORT("PW grid must be 3D!") + END IF + + CALL timestop(handle) + END SUBROUTINE pw_copy_to_array_c + +! ************************************************************************************************** +!> \brief ... +!> \param pw ... +!> \param array ... +! ************************************************************************************************** + SUBROUTINE pw_copy_from_array_r(pw, array) + TYPE(pw_type), INTENT(IN) :: pw + REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: array + + CHARACTER(len=*), PARAMETER :: routineN = 'pw_copy_from_array_r' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + IF (pw%in_use == REALDATA3D) THEN +!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array) + pw%cr3d = array +!$OMP END PARALLEL WORKSHARE + ELSE IF (pw%in_use == COMPLEXDATA3D) THEN +!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array) + pw%cc3d = array +!$OMP END PARALLEL WORKSHARE + ELSE + CPABORT("PW grid must be 3D!") + END IF + + CALL timestop(handle) + END SUBROUTINE pw_copy_from_array_r + +! ************************************************************************************************** +!> \brief ... +!> \param pw ... +!> \param array ... +! ************************************************************************************************** + SUBROUTINE pw_copy_from_array_c(pw, array) + TYPE(pw_type), INTENT(IN) :: pw + COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: array + + CHARACTER(len=*), PARAMETER :: routineN = 'pw_copy_from_array_c' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + IF (pw%in_use == REALDATA3D) THEN +!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array) + pw%cr3d = REAL(array, KIND=dp) +!$OMP END PARALLEL WORKSHARE + ELSE IF (pw%in_use == COMPLEXDATA3D) THEN +!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(pw, array) + pw%cc3d = array +!$OMP END PARALLEL WORKSHARE + ELSE + CPABORT("PW grid must be 3D!") + END IF + + CALL timestop(handle) + END SUBROUTINE pw_copy_from_array_c + ! ************************************************************************************************** !> \brief multiplies pw coeffs with a number !> \param pw ... @@ -1731,7 +1887,7 @@ CONTAINS ASSOCIATE (mapl => pw%pw_grid%mapl%pos, mapm => pw%pw_grid%mapm%pos, mapn => pw%pw_grid%mapn%pos, & ghat => pw%pw_grid%g_hat, yzq => pw%pw_grid%para%yzq, ngpts => SIZE(pw%pw_grid%gsq)) - IF (.NOT. PRESENT(scale)) CALL zero_c(c) + IF (.NOT. PRESENT(scale)) c = z_zero IF (PRESENT(scale)) THEN !$OMP PARALLEL DO DEFAULT(NONE), & @@ -1922,7 +2078,7 @@ CONTAINS ! check if bitstream for the fft size is present ! if not, perform fft3d in CPU IF (pw_fpga_init_bitstream(n) == 1) THEN - CALL copy_rc(pw1%cr3d, c_out) + CALL pw_copy_to_array(pw1, c_out) #if (__PW_FPGA_SP && __PW_FPGA) CALL pw_fpga_r3dc1d_3d_sp(n, c_out) #else @@ -1931,14 +2087,14 @@ CONTAINS CALL zdscal(n(1)*n(2)*n(3), norm, c_out, 1) CALL pw_gather(pw2, c_out) ELSE - CALL copy_rc(pw1%cr3d, c_out) + CALL pw_copy_to_array(pw1, c_out) CALL fft3d(dir, n, c_out, scale=norm, debug=test) CALL pw_gather(pw2, c_out) END IF DEALLOCATE (c_out) #else ALLOCATE (c_out(n(1), n(2), n(3))) - CALL copy_rc(pw1%cr3d, c_out) + CALL pw_copy_to_array(pw1, c_out) CALL fft3d(dir, n, c_out, scale=norm, debug=test) CALL pw_gather(pw2, c_out) DEALLOCATE (c_out) @@ -1977,7 +2133,7 @@ CONTAINS #endif CALL zdscal(n(1)*n(2)*n(3), norm, c_out, 1) ! use real part only - CALL copy_cr(c_out, pw2%cr3d) + CALL pw_copy_from_array(pw2, c_out) ELSE IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " PW_SCATTER : 3d -> 1d " CALL pw_scatter(pw1, c_out) @@ -1985,7 +2141,7 @@ CONTAINS CALL fft3d(dir, n, c_out, scale=norm, debug=test) ! use real part only IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " REAL part " - CALL copy_cr(c_out, pw2%cr3d) + CALL pw_copy_from_array(pw2, c_out) END IF DEALLOCATE (c_out) #else @@ -1996,7 +2152,7 @@ CONTAINS CALL fft3d(dir, n, c_out, scale=norm, debug=test) ! use real part only IF (test .AND. out_unit > 0) WRITE (out_unit, '(A)') " REAL part " - CALL copy_cr(c_out, pw2%cr3d) + CALL pw_copy_from_array(pw2, c_out) DEALLOCATE (c_out) #endif END SELECT @@ -2038,7 +2194,7 @@ CONTAINS CASE ("FW_C3DC1D") !..prepare input c_in => pw1%cc3d - CALL zero_c(grays) + grays = z_zero !..transform IF (pw1%pw_grid%para%ray_distribution) THEN CALL fft3d(dir, n, c_in, grays, pw1%pw_grid%para%group, & @@ -2064,8 +2220,8 @@ CONTAINS !.. prepare input nloc = pw1%pw_grid%npts_local ALLOCATE (c_in(nloc(1), nloc(2), nloc(3))) - CALL copy_rc(pw1%cr3d, c_in) - CALL zero_c(grays) + CALL pw_copy_to_array(pw1, c_in) + grays = z_zero !..transform IF (pw1%pw_grid%para%ray_distribution) THEN CALL fft3d(dir, n, c_in, grays, pw1%pw_grid%para%group, & @@ -2089,7 +2245,7 @@ CONTAINS !..prepare input IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) & WRITE (out_unit, '(A)') " PW_SCATTER : 2d -> 1d " - CALL zero_c(grays) + grays = z_zero CALL pw_scatter(pw1, grays) c_in => pw2%cc3d !..transform @@ -2114,7 +2270,7 @@ CONTAINS !.. prepare input IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) & WRITE (out_unit, '(A)') " PW_SCATTER : 2d -> 1d " - CALL zero_c(grays) + grays = z_zero CALL pw_scatter(pw1, grays) nloc = pw2%pw_grid%npts_local ALLOCATE (c_in(nloc(1), nloc(2), nloc(3))) @@ -2131,7 +2287,7 @@ CONTAINS !..prepare output IF (test .AND. pw1%pw_grid%para%group_head .AND. out_unit > 0) & WRITE (out_unit, '(A)') " Real part " - CALL copy_cr(c_in, pw2%cr3d) + CALL pw_copy_from_array(pw2, c_in) DEALLOCATE (c_in) #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW) END IF diff --git a/src/pw/realspace_grid_types.F b/src/pw/realspace_grid_types.F index 54ba8519d7..a02ed1bb51 100644 --- a/src/pw/realspace_grid_types.F +++ b/src/pw/realspace_grid_types.F @@ -18,8 +18,6 @@ MODULE realspace_grid_types USE cp_array_utils, ONLY: cp_1d_r_p_type USE cp_log_handling, ONLY: cp_to_string - USE fast, ONLY: copy_cr,& - copy_rc USE kahan_sum, ONLY: accurate_sum USE kinds, ONLY: dp,& int_8 @@ -38,7 +36,9 @@ MODULE realspace_grid_types pw_grid_type USE pw_grids, ONLY: pw_grid_release,& pw_grid_retain - USE pw_methods, ONLY: pw_integrate_function + USE pw_methods, ONLY: pw_copy_from_array,& + pw_copy_to_array,& + pw_integrate_function USE pw_types, ONLY: COMPLEXDATA3D,& REALDATA3D,& pw_type @@ -659,7 +659,7 @@ CONTAINS SUBROUTINE transfer_rs2pw(rs, pw) TYPE(realspace_grid_type), INTENT(IN) :: rs - TYPE(pw_type), INTENT(IN) :: pw + TYPE(pw_type), INTENT(INOUT) :: pw CHARACTER(len=*), PARAMETER :: routineN = 'transfer_rs2pw' @@ -692,7 +692,7 @@ CONTAINS END IF ELSEIF (pw%in_use == COMPLEXDATA3D) THEN IF (rs%desc%border == 0) THEN - CALL copy_rc(rs%r, pw%cc3d) + CALL pw_copy_from_array(pw, rs%r) ELSE CPASSERT(LBOUND(pw%cr3d, 3) .EQ. rs%lb_real(3)) !$OMP PARALLEL DO DEFAULT(NONE) SHARED(pw,rs) @@ -779,7 +779,7 @@ CONTAINS END IF ELSEIF (pw%in_use == COMPLEXDATA3D) THEN IF (rs%desc%border == 0) THEN - CALL copy_cr(pw%cc3d, rs%r) + CALL pw_copy_to_array(pw, rs%r) ELSE !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP PRIVATE(i,im,j,jm,k,km) & @@ -831,7 +831,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE transfer_rs2pw_replicated(rs, pw) TYPE(realspace_grid_type), INTENT(IN) :: rs - TYPE(pw_type), INTENT(IN) :: pw + TYPE(pw_type), INTENT(INOUT) :: pw INTEGER :: dest, ii, ip, ix, iy, iz, nma, nn, s(3), & source @@ -888,7 +888,7 @@ CONTAINS IF (pw%in_use == REALDATA3D) THEN CALL dcopy(nn, sendbuf, 1, pw%cr3d, 1) ELSEIF (pw%in_use == COMPLEXDATA3D) THEN - CALL copy_rc(RESHAPE(sendbuf, SHAPE(pw%cc3d)), pw%cc3d) + CALL pw_copy_from_array(pw, RESHAPE(sendbuf, SHAPE(pw%cc3d))) ELSE CPABORT("PW type not compatible") END IF diff --git a/src/qs_collocate_density.F b/src/qs_collocate_density.F index fc822d3d26..ea1f175521 100644 --- a/src/qs_collocate_density.F +++ b/src/qs_collocate_density.F @@ -147,7 +147,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE calculate_rho_nlcc(rho_nlcc, qs_env) - TYPE(pw_type), INTENT(IN) :: rho_nlcc + TYPE(pw_type), INTENT(INOUT) :: rho_nlcc TYPE(qs_environment_type), POINTER :: qs_env CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_nlcc' @@ -326,7 +326,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE calculate_ppl_grid(vppl, qs_env) - TYPE(pw_type), INTENT(IN) :: vppl + TYPE(pw_type), INTENT(INOUT) :: vppl TYPE(qs_environment_type), POINTER :: qs_env CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ppl_grid'