   // For Upper/NoTrans and Lower/Trans process diagonal last
   k = min(n - bj, nb);
   while (k > 0) {
-
-    __syncthreads();
-
 #pragma unroll
     for (int l = 0; l < kb; l += by)
-      a[l + threadIdx.y][threadIdx.x] =
-          (diag != CBlasNonUnit && threadIdx.x == l + threadIdx.y) ? 1.0f : A[l * lda];
-
-    __syncthreads();
-
-    if (k < kb) break;
-
-    saxpy(16, B[0], &a[ 0][ 0], &x[ 0]); B += ldb;
-    saxpy(15, B[0], &a[ 1][ 1], &x[ 1]); B += ldb;
-    saxpy(14, B[0], &a[ 2][ 2], &x[ 2]); B += ldb;
-    saxpy(13, B[0], &a[ 3][ 3], &x[ 3]); B += ldb;
-    saxpy(12, B[0], &a[ 4][ 4], &x[ 4]); B += ldb;
-    saxpy(11, B[0], &a[ 5][ 5], &x[ 5]); B += ldb;
-    saxpy(10, B[0], &a[ 6][ 6], &x[ 6]); B += ldb;
-    saxpy( 9, B[0], &a[ 7][ 7], &x[ 7]); B += ldb;
-    saxpy( 8, B[0], &a[ 8][ 8], &x[ 8]); B += ldb;
-    saxpy( 7, B[0], &a[ 9][ 9], &x[ 9]); B += ldb;
-    saxpy( 6, B[0], &a[10][10], &x[10]); B += ldb;
-    saxpy( 5, B[0], &a[11][11], &x[11]); B += ldb;
-    saxpy( 4, B[0], &a[12][12], &x[12]); B += ldb;
-    saxpy( 3, B[0], &a[13][13], &x[13]); B += ldb;
-    saxpy( 2, B[0], &a[14][14], &x[14]); B += ldb;
-    saxpy( 1, B[0], &a[15][15], &x[15]); B += ldb;
+      a[l + threadIdx.y][threadIdx.x] = A[l * 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();
 


    __syncthreads();

    if (k < kb) break;

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


// y(1:n) += alpha * x(1:n)
__device__ void saxpy(int n, float alpha, const float * x, float * y) {
  if (n <= 0) return;
  y[ 0] += alpha * x[ 0]; if ( 1 >= n) return; y[ 1] += alpha * x[ 1]; if ( 2 >= n) return;
  y[ 2] += alpha * x[ 2]; if ( 3 >= n) return; y[ 3] += alpha * x[ 3]; if ( 4 >= n) return;
  y[ 4] += alpha * x[ 4]; if ( 5 >= n) return; y[ 5] += alpha * x[ 5]; if ( 6 >= n) return;
