From 66566a8ce790d5206760b2c1fe5d3b457c616c7d Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Tue, 26 Nov 2024 23:00:15 +0100 Subject: [PATCH] added Aph_dRPA_sing kernel --- src/RPA/phRRPA.f90 | 3 +- src/cuda/src/ph_drpa.c | 16 +++++------ src/cuda/src/phlr_drpa_a_sing.cu | 49 ++++++++++++++++++-------------- src/mod/cu_quack_module.f90 | 10 +++---- 4 files changed, 42 insertions(+), 36 deletions(-) diff --git a/src/RPA/phRRPA.f90 b/src/RPA/phRRPA.f90 index 94c6576..7d00061 100644 --- a/src/RPA/phRRPA.f90 +++ b/src/RPA/phRRPA.f90 @@ -80,9 +80,10 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO, ! DEBUG allocate(Om_gpu(nS), XpY_gpu(nS,nS), XmY_gpu(nS,nS)) - call ph_drpa(nO, nBas, eHF(1), ERI(1,1,1,1), Om_gpu(1), XpY_gpu(1,1), XmY_gpu(1,1)) + call ph_drpa(nO, nBas, nS, eHF(1), ERI(1,1,1,1), Om_gpu(1), XpY_gpu(1,1), XmY_gpu(1,1)) print *, ' CPU:', Aph(1,1) print *, ' GPU:', XpY_gpu(1,1) + print *, ' GPU:', XmY_gpu(1,1) stop call phLR(TDA,nS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY) diff --git a/src/cuda/src/ph_drpa.c b/src/cuda/src/ph_drpa.c index 61afcbe..0e0f671 100644 --- a/src/cuda/src/ph_drpa.c +++ b/src/cuda/src/ph_drpa.c @@ -7,7 +7,7 @@ #include "ph_drpa.h" -int ph_drpa(int nO, int nBas, double *h_eps, double *h_ERI, +int ph_drpa(int nO, int nBas, int nS, double *h_eps, double *h_ERI, double *h_Omega, double *h_XpY, double *h_XmY) { double *d_eps; @@ -16,25 +16,23 @@ int ph_drpa(int nO, int nBas, double *h_eps, double *h_ERI, int nBas2 = nBas * nBas; int nBas4 = nBas2 * nBas2; - int ia, jb; - int nS = nO * (nBas - nO); for (ia = 0; ia < nS; ia++) { h_Omega[ia] = 0.0; for (jb = 0; jb < nS; jb++) { - h_XmY[jb + nO * nBas * ia] = 0.0; - h_XpY[jb + nO * nBas * ia] = 0.0; + h_XmY[jb + ia * nS] = 4.0; + //h_XpY[jb + ia * nS] = 5.0; } } - check_Cuda_Errors(cudaMalloc((void**)&d_eps, nO * sizeof(double)), + check_Cuda_Errors(cudaMalloc((void**)&d_eps, nBas * sizeof(double)), "cudaMalloc", __FILE__, __LINE__); check_Cuda_Errors(cudaMalloc((void**)&d_ERI, nBas4 * sizeof(double)), "cudaMalloc", __FILE__, __LINE__); - check_Cuda_Errors(cudaMemcpy(d_eps, h_eps, nO * sizeof(double), cudaMemcpyHostToDevice), + check_Cuda_Errors(cudaMemcpy(d_eps, h_eps, nBas * sizeof(double), cudaMemcpyHostToDevice), "cudaMemcpy", __FILE__, __LINE__); check_Cuda_Errors(cudaMemcpy(d_ERI, h_ERI, nBas4 * sizeof(double), cudaMemcpyHostToDevice), "cudaMemcpy", __FILE__, __LINE__); @@ -42,8 +40,8 @@ int ph_drpa(int nO, int nBas, double *h_eps, double *h_ERI, // construct A matrix double *d_A; check_Cuda_Errors(cudaMalloc((void**)&d_A, nS * nS * sizeof(double)), "cudaMalloc", __FILE__, __LINE__); -// phLR_dRPA_A_sing(nO, nBas, d_eps, d_ERI, d_A); -// check_Cuda_Errors(cudaGetLastError(), "cudaGetLastError", __FILE__, __LINE__); + phLR_dRPA_A_sing(nO, nBas, d_eps, d_ERI, d_A); + check_Cuda_Errors(cudaGetLastError(), "cudaGetLastError", __FILE__, __LINE__); check_Cuda_Errors(cudaMemcpy(h_XpY, d_A, nS * nS * sizeof(double), cudaMemcpyDeviceToHost), diff --git a/src/cuda/src/phlr_drpa_a_sing.cu b/src/cuda/src/phlr_drpa_a_sing.cu index 3cbc556..7844d66 100644 --- a/src/cuda/src/phlr_drpa_a_sing.cu +++ b/src/cuda/src/phlr_drpa_a_sing.cu @@ -4,39 +4,46 @@ __global__ void phLR_dRPA_A_sing_kernel(int nO, int nBas, double *eps, double *E int i, j, a, b; - int ia, jb, jb_off; + int aa, bb; + int nV, nS, nVS; + int nBas2, nBas3; + int i_A0, i_A1, i_A2; + int i_I0, i_I1, i_I2; - int ij_off0, ij_off; + nV = nBas - nO; + nS = nO * nV; + nVS = nV * nS; - int aa_max = nBas - nO; - int ia_max = aa_max * nO; + nBas2 = nBas * nBas; + nBas3 = nBas2 * nBas; - int nBas2 = nBas * nBas; - int nBas3 = nBas2 * nBas; + aa = blockIdx.x * blockDim.x + threadIdx.x; + bb = blockIdx.y * blockDim.y + threadIdx.y; - int aa = blockIdx.x * blockDim.x + threadIdx.x; - int bb = blockIdx.y * blockDim.y + threadIdx.y; - - while(aa < aa_max) { + while(aa < nV) { a = aa + nO; - ij_off0 = a * nBas2; + i_A0 = aa * nS; + i_I0 = a * nBas2; - while(bb < aa_max) { + while(bb < nV) { b = bb + nO; - ij_off = ij_off0 + b * nBas; + i_A1 = i_A0 + bb; + i_I1 = i_I0 + b * nBas; + i = 0; while(i < nO) { - ia = i * aa_max + aa; - jb_off = ia * ia_max; - - while(j < nO) { - jb = j * aa_max + bb; - A[jb + jb_off] = 2.0 * ERI[i + j * nBas3 + ij_off]; - if(a==b && i==j) { - A[jb + jb_off] += eps[a] - eps[i]; + i_A2 = i_A1 + i * nVS; + i_I2 = i_I1 + i; + + j = 0; + while(j < nO) { + + A[i_A2 + j * nV] = 2.0 * ERI[i_I2 + j * nBas3]; + if((a==b) && (i==j)) { + A[i_A2 + j * nV] += eps[a] - eps[i]; } j ++; diff --git a/src/mod/cu_quack_module.f90 b/src/mod/cu_quack_module.f90 index f73cc6c..a786587 100644 --- a/src/mod/cu_quack_module.f90 +++ b/src/mod/cu_quack_module.f90 @@ -8,16 +8,16 @@ module cu_quack_module interface - subroutine ph_drpa(nO, nBas, eps, ERI, & + subroutine ph_drpa(nO, nBas, nS, eps, ERI, & Omega, XpY, XmY) bind(C, name = "ph_drpa") import c_int, c_double - integer(c_int), intent(in), value :: nO, nBas + integer(c_int), intent(in), value :: nO, nBas, nS real(c_double), intent(in) :: eps(nBas) real(c_double), intent(in) :: ERI(nBas,nBas,nBas,nBas) - real(c_double), intent(out) :: Omega(nO*nBas) - real(c_double), intent(out) :: XpY(nO*nBas,nO*nBas) - real(c_double), intent(out) :: XmY(nO*nBas,nO*nBas) + real(c_double), intent(out) :: Omega(nS) + real(c_double), intent(out) :: XpY(nS,nS) + real(c_double), intent(out) :: XmY(nS,nS) end subroutine ph_drpa