title: CUDA Kernel Optimization regarding register
body: <p>I'm quite new to CUDA and GPU programming. I'm trying to write a Kernel for an application in physics. The parallelization is made over a quadrature of directions, each direction resulting in a sweep of a 2D cartesian domain. Here is the kernel. it actually works well, giving good results.</p>
  
  <p>However, a very high number of registers per blocks leads to a spill to local memory that harshly slow down the code performance.</p>
  
  <pre><code>__global__ void KERNEL (int imax, int jmax, int mmax, int lg, int lgmax, 
                          double *x,   double *y, double *qd, double *kappa,
                          double *S, double *G, double *qw, double *SkG, 
                          double *Ska,double *a, double *Ljm, int *data)
  
  {
  int m   = 1+blockIdx.x*blockDim.x + threadIdx.x ; 
  int tid = threadIdx.x ;
  
  //Var needed for thread execution
  ...
  
  extern __shared__ double shared[] ;
  
  //Read some data from Global mem
  mu  = qd[        (m-1)];
  eta = qd[  MSIZE+(m-1)];
  wm  = qd[3*MSIZE+(m-1)];
  amu = fabs(mu);
  aeta= fabs(eta);
  ista = data[        (m-1)] ;
  iend = data[1*MSIZE+(m-1)] ;
  istp = data[2*MSIZE+(m-1)] ;
  jsta = data[3*MSIZE+(m-1)] ;
  jend = data[4*MSIZE+(m-1)] ;
  jstp = data[5*MSIZE+(m-1)] ;
  
  j1 = (1-jstp)   ;
  j2 = (1+jstp)/2 ;
  i1 = (1-istp)   ;
  i2 = (1+istp)/2 ;
  
  isw = ista-istp ;
  jsw = jsta-jstp ;
  
  dy = dx = 1.0e-2 ;
  
  for(i=1 ; i&lt;=imax; i++) Ljm[MSIZE*(i-1)+m] = S[jsw*(imax+2)+i] ;
  
  //Beginning of the vertical Sweep, can be from left to right, 
  // or opposite depending on the thread
  
  for(j=jsta ; j1*jend + j2*j&lt;=j2*jend + j1*j ; j=j+jstp) {
  
  Lw = S[j*(imax+2)+isw] ;
  
  //Beginning of the horizontal Sweep, can be from left to right, 
  // or opposite depending on the thread
  
     for(i=ista ; i1*iend + i2*i&lt;=i2*iend + i1*i ; i=i+istp) {
  
              ax = dy ;
              Lx = ax*amu/ex ;
              ay = dx ;
              Ly = ay*aeta/ey ;
  
              dv = ax*ay ;
              L0 = dv*kappaij ;
              Sp = S[j*(imax+2)+i]*dv ;
              Ls = Ljm[MSIZE*(i-1)+m] ;
  
              Lp = (Lx*Lw+Ly*Ls+Sp)/(Lx+Ly+L0) ;
  
              Lw = Lw+(Lp-Lw)/ex ;
              Ls = Ls+(Lp-Ls)/ey ;
  
              Ljm[MSIZE*(i-1)+m] = Ls ;
  
              shared[tid] = wm*Lp ;
              __syncthreads();
  
              for (s=16; s&gt;0; s&gt;&gt;=1) {
                  if (tid &lt; s) {
                     shared[tid] += shared[tid + s] ;
                  }
              }
  
              if(tid==0) atomicAdd(&amp;SkG[imax*(j-1)+(i-1)],shared[tid]*kappaij);
  
      }
      // End of horizontal sweep
   }
   // End of vertical sweep
  </code></pre>
  
  <p>}</p>
  
  <p>How can i optimize the execution of this code ? I run it over 8 blocks of 32 threads.
  The occupancy for this kernel is really low, limited by the registers according to the Visual profiler.</p>
  
  <p>I have no idea on how to improve it.</p>
  
  <p>Thanks !</p>