Slow kernel if I don't specify bounds at compile time

Accelerated Computing CUDA CUDA Programming and Performance
,
1

Hello,

I am trying to port my CFD solver on GPU exploiting both CUDA-fortran and openMP.

I started “isolating” the single subroutines (and creating simple models) just to understand how to map thread to DoFs. In the past days I have isolated a subroutine which serves for computing the numeric fluxes in my solver and I have noticed that it is really slower compared to the other subroutines (interpolation and other stuffs).

So, I created a specific module (not in the solver), but a “dummy module” which serves for testing. I have implemented a “simplified version”, basically a version of the subroutine that accesses only to the left state for computing the numeric fluxes and it is really slow. This is the subroutine:

!>@brief The subroutine used for computing the interface fluxes in cuda.
attributes(global) SUBROUTINE dummy_interfaceFluxes_CUDA_v0(nEQ,ninterf,order,nc,qfi,qfj,qfk,f1,g2,h3)
USE cudafor
IMPLICIT NONE
!*
INTEGER,VALUE :: nEQ,ninterf,order,nc
REAL(KIND=RK),DEVICE :: qfi(neq,order+1,order,order,nc),qfj(neq,order,order+1,order,nc),&
qfk(neq,order,order,order+1,nc),&
F1(neq,order+1,order,order,nc),G2(neq,order,order+1,order,nc),H3(neq,order,order,order+1,nc)
!*
INTEGER :: ig
INTEGER :: ib,k,icl,icr,faml,famr,nfp,mfp,ifpl,jfpl,kfpl,k_face
REAL(KIND=RK) :: sign_l,Qfl(nEQ)

ig = threadIdx%x+(blockIdx%x-1)*blockDim%x
IF(ig.le.ninterf*order*order) THEN

!* DOFS on the faces:
k_face = tmap_fint(1,ig)
nfp = tmap_fint(2,ig)
mfp = tmap_fint(3,ig)
icl = if2c_d(k_face,1);icl=1
icr = if2c_d(k_face,2);icr=2
faml = if2fam_d(1,k_face);faml=1
famr = if2fam_d(2,k_face);famr=2

ifpl = Gfp2Lfp_d(5,mfp,nfp,k_face);ifpl=1
jfpl = Gfp2Lfp_d(6,mfp,nfp,k_face);jfpl=2
kfpl = Gfp2Lfp_d(7,mfp,nfp,k_face);kfpl=3
sign_l = 1
IF(faml.eq.1) THEN
if(ifpl.eq.1) sign_l = -1
Qfl(1:neq) = Qfi(1:neq,ifpl,jfpl,kfpl,icl)
else if(faml.eq.2) THEN
if(jfpl.eq.1) sign_l = -1
Qfl(1:nEQ) = Qfj(1:nEQ,ifpl,jfpl,kfpl,icl)
else if(faml.eq.3) THEN
if(kfpl.eq.1) sign_l = -1
Qfl(1:nEQ) = Qfk(1:nEQ,ifpl,jfpl,kfpl,icl)
ENDIF
ENDIF

END SUBROUTINE dummy_interfaceFluxes_CUDA_v0
!/~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~/

The time employed by this subroutine is about 5000ms (measured with cudaEvent).

Now, I tried to specify the array bounds at compile time trough a directive, resulting in this code:

!/~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~/
define ncells 166375
define neq 5
define NORDER 6

attributes(global) SUBROUTINE dummy_interfaceFluxes_CUDA_ctime(ninterf,order,qfi,qfj,qfk,f1,g2,h3)
USE cudafor
IMPLICIT NONE
!*
INTEGER,VALUE :: ninterf,order
REAL(KIND=RK),DEVICE :: qfi(nEQ,NORDER+1,NORDER,NORDER,ncells),qfj(nEQ,NORDER,NORDER+1,NORDER,ncells),&
qfk(nEQ,NORDER,NORDER,NORDER+1,ncells),&
F1(nEQ,NORDER+1,NORDER,NORDER,ncells),&
G2(nEQ,NORDER,NORDER+1,NORDER,ncells),&
H3(nEQ,NORDER,NORDER,NORDER+1,ncells)
!*
INTEGER :: ig
INTEGER :: ib,k,icl,icr,faml,famr,nfp,mfp,ifpl,jfpl,kfpl,k_face
REAL(KIND=RK) :: sign_l,Qfl(nEQ)

ig = threadIdx%x+(blockIdx%x-1)*blockDim%x
IF(ig.le.ninterf*order*order) THEN

!* DOFS on the faces:
k_face = tmap_fint(1,ig)
nfp = tmap_fint(2,ig)
mfp = tmap_fint(3,ig)
icl = if2c_d(k_face,1);icl=1
icr = if2c_d(k_face,2);icr=2
faml = if2fam_d(1,k_face);faml=1
famr = if2fam_d(2,k_face);famr=2

ifpl = Gfp2Lfp_d(5,mfp,nfp,k_face);ifpl=1
jfpl = Gfp2Lfp_d(6,mfp,nfp,k_face);jfpl=2
kfpl = Gfp2Lfp_d(7,mfp,nfp,k_face);kfpl=3
sign_l = 1
IF(faml.eq.1) THEN
if(ifpl.eq.1) sign_l = -1
Qfl(1:neq) = Qfi(1:neq,ifpl,jfpl,kfpl,icl)

!normf(1:3) = S1(1,1:3,ifpl,jfpl,kfpl,icl)

else if(faml.eq.2) THEN
if(jfpl.eq.1) sign_l = -1
Qfl(1:nEQ) = Qfj(1:nEQ,ifpl,jfpl,kfpl,icl)
!normf(1:3) = S2(2,1:3,ifpl,jfpl,kfpl,icl)
else if(faml.eq.3) THEN
if(kfpl.eq.1) sign_l = -1
Qfl(1:nEQ) = Qfk(1:nEQ,ifpl,jfpl,kfpl,icl)
!normf(1:3) = S3(3,1:3,ifpl,jfpl,kfpl,icl)
ENDIF
ENDIF

END SUBROUTINE dummy_interfaceFluxes_CUDA_ctime
!/~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~/

This version is really fast: it employs only 0.2 ms and I don’t know why.

I have also tried to flat the input arrays Qfi,Qfj,Qfk, but it does not work.

At this point, I don’t know if it is necessary to have array-bounds specified at compile time or I am missing something.

Thank you very much in advance,

Nicola.

2

Some suggestions, why sompile time array-bounds can speed up your kernel:

  • If your algorithms divide by the array size, it can be much slower than a precomputed multiplication with 1/array size
  • Some accesses may be put into registers instead of slower local memory (which actually is global memory, sometimes cached in L1). If loops access the arrays, the indices are dynamic, which is not compatible with register access