 
     g_ybar += c_width;
     y = *g_ybar += y;
-
   }
 }
 

 *  @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_inout,
    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_inout += (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_inout;
    bdata += SOW;
    g_inout += SOW * c_width;
  }
  if( ty < WS%SOW ) {
    **bdata = *g_inout;
  }

  __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_inout -= (WS-(WS%SOW))*c_width;

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

 *  @param[out] g_ybar All \f$P_{m,n}(\bar{Y})\f$
 *  @param[out] g_vhat All \f$P^T_{m,n}(\hat{V})\f$
 */
__global__ void algSAT_stage1(const int c_width, const int c_height,
                              const float *g_in, float *g_ybar, float *g_vhat)
{

  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;
  g_ybar += by * c_width + col;
  g_vhat += bx * 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 ybar -----------------------
      float(*bdata)[WS + 1] = (float(*)[WS + 1]) & s_block[0][tx];

      float prev = **bdata;
      ++bdata;
