int2_grad1_u12_ao is computed correctly on CUDA

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
+

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
+  PROVIDE nxBlocks nyBlocks nzBlocks
-  blockSize = 32
+  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,                                                                        &
+  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,                &
+                             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, &
+                             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

plugins/local/tc_int/gpu_module.F90

@@ -38,61 +38,23 @@ module gpu_module
 
     ! ---
 
-    subroutine deb_int_long_range(nBlocks, blockSize,            &
+    subroutine deb_int2_grad1_u12_ao(nxBlocks, nyBlocks, nzBlocks, blockxSize, blockySize, blockzSize, &
-                                  n_grid2, n_ao, wr2, aos_data2, &
+                                     n_grid1, n_grid2, n_ao, n_nuc, size_bh,             &
-                                  int_fct_long_range) bind(C, name = "deb_int_long_range")
+                                     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
-      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")
 
       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)        :: r1(3,n_grid1)
-      real(c_double), intent(in)        :: r2(n_grid2,3)
+      real(c_double), intent(in)        :: r2(3,n_grid2)
-      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)        :: 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)