From f43ee8cf61e6a59a1e6d8b013cd5fe75b1724b42 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Thu, 1 Aug 2024 18:29:50 +0200 Subject: [PATCH] int2_grad1_u12_ao is computed correctly on CUDA --- plugins/local/tc_int/EZFIO.cfg | 36 +++ plugins/local/tc_int/deb_tc_int_cuda.irp.f | 283 +++------------------ plugins/local/tc_int/gpu_module.F90 | 58 +---- 3 files changed, 75 insertions(+), 302 deletions(-) create mode 100644 plugins/local/tc_int/EZFIO.cfg diff --git a/plugins/local/tc_int/EZFIO.cfg b/plugins/local/tc_int/EZFIO.cfg new file mode 100644 index 00000000..12366f01 --- /dev/null +++ b/plugins/local/tc_int/EZFIO.cfg @@ -0,0 +1,36 @@ +[nxBlocks] +type: integer +doc: nb of x blocks in the Grid +interface: ezfio,provider,ocaml +default: 10 + +[nyBlocks] +type: integer +doc: nb of y blocks in the Grid +interface: ezfio,provider,ocaml +default: 10 + +[nzBlocks] +type: integer +doc: nb of z blocks in the Grid +interface: ezfio,provider,ocaml +default: 1 + +[blockxSize] +type: integer +doc: size of x blocks +interface: ezfio,provider,ocaml +default: 32 + +[blockySize] +type: integer +doc: size of y blocks +interface: ezfio,provider,ocaml +default: 32 + +[blockzSize] +type: integer +doc: size of z blocks +interface: ezfio,provider,ocaml +default: 1 + diff --git a/plugins/local/tc_int/deb_tc_int_cuda.irp.f b/plugins/local/tc_int/deb_tc_int_cuda.irp.f index 6eb3a8b5..a7b5b125 100644 --- a/plugins/local/tc_int/deb_tc_int_cuda.irp.f +++ b/plugins/local/tc_int/deb_tc_int_cuda.irp.f @@ -36,8 +36,6 @@ subroutine main() implicit none - !call deb_int_long_range_gpu() - !call deb_int_bh_kernel_gpu() call deb_int2_grad1_u12_ao_gpu() return @@ -45,246 +43,6 @@ end ! --- -subroutine deb_int_long_range_gpu() - - use gpu_module - - implicit none - - integer :: i, j, k - integer :: ipoint, jpoint - - integer :: nBlocks, blockSize - - double precision :: acc_thr, err_tot, nrm_tot, err_loc - - double precision :: time0, time1 - double precision :: cuda_time0, cuda_time1 - double precision :: cpu_time0, cpu_time1 - - double precision, allocatable :: aos_data2(:,:,:) - double precision, allocatable :: int_fct_long_range(:,:,:) - double precision, allocatable :: int_fct_long_range_gpu(:,:,:) - - - - call wall_time(time0) - print*, ' start deb_int_long_range_gpu' - - - ! --- - - nBlocks = 256 - blockSize = 32 - - allocate(aos_data2(n_points_extra_final_grid,ao_num,4)) - allocate(int_fct_long_range_gpu(n_points_extra_final_grid,ao_num,ao_num)) - - do k = 1, ao_num - do ipoint = 1, n_points_extra_final_grid - aos_data2(ipoint,k,1) = aos_in_r_array_extra(k,ipoint) - aos_data2(ipoint,k,2) = aos_grad_in_r_array_extra(k,ipoint,1) - aos_data2(ipoint,k,3) = aos_grad_in_r_array_extra(k,ipoint,2) - aos_data2(ipoint,k,4) = aos_grad_in_r_array_extra(k,ipoint,3) - enddo - enddo - - ! --- - - call wall_time(cuda_time0) - - call deb_int_long_range(nBlocks, blockSize, & - n_points_extra_final_grid, ao_num, final_weight_at_r_vector_extra, aos_data2, & - int_fct_long_range_gpu) - - call wall_time(cuda_time1) - print*, ' wall time for CUDA kernel (min) = ', (cuda_time1-cuda_time0) / 60.d0 - - deallocate(aos_data2) - - ! --- - - allocate(int_fct_long_range(n_points_extra_final_grid,ao_num,ao_num)) - - call wall_time(cpu_time0) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (j, i, jpoint) & - !$OMP SHARED (int_fct_long_range, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp) - !$OMP DO SCHEDULE (static) - do j = 1, ao_num - do i = 1, ao_num - do jpoint = 1, n_points_extra_final_grid - int_fct_long_range(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - call wall_time(cpu_time1) - print*, ' wall time on CPU (min) = ', (cpu_time1-cpu_time0) / 60.d0 - - ! --- - - acc_thr = 1d-12 - err_tot = 0.d0 - nrm_tot = 0.d0 - - do j = 1, ao_num - do i = 1, ao_num - do jpoint = 1, n_points_extra_final_grid - err_loc = dabs(int_fct_long_range(jpoint,i,j) - int_fct_long_range_gpu(jpoint,i,j)) - if(err_loc > acc_thr) then - print*, " error on", jpoint, i, j - print*, " CPU res", int_fct_long_range (jpoint,i,j) - print*, " GPU res", int_fct_long_range_gpu(jpoint,i,j) - stop - endif - err_tot = err_tot + err_loc - nrm_tot = nrm_tot + dabs(int_fct_long_range(jpoint,i,j)) - enddo - enddo - enddo - - print *, ' absolute accuracy (%) =', 100.d0 * err_tot / nrm_tot - - ! --- - - deallocate(int_fct_long_range) - deallocate(int_fct_long_range_gpu) - - call wall_time(time1) - print*, ' wall time for deb_int_long_range_gpu (min) = ', (time1-time0) / 60.d0 - - return -end - -! --- - -subroutine deb_int_bh_kernel_gpu() - - use gpu_module - - implicit none - - integer :: m - integer :: ipoint, jpoint - - integer :: nBlocks, blockSize - - double precision :: acc_thr, err_tot, nrm_tot, err_loc - - double precision :: time0, time1 - double precision :: cuda_time0, cuda_time1 - double precision :: cpu_time0, cpu_time1 - - double precision, allocatable :: r1(:,:), r2(:,:) - double precision, allocatable :: grad1_u12(:,:,:) - double precision, allocatable :: grad1_u12_gpu(:,:,:) - - - - call wall_time(time0) - print*, ' start deb_int_bh_kernel_gpu' - - - ! --- - - allocate(r1(n_points_final_grid,3)) - allocate(r2(n_points_extra_final_grid,3)) - - do ipoint = 1, n_points_final_grid - r1(ipoint,1) = final_grid_points(1,ipoint) - r1(ipoint,2) = final_grid_points(2,ipoint) - r1(ipoint,3) = final_grid_points(3,ipoint) - enddo - - do ipoint = 1, n_points_extra_final_grid - r2(ipoint,1) = final_grid_points_extra(1,ipoint) - r2(ipoint,2) = final_grid_points_extra(2,ipoint) - r2(ipoint,3) = final_grid_points_extra(3,ipoint) - enddo - - ! --- - - nBlocks = 256 - blockSize = 32 - - allocate(grad1_u12_gpu(n_points_extra_final_grid,n_points_final_grid,4)) - - call wall_time(cuda_time0) - - call deb_int_bh_kernel(nBlocks, blockSize, & - n_points_final_grid, n_points_extra_final_grid, ao_num, nucl_num, jBH_size, & - r1, r2, nucl_coord, jBH_c, jBH_m, jBH_n, jBH_o, & - grad1_u12_gpu) - - call wall_time(cuda_time1) - print*, ' wall time for CUDA kernel (min) = ', (cuda_time1-cuda_time0) / 60.d0 - - ! --- - - allocate(grad1_u12(n_points_extra_final_grid,n_points_final_grid,4)) - - call wall_time(cpu_time0) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint) & - !$OMP SHARED (n_points_final_grid, n_points_extra_final_grid, grad1_u12) - !$OMP DO - do ipoint = 1, n_points_final_grid - call get_grad1_u12_for_tc(ipoint, n_points_extra_final_grid, grad1_u12(1,ipoint,1) & - , grad1_u12(1,ipoint,2) & - , grad1_u12(1,ipoint,3) & - , grad1_u12(1,ipoint,4) ) - enddo - !$OMP END DO - !$OMP END PARALLEL - - call wall_time(cpu_time1) - print*, ' wall time on CPU (min) = ', (cpu_time1-cpu_time0) / 60.d0 - - ! --- - - acc_thr = 1d-12 - err_tot = 0.d0 - nrm_tot = 0.d0 - - do m = 1, 4 - do ipoint = 1, n_points_final_grid - do jpoint = 1, n_points_extra_final_grid - err_loc = dabs(grad1_u12(jpoint,ipoint,m) - grad1_u12_gpu(jpoint,ipoint,m)) - if(err_loc > acc_thr) then - print*, " error on", jpoint, ipoint, m - print*, " CPU res", grad1_u12 (jpoint,ipoint,m) - print*, " GPU res", grad1_u12_gpu(jpoint,ipoint,m) - stop - endif - err_tot = err_tot + err_loc - nrm_tot = nrm_tot + dabs(grad1_u12(jpoint,ipoint,m)) - enddo - enddo - enddo - - print *, ' absolute accuracy (%) =', 100.d0 * err_tot / nrm_tot - - ! --- - - deallocate(r1, r2) - deallocate(grad1_u12) - deallocate(grad1_u12_gpu) - - call wall_time(time1) - print*, ' wall time for deb_int_bh_kernel_gpu (min) = ', (time1-time0) / 60.d0 - - return -end - -! --- - subroutine deb_int2_grad1_u12_ao_gpu() use gpu_module @@ -295,15 +53,14 @@ subroutine deb_int2_grad1_u12_ao_gpu() integer :: i, j, k integer :: ipoint, jpoint - integer :: nBlocks, blockSize - double precision :: acc_thr, err_tot, nrm_tot, err_loc double precision :: time0, time1 double precision :: cuda_time0, cuda_time1 double precision :: cpu_time0, cpu_time1 + double precision :: cpu_ttime0, cpu_ttime1 - double precision, allocatable :: r1(:,:), r2(:,:), aos_data2(:,:,:) + double precision, allocatable :: r1(:,:), r2(:,:), rn(:,:), aos_data2(:,:,:) double precision, allocatable :: grad1_u12(:,:,:), int_fct_long_range(:,:,:) double precision, allocatable :: int2_grad1_u12_ao(:,:,:,:) double precision, allocatable :: int2_grad1_u12_ao_gpu(:,:,:,:) @@ -318,6 +75,7 @@ subroutine deb_int2_grad1_u12_ao_gpu() allocate(r1(n_points_final_grid,3)) allocate(r2(n_points_extra_final_grid,3)) + allocate(rn(3,nucl_num)) allocate(aos_data2(n_points_extra_final_grid,ao_num,4)) do ipoint = 1, n_points_final_grid @@ -332,6 +90,12 @@ subroutine deb_int2_grad1_u12_ao_gpu() r2(ipoint,3) = final_grid_points_extra(3,ipoint) enddo + do k = 1, nucl_num + rn(1,k) = nucl_coord(k,1) + rn(2,k) = nucl_coord(k,2) + rn(3,k) = nucl_coord(k,3) + enddo + do k = 1, ao_num do ipoint = 1, n_points_extra_final_grid aos_data2(ipoint,k,1) = aos_in_r_array_extra(k,ipoint) @@ -343,20 +107,21 @@ subroutine deb_int2_grad1_u12_ao_gpu() ! --- - nBlocks = 256 - blockSize = 32 + PROVIDE nxBlocks nyBlocks nzBlocks + PROVIDE blockxSize blockySize blockzSize allocate(int2_grad1_u12_ao_gpu(ao_num,ao_num,n_points_final_grid,4)) call wall_time(cuda_time0) - call deb_int2_grad1_u12_ao(nBlocks, blockSize, & - n_points_final_grid, n_points_extra_final_grid, ao_num, nucl_num, jBH_size, & - r1, r2, final_weight_at_r_vector_extra, nucl_coord, aos_data2, jBH_c, jBH_m, jBH_n, jBH_o, & + call deb_int2_grad1_u12_ao(nxBlocks, nyBlocks, nzBlocks, blockxSize, blockySize, blockzSize, & + n_points_final_grid, n_points_extra_final_grid, ao_num, nucl_num, jBH_size, & + final_grid_points, final_grid_points_extra, final_weight_at_r_vector_extra, rn, & + aos_data2, jBH_c, jBH_m, jBH_n, jBH_o, & int2_grad1_u12_ao_gpu) call wall_time(cuda_time1) - print*, ' wall time for CUDA kernel (min) = ', (cuda_time1-cuda_time0) / 60.d0 + write(*,"(A,2X,F15.7)") ' wall time for CUDA kernel (sec) = ', (cuda_time1 - cuda_time0) ! --- @@ -367,6 +132,7 @@ subroutine deb_int2_grad1_u12_ao_gpu() call wall_time(cpu_time0) + call wall_time(cpu_ttime0) !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (j, i, jpoint) & @@ -381,8 +147,11 @@ subroutine deb_int2_grad1_u12_ao_gpu() enddo !$OMP END DO !$OMP END PARALLEL + call wall_time(cpu_ttime1) + write(*,"(A,2X,F15.7)") ' wall time for int_long_range (sec) = ', (cpu_ttime1 - cpu_ttime0) + call wall_time(cpu_ttime0) !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (ipoint) & @@ -396,15 +165,21 @@ subroutine deb_int2_grad1_u12_ao_gpu() enddo !$OMP END DO !$OMP END PARALLEL + call wall_time(cpu_ttime1) + write(*,"(A,2X,F15.7)") ' wall time for tc_int_bh (sec) = ', (cpu_ttime1 - cpu_ttime0) + + call wall_time(cpu_ttime0) do m = 1, 4 call dgemm("T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, 1.d0 & , int_fct_long_range(1,1,1), n_points_extra_final_grid, grad1_u12(1,1,m), n_points_extra_final_grid & , 0.d0, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num) enddo + call wall_time(cpu_ttime1) + write(*,"(A,2X,F15.7)") ' wall time for DGEMM (sec) = ', (cpu_ttime1 - cpu_ttime0) call wall_time(cpu_time1) - print*, ' wall time on CPU (min) = ', (cpu_time1-cpu_time0) / 60.d0 + write(*,"(A,2X,F15.7)") ' wall time on cpu (sec) = ', (cpu_time1 - cpu_time0) ! --- @@ -434,13 +209,13 @@ subroutine deb_int2_grad1_u12_ao_gpu() ! --- - deallocate(r1, r2, aos_data2) + deallocate(r1, r2, rn, aos_data2) deallocate(int_fct_long_range, grad1_u12) deallocate(int2_grad1_u12_ao) deallocate(int2_grad1_u12_ao_gpu) call wall_time(time1) - print*, ' wall time for deb_int2_grad1_u12_ao_gpu (min) = ', (time1-time0) / 60.d0 + write(*,"(A,2X,F15.7)") ' wall time for deb_int2_grad1_u12_ao_gpu (sec) = ', (time1 - time0) return end diff --git a/plugins/local/tc_int/gpu_module.F90 b/plugins/local/tc_int/gpu_module.F90 index f625e57d..1f07a2ea 100644 --- a/plugins/local/tc_int/gpu_module.F90 +++ b/plugins/local/tc_int/gpu_module.F90 @@ -38,61 +38,23 @@ module gpu_module ! --- - subroutine deb_int_long_range(nBlocks, blockSize, & - n_grid2, n_ao, wr2, aos_data2, & - int_fct_long_range) bind(C, name = "deb_int_long_range") - - import c_int, c_double - integer(c_int), intent(in), value :: nBlocks, blockSize - integer(c_int), intent(in), value :: n_grid2 - integer(c_int), intent(in), value :: n_ao - real(c_double), intent(in) :: wr2(n_grid2) - real(c_double), intent(in) :: aos_data2(n_grid2,n_ao,4) - real(c_double), intent(out) :: int_fct_long_range(n_grid2,n_ao,n_ao) - - end subroutine deb_int_long_range - - ! --- - - subroutine deb_int_bh_kernel(nBlocks, blockSize, & - n_grid1, n_grid2, n_ao, n_nuc, size_bh, & - r1, r2, rn, c_bh, m_bh, n_bh, o_bh, & - grad1_u12) bind(C, name = "deb_int_bh_kernel") + subroutine deb_int2_grad1_u12_ao(nxBlocks, nyBlocks, nzBlocks, blockxSize, blockySize, blockzSize, & + n_grid1, n_grid2, n_ao, n_nuc, size_bh, & + r1, r2, wr2, rn, aos_data2, c_bh, m_bh, n_bh, o_bh, & + int2_grad1_u12_ao) bind(C, name ="deb_int2_grad1_u12_ao") import c_int, c_double, c_ptr - integer(c_int), intent(in), value :: nBlocks, blockSize + integer(c_int), intent(in), value :: nxBlocks, blockxSize + integer(c_int), intent(in), value :: nyBlocks, blockySize + integer(c_int), intent(in), value :: nzBlocks, blockzSize integer(c_int), intent(in), value :: n_grid1, n_grid2 integer(c_int), intent(in), value :: n_ao integer(c_int), intent(in), value :: n_nuc integer(c_int), intent(in), value :: size_bh - real(c_double), intent(in) :: r1(n_grid1,3) - real(c_double), intent(in) :: r2(n_grid2,3) - real(c_double), intent(in) :: rn(n_nuc,3) - real(c_double), intent(in) :: c_bh(size_bh,n_nuc) - 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(out) :: grad1_u12(n_grid2,n_grid1,4) - - end subroutine deb_int_bh_kernel - - ! --- - - subroutine deb_int2_grad1_u12_ao(nBlocks, blockSize, & - n_grid1, n_grid2, n_ao, n_nuc, size_bh, & - r1, r2, wr2, rn, aos_data2, c_bh, m_bh, n_bh, o_bh, & - int2_grad1_u12_ao) bind(C, name ="deb_int2_grad1_u12_ao") - - import c_int, c_double, c_ptr - integer(c_int), intent(in), value :: nBlocks, blockSize - integer(c_int), intent(in), value :: n_grid1, n_grid2 - integer(c_int), intent(in), value :: n_ao - integer(c_int), intent(in), value :: n_nuc - integer(c_int), intent(in), value :: size_bh - real(c_double), intent(in) :: r1(n_grid1,3) - real(c_double), intent(in) :: r2(n_grid2,3) + real(c_double), intent(in) :: r1(3,n_grid1) + real(c_double), intent(in) :: r2(3,n_grid2) real(c_double), intent(in) :: wr2(n_grid2) - real(c_double), intent(in) :: rn(n_nuc,3) + real(c_double), intent(in) :: rn(3,n_nuc) real(c_double), intent(in) :: aos_data2(n_grid2,n_ao,4) real(c_double), intent(in) :: c_bh(size_bh,n_nuc) integer(c_int), intent(in) :: m_bh(size_bh,n_nuc)