Optimize kernels in_perm_2d_k and out_perm_2d_k.

include/cuda_func.hpp

@@ -14,9 +14,10 @@
 #ifdef __cplusplus
 extern "C" {
 #endif
-  void test_d(uintptr_t, uintptr_t, uintptr_t, uintptr_t, int, cudaStream_t*);
+  void texture_bind_d (char*, size_t);
+  void texture_unbind_d ();
   void in_perm_d (char*, size_t*, char*, char*, size_t, size_t, size_t, cudaStream_t*);
-  void out_perm_d (double*, char*, size_t*, char*, char*, size_t, size_t, size_t, cudaStream_t*);
+  void out_perm_d (double*, char*, size_t*, char*, char*, size_t, size_t, size_t, cudaStream_t*, size_t*, size_t*, size_t);
 #ifdef __cplusplus
 }
 #endif
@@ -24,13 +25,21 @@ extern "C" {
 template <class Real_t>
 class cuda_func {
   public:
+    static void texture_bind_h (char *ptr_d, size_t len);
+    static void texture_unbind_h ();
     static void in_perm_h (char *precomp_data, char *input_perm, char *input_data, char *buff_in, 
         size_t interac_indx, size_t M_dim0, size_t vec_cnt);
     static void out_perm_h (char *scaling, char *precomp_data, char *output_perm, char *output_data, char *buff_out, 
-        size_t interac_indx, size_t M_dim0, size_t vec_cnt);
+        size_t interac_indx, size_t M_dim0, size_t vec_cnt, size_t *tmp_a, size_t *tmp_b, size_t counter);
     static void compare_h (Real_t *gold, Real_t *mine, size_t n);
 };
 
+template <class Real_t>
+void cuda_func<Real_t>::texture_bind_h (char *ptr_d, size_t len) { texture_bind_d(ptr_d, len); };
+
+template <class Real_t>
+void cuda_func<Real_t>::texture_unbind_h () { texture_unbind_d(); };
+
 template <class Real_t>
 void cuda_func<Real_t>::in_perm_h (
   char *precomp_data,
@@ -56,13 +65,23 @@ void cuda_func<Real_t>::out_perm_h (
   char *buff_out,
   size_t interac_indx,
   size_t M_dim1,
-  size_t vec_cnt )
+  size_t vec_cnt,
+  size_t *tmp_a,
+  size_t *tmp_b,
+  size_t counter )
 {
   cudaStream_t *stream;
   stream = pvfmm::CUDA_Lock::acquire_stream(0);
+  size_t *a_d, *b_d;
+  cudaMalloc((void**)&a_d, sizeof(size_t)*counter);
+  cudaMalloc((void**)&b_d, sizeof(size_t)*counter);
+  cudaMemcpy(a_d, tmp_a, sizeof(size_t)*counter, cudaMemcpyHostToDevice);
+  cudaMemcpy(b_d, tmp_b, sizeof(size_t)*counter, cudaMemcpyHostToDevice);
   out_perm_d((double *) scaling, precomp_data, (size_t *) output_perm, output_data, buff_out, 
-	  interac_indx, M_dim1, vec_cnt, stream);
-}
+	  interac_indx, M_dim1, vec_cnt, stream, a_d, b_d, counter);
+  cudaFree(a_d);
+  cudaFree(b_d);
+};
 
 template <class Real_t>
 void cuda_func<Real_t>::compare_h (
@@ -84,6 +103,6 @@ void cuda_func<Real_t>::compare_h (
 	  } 
   }
   free(mine_h);
-}
+};
 
 #endif //_CUDA_FUNC_HPP_

include/fmm_pts.txx

@@ -1574,6 +1574,8 @@ void FMM_Pts<FMMNode>::EvalList_cuda(SetupData<Real_t>& setup_data) {
 
       //std::cout << "Before Loop. " << '\n';
 
+      //cuda_func<Real_t>::texture_bind_h (input_perm_d, input_perm.Dim());
+
       for (size_t k = 0; k < interac_blk.Dim(); k++) {
         size_t vec_cnt = 0;
 
@@ -1582,6 +1584,7 @@ void FMM_Pts<FMMNode>::EvalList_cuda(SetupData<Real_t>& setup_data) {
           vec_cnt += interac_cnt[j];
 
         /* CPU Version */
+        /*
         char* buff_in =&buff[0];
         char* buff_out=&buff[vec_cnt*dof*M_dim0*sizeof(Real_t)];
 
@@ -1597,9 +1600,10 @@ void FMM_Pts<FMMNode>::EvalList_cuda(SetupData<Real_t>& setup_data) {
             << input_perm[(interac_indx + i + 1)*4 + 2] 
             << '\n';
           for(size_t j = 0; j < M_dim0; j++) {
-            //v_out[j] = v_in[perm[j]]*scal[j];
+            v_out[j] = v_in[perm[j]]*scal[j];
           }
         }
+        */
 
         /* GPU Kernel call */	
         char *buff_in_d = (char *) buff_d.dev_ptr;
@@ -1621,7 +1625,7 @@ void FMM_Pts<FMMNode>::EvalList_cuda(SetupData<Real_t>& setup_data) {
             vec_cnt1 += interac_cnt[j];
 
           /* Compare buff_in */
-          std::cout << M_dim0*vec_cnt0*dof*sizeof(Real_t) << '\n';
+          //std::cout << M_dim0*vec_cnt0*dof*sizeof(Real_t) << '\n';
           /*
           cuda_func<Real_t>::compare_h(
               (Real_t *) (buff_in   + M_dim0*vec_cnt0*dof*sizeof(Real_t)), 
@@ -1649,7 +1653,7 @@ void FMM_Pts<FMMNode>::EvalList_cuda(SetupData<Real_t>& setup_data) {
               (Real_t *) (buff_out_d + M_dim1*vec_cnt0*dof*sizeof(Real_t)), 
               dof*vec_cnt1*M_dim1);
 */
-          std::cout << '\n';
+          //std::cout << '\n';
 
           vec_cnt0 += vec_cnt1;
         }
@@ -1659,6 +1663,34 @@ void FMM_Pts<FMMNode>::EvalList_cuda(SetupData<Real_t>& setup_data) {
         if (error != cudaSuccess) 
           std::cout << "cublasXgemm(): " << cudaGetErrorString(error) << '\n';
 
+        size_t *tmp_a, *tmp_b;
+        size_t counter = 0;
+        size_t last = -1;
+        if (vec_cnt > 0) {
+          int i;
+          cudaMallocHost((void**)&tmp_a, sizeof(size_t)*vec_cnt);
+          cudaMallocHost((void**)&tmp_b, sizeof(size_t)*vec_cnt);
+          tmp_a[0] = 0;
+          for (i = 1; i < vec_cnt; i++) {
+            if (output_perm[(interac_indx + i)*4 + 3] != last) {
+              last = output_perm[(interac_indx + i)*4 + 3];
+              tmp_b[counter] = i;
+              counter ++;
+              tmp_a[counter] = i;
+            }
+          }
+          tmp_b[counter] = i;
+          for (i = 0; i < 12; i++) std::cout << tmp_a[i] << ", ";
+          std::cout << '\n';
+          for (i = 0; i < 12; i++) std::cout << tmp_b[i] << ", ";
+          std::cout << '\n';
+          for (i = 0; i < counter; i++) {
+            if (tmp_a[i] == tmp_b[i]) std::cout << tmp_a[i] << ", " << tmp_b[i] << '\n';
+          }
+          std::cout << counter << '\n';
+        }
+
+        /*
         for(int i = 0; i < vec_cnt; i++) {
           Real_t scaling_factor=scaling[interac_indx+i];
           const PERM_INT_T*  perm=(PERM_INT_T*)(precomp_data[0]+output_perm[(interac_indx+i)*4+0]);
@@ -1666,11 +1698,18 @@ void FMM_Pts<FMMNode>::EvalList_cuda(SetupData<Real_t>& setup_data) {
           const     Real_t* v_in =(    Real_t*)(    buff_out   +output_perm[(interac_indx+i)*4+2]);
           Real_t*           v_out=(    Real_t*)( output_data[0]+output_perm[(interac_indx+i)*4+3]);
           for(size_t j=0;j<M_dim1;j++ ){
-            //v_out[j]+=v_in[perm[j]]*scal[j]*scaling_factor;
+            v_out[j]+=v_in[perm[j]]*scal[j]*scaling_factor;
           }
         }
+        */
+
         cuda_func<Real_t>::out_perm_h (scaling_d, (char *) precomp_data_d.dev_ptr, output_perm_d, 
-            (char *) output_data_d.dev_ptr, buff_out_d, interac_indx, M_dim1, vec_cnt);
+            (char *) output_data_d.dev_ptr, buff_out_d, interac_indx, M_dim1, vec_cnt, tmp_a, tmp_b, counter);
+        if (vec_cnt > 0) {
+          cudaFreeHost(tmp_a);
+          cudaFreeHost(tmp_b);
+        }
+
         //CUDA_Lock::wait(0); 
         error = cudaGetLastError();
         if (error != cudaSuccess) 
@@ -1679,6 +1718,8 @@ void FMM_Pts<FMMNode>::EvalList_cuda(SetupData<Real_t>& setup_data) {
         interac_indx += vec_cnt;
         interac_blk_dsp += interac_blk[k];
       }
+
+      //cuda_func<Real_t>::texture_unbind_h ();
     }
   }
   dummy.Device2Host();

m4/ax_check_cuda.m4

@@ -79,7 +79,7 @@ then
    then
       NVCCFLAGS+=" -g"
    else
-      NVCCFLAGS+=" "
+      NVCCFLAGS+=" -O3"
    fi
    AC_ARG_ENABLE([emu],
       AC_HELP_STRING([--enable-emu],[turn on device emulation for CUDA]),

src/cuda_func.cu

@@ -4,6 +4,32 @@
 
 #define DEFAULT_NUM_THREAD 128
 
+texture<int, 1, cudaReadModeElementType> tex_in_perm;
+texture<double, 1, cudaReadModeElementType> tex_out_perm;
+
+
+__global__ void in_perm_2d_k (
+  char *precomp_data,
+  size_t *input_perm,
+  char *input_data,
+  char *buff_in,
+  size_t interac_indx,
+  size_t M_dim0,
+  size_t vec_cnt )
+{
+  int tidx = blockIdx.x*blockDim.x + threadIdx.x;
+  int tidy = blockIdx.y*blockDim.y + threadIdx.y;
+
+  if (tidy< vec_cnt && tidx < M_dim0) {
+    size_t k = (interac_indx + tidy)*4;
+    size_t *perm  = (size_t *) (precomp_data + input_perm[k + 0]);
+    double *scal  = (double *) (precomp_data + input_perm[k + 1]);
+    double *v_in  = (double *) (input_data ) + input_perm[k + 3];
+    double *v_out = (double *) (buff_in      + input_perm[k + 2]);
+    v_out[tidx] = v_in[perm[tidx]]*scal[tidx];
+  }
+}
+
 /* Case: double */
 __global__ void in_perm_k (
   char *precomp_data,
@@ -19,11 +45,19 @@ __global__ void in_perm_k (
 
   if (tid < vec_cnt) {
     /* Convert to ptr. */
+    
     size_t *perm  = (size_t *) (precomp_data + input_perm[(interac_indx + tid)*4 + 0]);
     double *scal  = (double *) (precomp_data + input_perm[(interac_indx + tid)*4 + 1]);
     double *v_in  = (double *) (input_data ) + input_perm[(interac_indx + tid)*4 + 3];
     double *v_out = (double *) (buff_in      + input_perm[(interac_indx + tid)*4 + 2]);
 
+    /*
+    size_t *perm  = (size_t *) (precomp_data + tex1Dfetch(tex_in_perm, (interac_indx + tid)*4 + 0));
+    double *scal  = (double *) (precomp_data + tex1Dfetch(tex_in_perm, (interac_indx + tid)*4 + 1));
+    double *v_in  = (double *) (input_data ) + tex1Dfetch(tex_in_perm, (interac_indx + tid)*4 + 3);
+    double *v_out = (double *) (buff_in      + tex1Dfetch(tex_in_perm, (interac_indx + tid)*4 + 2));
+    */
+
     /* PRAM Model: assuming as many threads as we need. */
     for (size_t j = 0; j < M_dim0; j++) {
       v_out[j] = v_in[perm[j]]*scal[j];
@@ -39,7 +73,73 @@ __global__ void in_perm_k (
     */
   }
   //__syncthreads();
-}
+};
+
+__global__ void out_perm_2d_k (
+  double *scaling,
+  char *precomp_data,
+  size_t *output_perm,
+  char *output_data,
+  char *buff_out,
+  size_t interac_indx,
+  size_t M_dim1,
+  size_t vec_cnt,
+  size_t *a_d,
+  size_t *b_d,
+  size_t counter )
+{
+  /* 2-dim thread Id. */
+  int tidx = blockIdx.x*blockDim.x + threadIdx.x;
+  int tidy = blockIdx.y*blockDim.y + threadIdx.y;
+
+  int a = a_d[tidy];
+  int b = b_d[tidy];
+  /*
+  int a = (tidy*vec_cnt)/vec_cnt;
+  int b = ((tidy + 1)*vec_cnt)/vec_cnt;
+
+  if (tidy > 0 && a < vec_cnt) {
+    size_t out_ptr = output_perm[(interac_indx + a)*4 + 3];
+    if (tidy > 0) while (a < vec_cnt && out_ptr == output_perm[(interac_indx + a)*4 + 3]) a++;
+  }
+  if (tidy < vec_cnt - 1 &&  - 1 && b < vec_cnt) {
+    size_t out_ptr = output_perm[(interac_indx + b)*4 + 3];
+    if (tidy < vec_cnt - 1) while (b < vec_cnt && out_ptr == output_perm[(interac_indx + b)*4 + 3]) b++;
+  }
+*/
+
+  //if (tidy < vec_cnt && tidx < M_dim1) {
+  if (tidy < counter && tidx < M_dim1) {
+    double v = 0.0;
+    
+    for(int i = a; i < b; i++) { // Compute permutations.
+      double scaling_factor = scaling[interac_indx + i];
+      size_t k = (interac_indx + i)*4;
+      size_t *perm = (size_t*) (precomp_data + output_perm[k + 0]);
+      double *scal = (double*) (precomp_data + output_perm[k + 1]);
+      double *v_in = (double*) (buff_out     + output_perm[k + 2]);
+      //double *v_out = (double*) (output_data)+ output_perm[k + 3];
+      //v_out[tidx] += v_in[perm[tidx]]*scal[tidx]*scaling_factor;
+      v += v_in[perm[tidx]]*scal[tidx]*scaling_factor;
+    }
+    double *v_out = (double*) (output_data)+ output_perm[(interac_indx + a)*4 + 3];
+    v_out[tidx] += v;
+
+    
+    /*
+    double *v_out = (double*) (output_data)+ output_perm[(interac_indx + a)*4 + 3];
+    for(int i = a; i < b; i++) { // Compute permutations.
+      double scaling_factor = scaling[interac_indx + i];
+      size_t *perm = (size_t*) (precomp_data + output_perm[(interac_indx + i)*4 + 0]);
+      double *scal = (double*) (precomp_data + output_perm[(interac_indx + i)*4 + 1]);
+      double *v_in = (double*) (buff_out     + output_perm[(interac_indx + i)*4 + 2]);
+      double *v_out = (double*) (output_data)+ output_perm[(interac_indx + a)*4 + 3];
+      v += v_in[perm[tidx]]*scal[tidx]*scaling_factor;
+    }
+    if (a < b) v_out[tidx] = v;
+    */
+  }
+};
 
 __global__ void out_perm_k (
   double *scaling,
@@ -49,25 +149,33 @@ __global__ void out_perm_k (
   char *buff_out,
   size_t interac_indx,
   size_t M_dim1,
-  size_t vec_cnt )
+  size_t vec_cnt,
+  size_t *a_d,
+  size_t *b_d,
+  size_t counter )
 {
   /* 1-dim thread Id. */
   int tid = blockIdx.x*blockDim.x + threadIdx.x;
 
   /* Specifing range. */
+  int a = a_d[tid];
+  int b = b_d[tid];
+  /*
   int a = (tid*vec_cnt)/vec_cnt;
   int b = ((tid + 1)*vec_cnt)/vec_cnt;
 
-  if (tid > 0 && a < vec_cnt) { // Find 'a' independent of other threads.
+  if (tid > 0 && a < vec_cnt) {
     size_t out_ptr = output_perm[(interac_indx + a)*4 + 3];
     if (tid > 0) while (a < vec_cnt && out_ptr == output_perm[(interac_indx + a)*4 + 3]) a++;
   }
-  if (tid < vec_cnt - 1 &&  - 1 && b < vec_cnt) { // Find 'b' independent of other threads.
+  if (tid < vec_cnt - 1 &&  - 1 && b < vec_cnt) {
     size_t out_ptr = output_perm[(interac_indx + b)*4 + 3];
     if (tid < vec_cnt - 1) while (b < vec_cnt && out_ptr == output_perm[(interac_indx + b)*4 + 3]) b++;
   }
+  */
 
-  if (tid < vec_cnt) {
+  //if (tid < vec_cnt) {
+  if (tid < counter) {
     /* PRAM Model: assuming as many threads as we need. */
     for(int i = a; i < b; i++) { // Compute permutations.
       double scaling_factor = scaling[interac_indx + i];
@@ -79,11 +187,16 @@ __global__ void out_perm_k (
     }
   }
   //__syncthreads();
-}
+};
 
 extern "C" { 
-void test_d (uintptr_t precomp_data, uintptr_t input_perm, uintptr_t input_data, uintptr_t buff_in, 
-  int interac_indx, cudaStream_t *stream) {};
+void texture_bind_d (char *input_perm, size_t len) {
+  cudaBindTexture(0, tex_in_perm, input_perm, len);
+};
+
+void texture_unbind_d () {
+  cudaUnbindTexture(tex_in_perm);
+};
 
 void in_perm_d (
     char *precomp_data,
@@ -95,18 +208,37 @@ void in_perm_d (
     size_t vec_cnt, 
     cudaStream_t *stream )
 {
+  if (vec_cnt == 0) return;
+
+  float time_ms;
   int n_thread, n_block;
   n_thread = DEFAULT_NUM_THREAD;
   n_block = vec_cnt/n_thread + 1;
+
+  cudaEvent_t beg, end;
+  cudaEventCreate(&beg);
+  cudaEventCreate(&end);
+  //cudaBindTexture(0, tex_in_perm, input_perm, );
   /*
-     in_perm_k<<<n_thread, n_block, 0, *stream>>>(precomp_data, input_perm, input_data, buff_in, 
-     interac_indx, M_dim0, vec_cnt);
-     */
-  if (vec_cnt == 0) return;
   printf("in_perm_k : vec_cnt: %d, M_dim0: %d, n_block: %d, n_thread: %d\n", 
       (int) vec_cnt, (int) M_dim0, n_block, n_thread);
+  */
+  cudaEventRecord(beg, 0);
+  /*
   in_perm_k<<<n_block, n_thread>>>(precomp_data, input_perm, input_data, buff_in, 
       interac_indx, M_dim0, vec_cnt);
+      */
+  dim3 dimBlock(16, 32);
+  dim3 dimGrid(M_dim0/16 + 1, vec_cnt/32 + 1);
+  in_perm_2d_k<<<dimGrid, dimBlock>>>(precomp_data, input_perm, input_data, buff_in, 
+      interac_indx, M_dim0, vec_cnt);
+  cudaEventRecord(end, 0);
+  cudaEventSynchronize(end);
+  cudaEventElapsedTime(&time_ms, beg, end);
+  printf("in_perm_d : %f ms\n", time_ms);	
+    
+  cudaEventDestroy(beg);
+  cudaEventDestroy(end);
 };
 
 void out_perm_d (
@@ -118,20 +250,44 @@ void out_perm_d (
   size_t interac_indx,
   size_t M_dim1,
   size_t vec_cnt,
-  cudaStream_t *stream )
+  cudaStream_t *stream,
+  size_t *a_d,
+  size_t *b_d,
+  size_t counter )
 {
+  if (vec_cnt == 0) return;
+
+  float time_ms;
   int n_thread, n_block;
   n_thread = DEFAULT_NUM_THREAD;
-  n_block = vec_cnt/n_thread + 1;
-/*
-  out_perm_k<<<n_block, n_thread, 0, *stream>>>(scaling, precomp_data, output_perm, output_data, buff_out, 
-    interac_indx, M_dim0, vec_cnt);
-*/
-  if (vec_cnt == 0) return;
+  //n_block = vec_cnt/n_thread + 1;
+  n_block = counter/n_thread + 1;
+
+  cudaEvent_t beg, end;
+  cudaEventCreate(&beg);
+  cudaEventCreate(&end);
+  /*
   printf("out_perm_k : vec_cnt: %d, M_dim0: %d, n_block: %d, n_thread: %d\n", 
       (int) vec_cnt, (int) M_dim1, n_block, n_thread);
+      */
+  cudaEventRecord(beg, 0);
+  /*
   out_perm_k<<<n_block, n_thread>>>(scaling, precomp_data, output_perm, output_data, buff_out, 
-    interac_indx, M_dim1, vec_cnt);
+    interac_indx, M_dim1, vec_cnt, a_d, b_d, counter);
+    */
+  
+  dim3 dimBlock(16, 32);
+  dim3 dimGrid(M_dim1/8 + 1, vec_cnt/64 + 1);
+  out_perm_2d_k<<<dimGrid, dimBlock>>>(scaling, precomp_data, output_perm, output_data, buff_out, 
+    interac_indx, M_dim1, vec_cnt, a_d, b_d, counter);
+    
+  cudaEventRecord(end, 0);
+  cudaEventSynchronize(end);
+  cudaEventElapsedTime(&time_ms, beg, end);
+  printf("out_perm_d : %f ms\n", time_ms);	
+    
+  cudaEventDestroy(beg);
+  cudaEventDestroy(end);
 };
 
 }