  *  @param[in,out] g_ybar All \f$P_{m,n}(\bar{Y})\f$ fixed to \f$P_{m,n}(Y)\f$
  *  @param[out] g_ysum All \f$s(P_{m,n}(Y))\f$
  */
-__global__
-void algSAT_stage2(const int c_n_size, const int c_m_size, const int c_width, float *g_ybar,
-    float *g_ysum )
-{
-  const int tx = threadIdx.x, ty = threadIdx.y,
-        bx = blockIdx.x, col0 = bx*MW+ty, col = col0*WS+tx;
-
-  if( col >= c_width ) return;
+__global__ void algSAT_stage2(const int c_n_size, const int c_m_size,
+                              const int c_width, float *g_ybar, float *g_ysum)
+{
+  const int tx = threadIdx.x, ty = threadIdx.y, bx = blockIdx.x,
+            col0 = bx * MW + ty, col = col0 * WS + tx;
+
+  if (col >= c_width)
+    return;
 
   g_ybar += col;
   float y = *g_ybar;
-  int ln = HWS+tx;
-
-  if( tx == WS-1 )
+  int ln = HWS + tx;
+
+  if (tx == WS - 1)
     g_ysum += col0;
 
-  volatile __shared__ float s_block[ MW ][ HWS+WS+1 ];
-
-  if( tx < HWS ) s_block[ty][tx] = 0.f;
-  else s_block[ty][ln] = 0.f;
-
-  for (int n = 1; n < c_n_size; ++n) {
+  volatile __shared__ float s_block[MW][HWS + WS + 1];
+
+  if (tx < HWS)
+    s_block[ty][tx] = 0.f;
+  else
+    s_block[ty][ln] = 0.f;
+
+  for (int n = 1; n < c_n_size; ++n)
+  {
 
     // calculate ysum -----------------------
 
     s_block[ty][ln] = y;
 
-    s_block[ty][ln] += s_block[ty][ln-1];
-    s_block[ty][ln] += s_block[ty][ln-2];
-    s_block[ty][ln] += s_block[ty][ln-4];
-    s_block[ty][ln] += s_block[ty][ln-8];
-    s_block[ty][ln] += s_block[ty][ln-16];
-
-    if( tx == WS-1 ) {
+    s_block[ty][ln] += s_block[ty][ln - 1];
+    s_block[ty][ln] += s_block[ty][ln - 2];
+    s_block[ty][ln] += s_block[ty][ln - 4];
+    s_block[ty][ln] += s_block[ty][ln - 8];
+    s_block[ty][ln] += s_block[ty][ln - 16];
+
+    if (tx == WS - 1)
+    {
       *g_ysum = s_block[ty][ln];
       g_ysum += c_m_size;
     }


texture<float, 2, cudaReadModeElementType> textureReadDepth;

inline unsigned int divUp(unsigned int a, unsigned int b) {
  return (a % b != 0) ? (a / b + 1) : (a / b);
}

#define WS 32 // Warp size (defines b x b block size where b = WS)
#define HWS 16 // Half Warp Size
#define DW 8 // Default number of warps (computational block height)
#define CHW 7 // Carry-heavy number of warps (computational block height for some kernels)
#define OW 6 // Optimized number of warps (computational block height for some kernels)
#define DNB 6 // Default number of blocks per SM (minimum blocks per SM launch bounds)
#define ONB 5 // Optimized number of blocks per SM (minimum blocks per SM for some kernels)
#define MTS 192 // Maximum number of threads per block with 8 blocks per SM
#define MBO 8 // Maximum number of blocks per SM using optimize or maximum warps
#define CHB 7 // Carry-heavy number of blocks per SM using default number of warps
#define MW 6 // Maximum number of warps per block with 8 blocks per SM (with all warps computing)
#define SOW 5 // Dual-scheduler optimized number of warps per block (with 8 blocks per SM and to use the dual scheduler with 1 computing warp)
#define MBH 3 // Maximum number of blocks per SM using half-warp size

/**
 *  @ingroup gpu

 *  @param[in] g_y All \f$P_{m,n}(Y)\f$
 *  @param[in] g_v All \f$P^T_{m,n}(V)\f$
 */
__global__ void algSAT_stage4(const int c_width, const int c_height,
                              float *g_out, const float *g_in, const float *g_y,
                              const float *g_v)
{
  const int tx = threadIdx.x, ty = threadIdx.y, bx = blockIdx.x,
            by = blockIdx.y, col = bx * WS + tx, row0 = by * WS;

  __shared__ float s_block[WS][WS + 1];

  float(*bdata)[WS + 1] = (float(*)[WS + 1]) & s_block[ty][tx];

  g_in += (row0 + ty) * c_width + col;
  if (by > 0)
    g_y += (by - 1) * c_width + col;
  if (bx > 0)
    g_v += (bx - 1) * c_height + row0 + tx;

#pragma unroll
  for (int i = 0; i < WS - (WS % SOW); i += SOW)
  {
    **bdata = *g_in;
    bdata += SOW;
    g_in += SOW * c_width;
  }
  if (ty < WS % SOW)
  {
    **bdata = *g_in;
  }

  __syncthreads();

  if (ty == 0)
  {

    {  // calculate y -----------------------
      float(*bdata)[WS + 1] = (float(*)[WS + 1]) & s_block[0][tx];

      float prev;
      if (by > 0)
        prev = *g_y;
      else
        prev = 0.f;

#pragma unroll
      for (int i = 0; i < WS; ++i, ++bdata)
        **bdata = prev = **bdata + prev;
    }

    {  // calculate x -----------------------
      float *bdata = s_block[tx];

      float prev;
      if (bx > 0)
        prev = *g_v;
      else
        prev = 0.f;

#pragma unroll
      for (int i = 0; i < WS; ++i, ++bdata)
        *bdata = prev = *bdata + prev;
    }
  }

  __syncthreads();

  bdata = (float(*)[WS + 1]) & s_block[ty][tx];

  g_out += (row0 + ty) * c_width + col;

#pragma unroll
  for (int i = 0; i < WS - (WS % SOW); i += SOW)
  {
    *g_out = **bdata;
    bdata += SOW;
    g_out += SOW * c_width;
  }
  if (ty < WS % SOW)
  {
    *g_out = **bdata;
  }
}
