From 73066b4ac54eeb99ece66947d7827b0482da4e12 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sat, 29 Jun 2024 01:15:47 +0200 Subject: [PATCH] issue with linking with CUDA --- plugins/local/tc_int/compute_tc_int_gpu.irp.f | 157 ++++++++++++++++++ plugins/local/tc_int/gpu.c | 2 + plugins/local/tc_int/gpu_module.F90 | 38 +++++ plugins/local/tc_int/write_tc_int_gpu.irp.f | 56 +++++++ src/dft_utils_in_r/ao_in_r.irp.f | 113 ++++++++----- 5 files changed, 329 insertions(+), 37 deletions(-) create mode 100644 plugins/local/tc_int/compute_tc_int_gpu.irp.f create mode 100644 plugins/local/tc_int/gpu.c create mode 100644 plugins/local/tc_int/gpu_module.F90 create mode 100644 plugins/local/tc_int/write_tc_int_gpu.irp.f diff --git a/plugins/local/tc_int/compute_tc_int_gpu.irp.f b/plugins/local/tc_int/compute_tc_int_gpu.irp.f new file mode 100644 index 00000000..146574a6 --- /dev/null +++ b/plugins/local/tc_int/compute_tc_int_gpu.irp.f @@ -0,0 +1,157 @@ + +! --- + +subroutine provide_int2_grad1_u12_ao_gpu() + + use gpu_module + + BEGIN_DOC + ! + ! int2_grad1_u12_ao(i,j,ipoint,1) = \int dr2 [\grad1 u(r1,r2)]_x1 \chi_i(r2) \chi_j(r2) + ! int2_grad1_u12_ao(i,j,ipoint,2) = \int dr2 [\grad1 u(r1,r2)]_y1 \chi_i(r2) \chi_j(r2) + ! int2_grad1_u12_ao(i,j,ipoint,3) = \int dr2 [\grad1 u(r1,r2)]_z1 \chi_i(r2) \chi_j(r2) + ! int2_grad1_u12_ao(i,j,ipoint,4) = \int dr2 [-(1/2) [\grad1 u(r1,r2)]^2] \chi_i(r2) \chi_j(r2) + ! + ! + ! tc_int_2e_ao(k,i,l,j) = (ki|V^TC(r_12)|lj) + ! = where V^TC(r_12) is the total TC operator + ! = tc_grad_and_lapl_ao(k,i,l,j) + tc_grad_square_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j) + ! where: + ! + ! tc_grad_and_lapl_ao(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) . \grad_1 | ij > + ! = -1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) + ! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 (-1) \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) + ! + ! tc_grad_square_ao(k,i,l,j) = -1/2 + ! + ! ao_two_e_coul(k,i,l,j) = < l k | 1/r12 | j i > = ( k i | 1/r12 | l j ) + ! + END_DOC + + implicit none + + integer :: i, j, k, l, m, ipoint, jpoint + integer :: n_blocks, n_rest, n_pass + integer :: i_blocks, i_rest, i_pass, ii + double precision :: mem, n_double + double precision :: weight1, ao_k_r, ao_i_r + double precision :: der_envsq_x, der_envsq_y, der_envsq_z, lap_envsq + double precision :: time0, time1, time2, tc1, tc2, tc + double precision, allocatable :: int2_grad1_u12_ao(:,:,:,:), tc_int_2e_ao(:,:,:,:) + double precision, allocatable :: tmp(:,:,:), c_mat(:,:,:), tmp_grad1_u12(:,:,:) + + double precision, external :: get_ao_two_e_integral + + + PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra + PROVIDE final_weight_at_r_vector aos_grad_in_r_array_transp_bis final_weight_at_r_vector aos_in_r_array_transp + + + + print*, ' start provide_int2_grad1_u12_ao_gpu ...' + call wall_time(time0) + + call total_memory(mem) + mem = max(1.d0, qp_max_mem - mem) + n_double = mem * 1.d8 + n_blocks = int(min(n_double / (n_points_extra_final_grid * 4.d0), 1.d0*n_points_final_grid)) + n_rest = int(mod(n_points_final_grid, n_blocks)) + n_pass = int((n_points_final_grid - n_rest) / n_blocks) + + call write_int(6, n_pass, 'Number of passes') + call write_int(6, n_blocks, 'Size of the blocks') + call write_int(6, n_rest, 'Size of the last block') + + ! --- + + allocate(int2_grad1_u12_ao(ao_num,ao_num,n_points_final_grid,4)) + allocate(tc_int_2e_ao(ao_num,ao_num,ao_num,ao_num)) + + double precision, allocatable :: aos_data1(:,:,:) + double precision, allocatable :: aos_data2(:,:,:) + allocate(aos_data1(n_points_final_grid,ao_num,4)) + allocate(aos_data2(n_points_extra_final_grid,ao_num,4)) + + do k = 1, ao_num + do ipoint = 1, n_points_final_grid + aos_data1(ipoint,k,1) = aos_in_r_array(i,ipoint) + aos_data1(ipoint,k,2) = aos_grad_in_r_array(i,ipoint,1) + aos_data1(ipoint,k,3) = aos_grad_in_r_array(i,ipoint,2) + aos_data1(ipoint,k,4) = aos_grad_in_r_array(i,ipoint,3) + enddo + + do ipoint = 1, n_points_extra_final_grid + aos_data1(ipoint,k,1) = aos_in_r_array_extra(i,ipoint) + aos_data1(ipoint,k,2) = aos_grad_in_r_array_extra(i,ipoint,1) + aos_data1(ipoint,k,3) = aos_grad_in_r_array_extra(i,ipoint,2) + aos_data1(ipoint,k,4) = aos_grad_in_r_array_extra(i,ipoint,3) + enddo + enddo + + call tc_int_bh(n_points_final_grid, n_points_extra_final_grid, ao_num, nucl_num, & + jBH_size, jBH_m, jBH_n, jBH_o, jBH_c, & + final_grid_points, final_grid_points_extra, nucl_coord, & + final_weight_at_r_vector, final_weight_at_r_vector_extra, & + aos_data1, aos_data2, int2_grad1_u12_ao, tc_int_2e_ao) + + ! --- + + call wall_time(time1) + + PROVIDE ao_integrals_map + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(ao_num, tc_int_2e_ao, ao_integrals_map) & + !$OMP PRIVATE(i, j, k, l) + !$OMP DO COLLAPSE(3) + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + ! < 1:i, 2:j | 1:k, 2:l > + tc_int_2e_ao(k,i,l,j) = tc_int_2e_ao(k,i,l,j) + get_ao_two_e_integral(i, j, k, l, ao_integrals_map) + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(time2) + print*, ' wall time of Coulomb part of tc_int_2e_ao (min) ', (time2 - time1) / 60.d0 + call print_memory_usage() + + ! --- + + print*, ' Writing int2_grad1_u12_ao in ', trim(ezfio_filename) // '/work/int2_grad1_u12_ao' + open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="write") + call ezfio_set_work_empty(.False.) + write(11) int2_grad1_u12_ao(:,:,:,1:3) + close(11) + + print*, ' Saving tc_int_2e_ao in ', trim(ezfio_filename) // '/work/ao_two_e_tc_tot' + open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/ao_two_e_tc_tot', action="write") + call ezfio_set_work_empty(.False.) + do i = 1, ao_num + write(11) tc_int_2e_ao(:,:,:,i) + enddo + close(11) + + ! ---- + + deallocate(int2_grad1_u12_ao) + deallocate(tc_int_2e_ao) + + call wall_time(time2) + print*, ' wall time for tc_int_2e_ao (min) = ', (time2-time1) / 60.d0 + call print_memory_usage() + + ! --- + + call wall_time(time1) + print*, ' wall time for TC-integrals (min) = ', (time1-time0) / 60.d0 + + return +end + +! --- + diff --git a/plugins/local/tc_int/gpu.c b/plugins/local/tc_int/gpu.c new file mode 100644 index 00000000..139597f9 --- /dev/null +++ b/plugins/local/tc_int/gpu.c @@ -0,0 +1,2 @@ + + diff --git a/plugins/local/tc_int/gpu_module.F90 b/plugins/local/tc_int/gpu_module.F90 new file mode 100644 index 00000000..cf7efad2 --- /dev/null +++ b/plugins/local/tc_int/gpu_module.F90 @@ -0,0 +1,38 @@ + +! --- + +module gpu_module + + use iso_c_binding + + implicit none + + interface + + subroutine tc_int_bh(n_grid1, n_grid2, ao_num, n_nuc, & + size_bh, m_bh, n_bh, o_bh, c_bh, & + r1, r2, rn, wr1, wr2, aos_data1, & + aos_data2, int2_grad1_u12, tc_int_2e_ao) bind(C) + + import c_int, c_double + + integer(c_int), intent(in), value :: n_grid1, n_grid2, ao_num, n_nuc, size_bh + integer(c_int), intent(in) :: m_bh(size_bh,n_nuc) + integer(c_int), intent(in) :: n_bh(size_bh,n_nuc) + integer(c_int), intent(in) :: o_bh(size_bh,n_nuc) + real(c_double), intent(in) :: c_bh(size_bh,n_nuc) + real(c_double), intent(in) :: r1(n_grid1,3), r2(n_grid2,3) + real(c_double), intent(in) :: rn(n_nuc,3) + real(c_double), intent(in) :: wr1(n_grid1), wr2(n_grid2) + real(c_double), intent(in) :: aos_data1(n_grid1,ao_num,4), aos_data2(n_grid2,ao_num,4) + real(c_double), intent(out) :: int2_grad1_u12(n_grid1,ao_num,ao_num,4) + real(c_double), intent(out) :: tc_int_2e_ao(ao_num,ao_num,ao_num,ao_num) + + end subroutine + + end interface + +end module + +! --- + diff --git a/plugins/local/tc_int/write_tc_int_gpu.irp.f b/plugins/local/tc_int/write_tc_int_gpu.irp.f new file mode 100644 index 00000000..c0dd9c90 --- /dev/null +++ b/plugins/local/tc_int/write_tc_int_gpu.irp.f @@ -0,0 +1,56 @@ +! --- + +program write_tc_int_gpu + + implicit none + + print *, ' j2e_type = ', j2e_type + print *, ' j1e_type = ', j1e_type + print *, ' env_type = ', env_type + + my_grid_becke = .True. + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + my_extra_grid_becke = .True. + PROVIDE tc_grid2_a tc_grid2_r + my_n_pt_r_extra_grid = tc_grid2_r + my_n_pt_a_extra_grid = tc_grid2_a + touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid + + call write_int(6, my_n_pt_r_grid, 'radial external grid over') + call write_int(6, my_n_pt_a_grid, 'angular external grid over') + + call write_int(6, my_n_pt_r_extra_grid, 'radial internal grid over') + call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over') + + call main() + +end + +! --- + +subroutine main() + + implicit none + + PROVIDE io_tc_integ + + print*, 'io_tc_integ = ', io_tc_integ + + if(io_tc_integ .ne. "Write") then + print*, 'io_tc_integ != Write' + print*, io_tc_integ + stop + endif + + call provide_int2_grad1_u12_ao_gpu() + + call ezfio_set_tc_keywords_io_tc_integ('Read') + +end + +! --- + diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f index 16414f39..e9c003d4 100644 --- a/src/dft_utils_in_r/ao_in_r.irp.f +++ b/src/dft_utils_in_r/ao_in_r.irp.f @@ -52,35 +52,39 @@ END_PROVIDER BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_grid,3)] - BEGIN_DOC - ! aos_grad_in_r_array(i,j,k) = value of the kth component of the gradient of ith ao on the jth grid point - ! - ! k = 1 : x, k= 2, y, k 3, z - END_DOC + BEGIN_DOC + ! + ! aos_grad_in_r_array(i,j,k) = value of the kth component of the gradient of ith ao on the jth grid point + ! + ! k = 1 : x, k= 2, y, k 3, z + ! + END_DOC - implicit none - integer :: i,j,m - double precision :: aos_array(ao_num), r(3) - double precision :: aos_grad_array(3,ao_num) - !$OMP PARALLEL DO & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,r,aos_array,aos_grad_array,j,m) & - !$OMP SHARED(aos_grad_in_r_array,n_points_final_grid,ao_num,final_grid_points) - do i = 1, n_points_final_grid - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array) - do m = 1, 3 - do j = 1, ao_num - aos_grad_in_r_array(j,i,m) = aos_grad_array(m,j) - enddo + implicit none + integer :: i, j, m + double precision :: aos_array(ao_num), r(3) + double precision :: aos_grad_array(3,ao_num) + + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,m,r,aos_array,aos_grad_array) & + !$OMP SHARED(aos_grad_in_r_array,n_points_final_grid,ao_num,final_grid_points) + do i = 1, n_points_final_grid + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array) + do m = 1, 3 + do j = 1, ao_num + aos_grad_in_r_array(j,i,m) = aos_grad_array(m,j) + enddo + enddo enddo - enddo - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO +END_PROVIDER - END_PROVIDER +! --- BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp, (3,ao_num,n_points_final_grid)] @@ -205,18 +209,53 @@ BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_gri END_PROVIDER - BEGIN_PROVIDER[double precision, aos_in_r_array_extra_transp, (n_points_extra_final_grid,ao_num)] - implicit none - BEGIN_DOC - ! aos_in_r_array_extra_transp(i,j) = value of the jth ao on the ith grid point - END_DOC - integer :: i,j - double precision :: aos_array(ao_num), r(3) - do i = 1, n_points_extra_final_grid - do j = 1, ao_num - aos_in_r_array_extra_transp(i,j) = aos_in_r_array_extra(j,i) +! --- + +BEGIN_PROVIDER[double precision, aos_in_r_array_extra_transp, (n_points_extra_final_grid,ao_num)] + + BEGIN_DOC + ! aos_in_r_array_extra_transp(i,j) = value of the jth ao on the ith grid point + END_DOC + + implicit none + integer :: i, j + double precision :: aos_array(ao_num), r(3) + + do i = 1, n_points_extra_final_grid + do j = 1, ao_num + aos_in_r_array_extra_transp(i,j) = aos_in_r_array_extra(j,i) + enddo enddo - enddo - END_PROVIDER +END_PROVIDER + +! --- + +BEGIN_PROVIDER[double precision, aos_grad_in_r_array_extra, (ao_num,n_points_extra_final_grid,3)] + + implicit none + integer :: i, j, m + double precision :: aos_array(ao_num), r(3) + double precision :: aos_grad_array(3,ao_num) + + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,m,r,aos_array,aos_grad_array) & + !$OMP SHARED(aos_grad_in_r_array_extra,n_points_extra_final_grid,ao_num,final_grid_points_extra) + do i = 1, n_points_extra_final_grid + r(1) = final_grid_points_extra(1,i) + r(2) = final_grid_points_extra(2,i) + r(3) = final_grid_points_extra(3,i) + call give_all_aos_and_grad_at_r(r, aos_array, aos_grad_array) + do m = 1, 3 + do j = 1, ao_num + aos_grad_in_r_array_extra(j,i,m) = aos_grad_array(m,j) + enddo + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + +! ---