                 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
-  if (uplo == CBlasUpper && trans != CBlasNoTrans ||
-      uplo == CBlasLower && trans == CBlasNoTrans) {
+//   if (uplo == CBlasUpper && trans != CBlasNoTrans ||
+//       uplo == CBlasLower && trans == CBlasNoTrans) {
+//     int k = min(n - bj, nb);
+//     int l = 0;
+//     while (k > 0) {
+//       if (trans == CBlasNoTrans) {
+// #pragma unroll
+//         for (int j = 0; j < nb; j += by)
+//           a[threadIdx.x][j + threadIdx.y] = A[j * lda];
+//       }
+//       else {
+// #pragma unroll
+//         for (int l = 0; l < kb; l += by)
+//           a[l + threadIdx.y][threadIdx.x] = A[l * lda];
+//       }
+//
+//       __syncthreads();
+//
+//       if (k < kb) break;
+//
+//       if (diag == CBlasNonUnit) {
+// #pragma unroll
+//         for (int ll = 0; ll < kb; ll++) {
+//           if (ti <= l++)
+//             saxpy(B[0], a[ll], x);
+//           B += ldb;
+//         }
+//       }
+//       else {
+// #pragma unroll
+//         for (int ll = 0; ll < kb; ll++) {
+//           if (ti == l)
+//             saxpy(1.0f, a[ll], x);
+//           else if (ti < l)
+//             saxpy(B[0], a[ll], x);
+//           B += ldb;
+//           l++;
+//         }
+//       }
+//
+//       __syncthreads();
+//
+//       A += (trans == CBlasNoTrans) ? kb : kb * lda;
+//       k -= kb;
+//     }
+//
+//     if (diag == CBlasNonUnit) {
+//       for (int ll = 0; ll < k; ll++) {
+//         if (ti <= l++)
+//           saxpy(B[0], a[ll], x);
+//         B += ldb;
+//       }
+//     }
+//     else {
+//       for (int ll = 0; ll < k; ll++) {
+//         if (ti == l)
+//           saxpy(1.0f, a[ll], x);
+//         else if (ti < l)
+//           saxpy(B[0], a[ll], x);
+//         B += ldb;
+//         l++;
+//       }
+//     }
+//
+//     __syncthreads();
+//   }
+//
+//   // Process non-diagonal blocks as for SGEMM
+//   int k = (trans == CBlasNoTrans) ? ((uplo == CBlasUpper) ? n - bj - nb : bj)
+//                                   : ((uplo == CBlasUpper) ? bj : n - bj - nb);
+//   while (k > 0) {
+//     if (trans == CBlasNoTrans) {
+// #pragma unroll
+//       for (int j = 0; j < nb; j += by)
+//         a[threadIdx.x][j + threadIdx.y] = A[j * lda];
+//     }
+//     else {
+// #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 += (trans == CBlasNoTrans) ? kb : kb * lda;
+//     k -= kb;
+//   }
+//
+//   for (int l = 0; l < k; l++) {
+//     saxpy(B[0], a[l], x);
+//     B += ldb;
+//   }
+//
+//   // For Upper/NoTrans and Lower/Trans process diagonal last
+//   if (uplo == CBlasUpper && trans == CBlasNoTrans ||
+//       uplo == CBlasLower && trans != CBlasNoTrans) {
+//
+//     __syncthreads();
+
     int k = min(n - bj, nb);
     int l = 0;
     while (k > 0) {


      if (k < kb) break;

      if (diag == CBlasNonUnit) {
#pragma unroll
        for (int ll = 0; ll < kb; ll++) {
          if (ti <= l++)
            saxpy(B[0], a[ll], x);
          B += ldb;
        }
      }
      else {
#pragma unroll
        for (int ll = 0; ll < kb; ll++) {
          if (ti == l)
            saxpy(1.0f, a[ll], x);
          else if (ti < l)
            saxpy(B[0], a[ll], x);
          B += ldb;
          l++;
        }
      }

      __syncthreads();

      k -= kb;
    }

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

  n -= bj;
  m -= bi + ti;
