     k -= kb;
   }
 
-  if (k > 0) { saxpy(16, B[0], &a[ 0][ 0], &x[ 0]); B += ldb;
-  if (k > 1) { saxpy(15, B[0], &a[ 1][ 1], &x[ 1]); B += ldb;
-  if (k > 2) { saxpy(14, B[0], &a[ 2][ 2], &x[ 2]); B += ldb;
-  if (k > 3) { saxpy(13, B[0], &a[ 3][ 3], &x[ 3]); B += ldb;
-  if (k > 4) { saxpy(12, B[0], &a[ 4][ 4], &x[ 4]); B += ldb;
-  if (k > 5) { saxpy(11, B[0], &a[ 5][ 5], &x[ 5]); B += ldb;
-  if (k > 6) { saxpy(10, B[0], &a[ 6][ 6], &x[ 6]); B += ldb;
-  if (k > 7) { saxpy( 9, B[0], &a[ 7][ 7], &x[ 7]); B += ldb;
-  if (k > 8) { saxpy( 8, B[0], &a[ 8][ 8], &x[ 8]); B += ldb;
-  if (k > 9) { saxpy( 7, B[0], &a[ 9][ 9], &x[ 9]); B += ldb;
-  if (k >10) { saxpy( 6, B[0], &a[10][10], &x[10]); B += ldb;
-  if (k >11) { saxpy( 5, B[0], &a[11][11], &x[11]); B += ldb;
-  if (k >12) { saxpy( 4, B[0], &a[12][12], &x[12]); B += ldb;
-  if (k >13) { saxpy( 3, B[0], &a[13][13], &x[13]); B += ldb;
-  if (k >14) { saxpy( 2, B[0], &a[14][14], &x[14]); B += ldb;
-  if (k >15) { saxpy( 1, B[0], &a[15][15], &x[15]); }}}}}}}}}}}}}}}}
+  if (k > 0) { INNER_RIGHT_LOOP_DEC( 0); B += ldb;
+  if (k > 1) { INNER_RIGHT_LOOP_DEC( 1); B += ldb;
+  if (k > 2) { INNER_RIGHT_LOOP_DEC( 2); B += ldb;
+  if (k > 3) { INNER_RIGHT_LOOP_DEC( 3); B += ldb;
+  if (k > 4) { INNER_RIGHT_LOOP_DEC( 4); B += ldb;
+  if (k > 5) { INNER_RIGHT_LOOP_DEC( 5); B += ldb;
+  if (k > 6) { INNER_RIGHT_LOOP_DEC( 6); B += ldb;
+  if (k > 7) { INNER_RIGHT_LOOP_DEC( 7); B += ldb;
+  if (k > 8) { INNER_RIGHT_LOOP_DEC( 8); B += ldb;
+  if (k > 9) { INNER_RIGHT_LOOP_DEC( 9); B += ldb;
+  if (k >10) { INNER_RIGHT_LOOP_DEC(10); B += ldb;
+  if (k >11) { INNER_RIGHT_LOOP_DEC(11); B += ldb;
+  if (k >12) { INNER_RIGHT_LOOP_DEC(12); B += ldb;
+  if (k >13) { INNER_RIGHT_LOOP_DEC(13); B += ldb;
+  if (k >14) { INNER_RIGHT_LOOP_DEC(14); B += ldb;
+  if (k >15) { INNER_RIGHT_LOOP_DEC(15); }}}}}}}}}}}}}}}}
 
   if (m - bi - ti > 0)
     sscal(n - bj, alpha, x, X, ldx);

  // Process non-diagonal blocks as for SGEMM
  k = n - bj - nb;
  while (k > 0) {
    __syncthreads();

#pragma unroll
    for (int l = 0; l < kb; l += by)
      a[l + threadIdx.y][threadIdx.x] = A[l * lda];

    __syncthreads();

    if (k < kb) break;

#pragma unroll
    for (int l = 0; l < kb; l++) {
      saxpy(B[0], a[l], x);
      B += ldb;
    }

    __syncthreads();

    A += kb * lda;
    k -= kb;
  }

  for (int l = 0; l < k; l++) {
    saxpy(B[0], a[l], x);
    B += ldb;
  }

  if (m - bi - ti > 0)
    sscal(n - bj, alpha, x, X, ldx);
}

template <CBlasDiag diag,
          unsigned int mb, unsigned int nb, unsigned int kb,
          unsigned int bx, unsigned int by>
__global__ void strmmRLN(int m, int n,
                        float alpha, const float * __restrict__ A, int lda,
                        const float * __restrict__ B, int ldb,
                        float * __restrict__ X, int ldx) {

  const int bi = blockIdx.x * mb;       // Starting row of block of X
  const int bj = blockIdx.y * nb;       // Starting column of block of X
  int ti = threadIdx.y * bx + threadIdx.x;

  A += (bj + threadIdx.y) * lda + bj + threadIdx.x;
  B += bj * ldb + bi + ti;
  X += bj * ldx + bi + ti;

  __shared__ float a[kb][nb];

  float x[] = { 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
                0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f };

  // For Upper/Trans and Lower/NoTrans process diagonal first
  int k = min(n - bj, nb);
  while (k > 0) {
#pragma unroll
    for (int j = 0; j < nb; j += by)
      a[threadIdx.x][j + threadIdx.y] =
        (diag != CBlasNonUnit && threadIdx.x == j + threadIdx.y) ? 1.0f : A[j * lda];

    __syncthreads();

    if (k < kb) break;

    saxpy( 1, B[0], a[ 0], x); B += ldb;
    saxpy( 2, B[0], a[ 1], x); B += ldb;
    saxpy( 3, B[0], a[ 2], x); B += ldb;
    saxpy( 4, B[0], a[ 3], x); B += ldb;
    saxpy( 5, B[0], a[ 4], x); B += ldb;
    saxpy( 6, B[0], a[ 5], x); B += ldb;
    saxpy( 7, B[0], a[ 6], x); B += ldb;
    saxpy( 8, B[0], a[ 7], x); B += ldb;
    saxpy( 9, B[0], a[ 8], x); B += ldb;
    saxpy(10, B[0], a[ 9], x); B += ldb;
    saxpy(11, B[0], a[10], x); B += ldb;
    saxpy(12, B[0], a[11], x); B += ldb;
    saxpy(13, B[0], a[12], x); B += ldb;
    saxpy(14, B[0], a[13], x); B += ldb;
    saxpy(15, B[0], a[14], x); B += ldb;
    saxpy(16, B[0], a[15], x); B += ldb;

    __syncthreads();

    A += kb;
    k -= kb;
  }

  for (int ll = 0; ll < k; ll++) {
    saxpy(ll + 1, B[0], a[ll], x);
    B += ldb;
  }

  // Process non-diagonal blocks as for SGEMM
  k = n - bj - nb;
  while (k > 0) {

    __syncthreads();

#pragma unroll
    for (int j = 0; j < nb; j += by)
      a[threadIdx.x][j + threadIdx.y] = A[j * lda];

  // For Upper/NoTrans and Lower/Trans process diagonal last
  k = min(n - bj, nb);
  while (k > 0) {
#pragma unroll
    for (int j = 0; j < nb; j += by)
      a[threadIdx.x][j + threadIdx.y] = A[j * lda];

    __syncthreads();

    if (k < kb) break;

    INNER_RIGHT_LOOP_DEC( 0); B += ldb; INNER_RIGHT_LOOP_DEC( 1); B += ldb;
    INNER_RIGHT_LOOP_DEC( 2); B += ldb; INNER_RIGHT_LOOP_DEC( 3); B += ldb;
    INNER_RIGHT_LOOP_DEC( 4); B += ldb; INNER_RIGHT_LOOP_DEC( 5); B += ldb;
    INNER_RIGHT_LOOP_DEC( 6); B += ldb; INNER_RIGHT_LOOP_DEC( 7); B += ldb;
    INNER_RIGHT_LOOP_DEC( 8); B += ldb; INNER_RIGHT_LOOP_DEC( 9); B += ldb;
    INNER_RIGHT_LOOP_DEC(10); B += ldb; INNER_RIGHT_LOOP_DEC(11); B += ldb;
    INNER_RIGHT_LOOP_DEC(12); B += ldb; INNER_RIGHT_LOOP_DEC(13); B += ldb;
    INNER_RIGHT_LOOP_DEC(14); B += ldb; INNER_RIGHT_LOOP_DEC(15); B += ldb;

    __syncthreads();

    A += kb;
    k -= kb;
  }

  if (k > 0) { INNER_RIGHT_LOOP_DEC( 0); B += ldb;
  if (k > 1) { INNER_RIGHT_LOOP_DEC( 1); B += ldb;
  if (k > 2) { INNER_RIGHT_LOOP_DEC( 2); B += ldb;
  if (k > 3) { INNER_RIGHT_LOOP_DEC( 3); B += ldb;
  if (k > 4) { INNER_RIGHT_LOOP_DEC( 4); B += ldb;
  if (k > 5) { INNER_RIGHT_LOOP_DEC( 5); B += ldb;
  if (k > 6) { INNER_RIGHT_LOOP_DEC( 6); B += ldb;
  if (k > 7) { INNER_RIGHT_LOOP_DEC( 7); B += ldb;
  if (k > 8) { INNER_RIGHT_LOOP_DEC( 8); B += ldb;
  if (k > 9) { INNER_RIGHT_LOOP_DEC( 9); B += ldb;
  if (k >10) { INNER_RIGHT_LOOP_DEC(10); B += ldb;
  if (k >11) { INNER_RIGHT_LOOP_DEC(11); B += ldb;
  if (k >12) { INNER_RIGHT_LOOP_DEC(12); B += ldb;
  if (k >13) { INNER_RIGHT_LOOP_DEC(13); B += ldb;
  if (k >14) { INNER_RIGHT_LOOP_DEC(14); B += ldb;
  if (k >15) { INNER_RIGHT_LOOP_DEC(15); }}}}}}}}}}}}}}}}

  if (m - bi - ti > 0)
    sscal(n - bj, alpha, x, X, ldx);
