   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 };
 
-  int k = m - bi;
-  int l = 0;
-  while (k > 0) {
-#pragma unroll
-    for (int j = 0; j < nb; j += by)
-      b[threadIdx.x][j + threadIdx.y] = B[j * ldb];
-
-    __syncthreads();
-
-    if (k < kb) break;
+  // For CBlasUpper process diagonal first
+  if (uplo == CBlasUpper) {
+    int k = min(m, mb);
+    int l = 0;
+    while (k > 0) {
+#pragma unroll
+      for (int j = 0; j < nb; j += by)
+        b[threadIdx.x][j + threadIdx.y] = B[j * ldb];
+
+      __syncthreads();
+
+      if (k < kb) break;
+
+      if (diag == CBlasNonUnit) {
+#pragma unroll
+        for (int ll = 0; ll < kb; ll++) {
+          if (ti <= l++)
+            saxpy(A[0], b[ll], x);
+          A += lda;
+        }
+      }
+      else {
+#pragma unroll
+        for (int ll = 0; ll < kb; ll++) {
+          if (ti == l)
+            saxpy(1.0f, b[ll], x);
+          else if (ti < l)
+            saxpy(A[0], b[ll], x);
+          A += lda;
+          l++;
+        }
+      }
+
+      __syncthreads();
+
+      B += kb;
+      k -= kb;
+    }
 
     if (diag == CBlasNonUnit) {
-#pragma unroll
-      for (int ll = 0; ll < kb; ll++) {
+      for (int ll = 0; ll < k; ll++) {
         if (ti <= l++)
           saxpy(A[0], b[ll], x);
         A += lda;
       }
     }
     else {
-#pragma unroll
-      for (int ll = 0; ll < kb; ll++) {
+      for (int ll = 0; ll < k; ll++) {
         if (ti == l)
           saxpy(1.0f, b[ll], x);
         else if (ti < l)

  y[12] += alpha * x[12]; y[13] += alpha * x[13]; y[14] += alpha * x[14]; y[15] += alpha * x[15];
}

template <CBlasSide side, CBlasUplo uplo, CBlasTranspose trans, CBlasDiag diag,
          unsigned int mb, unsigned int nb, unsigned int kb,
          unsigned int bx, unsigned int by>

        l++;
      }
    }
  }

  // Process non-diagonal blocks as for SGEMM
  int k = (uplo == CBlasUpper) ? m - bi - mb : bi;
  while (k > 0) {
#pragma unroll
    for (int j = 0; j < nb; j += by)
      b[threadIdx.x][j + threadIdx.y] = B[j * ldb];

    __syncthreads();

    if (k < kb) break;

#pragma unroll
    for (int l = 0; l < kb; l++) {
      saxpy(A[0], b[l], x);
      A += lda;
    }

    __syncthreads();

