code: def _get_fprop_lrn(clss, compute_capability):
      """
      Local Response Normalization (LRN) layer.
      Implementation based on fprop_avg kernel.
  
      Implements the following operation:
      for each output pixel
          x' = x / response
      where the response is
          response = [1 + alpha/N * sum_neighbors x_neighbor**2 ]**beta
      so we compute the pooling output
      """
      code = r"""
  %(common)s
  
  __global__ void spool_fprop_lrn(
      const %(type)s* I, %(type)s* O, %(type)s* A,
      float alpha, float beta, float ascale, float bpower, int flags,
      int N, int W, int H, int D, int C,
      int WN, int HWN, int DHWN, int P, int Q,
      int magic_P, int shift_P, int QN, int PQN, int MPQN,
      int pad_c, int pad_d, int pad_h, int pad_w,
      int str_c, int str_d, int str_h, int str_w,
      int S, int RS, int RST, int JRST,
      int magic_S, int shift_S,
      int magic_RS, int shift_RS, int magic_RST, int shift_RST,
      int supP, int supQ, int shlP, int maskP, int shrP,
      int shlQ, int maskQ, int shrQ, int maskN, int shrN
      %(stats_args)s
      )
  {
      __shared__ float rcpWindowSize;
      extern __shared__ int lut[];
  
      int tid = threadIdx.x;
  
      // paralellism is over QMPK dimensions (output pixels and ofm's)
      int n  = tid;
      int q  = blockIdx.x;
      int mp = blockIdx.y;
      int k  = blockIdx.z;
  
      int m = mp * magic_P; m >>= shift_P;
      int p = mp - m*P;
  
      // zigzag q back and forth to improve L2 cache perf
      if (p & 1)
          q = Q - q - 1;
  
      const %(type)s* IonO = I;  // input pixel at output location
      I += n;
      IonO += k*MPQN + m*PQN + p*QN + q*N + n;
      O += k*MPQN + m*PQN + p*QN + q*N + n;
      A += k*MPQN + m*PQN + p*QN + q*N + n;
  
      float O_val = beta != 0.0f ? %(cvt)s(__ldg(O)) : 0.0f;
  
      if (tid < 32)
      {
          int kj = k * str_c - pad_c;
          int mt = m * str_d - pad_d;
          int pr = p * str_h - pad_h;
          int qs = q * str_w - pad_w;
  
          int window_size = 0;
          int jrst = tid;
          // this loop generates the LUT (same for pooling and normalization)
          while (jrst < JRST)
          {
              int j = jrst * magic_RST; j >>= shift_RST;
              int rst = jrst - j * RST;
  
              int t = rst * magic_RS; t >>= shift_RS;
              int rs = rst - t * RS;
  
              int r = rs * magic_S; r >>= shift_S;
              int s = rs - r*S;
  
              int x = qs + s;
              int y = pr + r;
              int z = mt + t;
              int c = kj + j;
  
              bool bounds_x  = x >= 0 && x < W;
              bool bounds_y  = y >= 0 && y < H;
              bool bounds_z  = z >= 0 && z < D;
              bool bounds_c  = c >= 0 && c < C;
              bool in_bounds = bounds_x && bounds_y && bounds_z && bounds_c;
  
              // Count the total valid slices
              window_size += __popc(__ballot(in_bounds));
  
              int sliceI  = c*DHWN + z*HWN + y*WN + x*N;
  
              lut[jrst] = in_bounds ? sliceI : -1;
              jrst += 32;
          }
  
          if(tid == 0)
          {
              //rcpWindowSize = 1.0f / (float)window_size;
              rcpWindowSize = (float)RST/(float)JRST;
          }
      }
      __syncthreads();
  
      float out = 0.0f;
      float denom;
      float sumsquare = 0.0f;
      float input = 0.0f;
      int jrst = 0;
      while (jrst < JRST)
      {
          int slice0 = lut[jrst + 0];
          int slice1 = lut[jrst + 1];
          int slice2 = lut[jrst + 2];
          int slice3 = lut[jrst + 3];
  
          // TODO: May not need to load all slices if they are not used.
          input =      jrst + 0 < JRST && slice0 >= 0 ? %(cvt)s(__ldg(I + slice0)) : 0.0f;
          sumsquare += jrst + 0 < JRST && slice0 >= 0 ? input * input: 0.0f;
          input =      jrst + 1 < JRST && slice1 >= 0 ? %(cvt)s(__ldg(I + slice1)) : 0.0f;
          sumsquare += jrst + 1 < JRST && slice1 >= 0 ? input * input: 0.0f;
          input =      jrst + 2 < JRST && slice2 >= 0 ? %(cvt)s(__ldg(I + slice2)) : 0.0f;
          sumsquare += jrst + 2 < JRST && slice2 >= 0 ? input * input: 0.0f;
          input =      jrst + 3 < JRST && slice3 >= 0 ? %(cvt)s(__ldg(I + slice3)) : 0.0f;
          sumsquare += jrst + 3 < JRST && slice3 >= 0 ? input * input: 0.0f;
  
          jrst += 4;
      }
  
      denom = (1 + ascale*sumsquare*rcpWindowSize);
      out = %(cvt)s(__ldg(IonO)) / powf(denom, bpower);
  
  
      // convert back to fp to write out
      %(type)s temp_out = %(cvt_out)s( %(mul_by_scale)s (out*alpha + O_val*beta));
  
      // predicate write with no-op flag
      if (!(flags & 1)) {
          *O = temp_out;
          *A = %(cvt_out)s( %(mul_by_scale)s denom );  // write the denomiantor to address
      }
  
      // collect max abs stats
      int intermediate_max = max_abs(0, temp_out); // compute abs
      %(atomic_max)s
  }
  """
  
      template_vals = prepare_template_vals(clss, compute_capability)
      code = code % template_vals
      module = SourceModule(code)
      kernel = module.get_function("spool_fprop_lrn")
      kernel.prepare(
          "3P 4f 34I 10I" + ("Pf" if (clss[0] == "x") else "")
      )  # add superblocking parameter
      return kernel
  
file_path: neon/neon/backends/kernels/cuda/pooling.py
imports: from builtins import str
  from pycuda.tools import context_dependent_memoize
  from neon.backends import cuda_templates
  from neon.backends.cuda_templates import _common_fp16_to_fp32, _common_round, _common_max_abs, _common_kepler, _ew_types
  from neon.backends.util.source_module import SourceModule
label_desc: Use image node_id 5 for calling a global function with example usage: _get_fprop_lrn(clss, compute_capability) and returns: kernel