 
       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++;
-        }
+#pragma unroll
+      for (int ll = 0; ll < kb; ll++) {
+        saxpy(nb - ll, B[0], &a[ll][ll], &x[ll]);
+        B += ldb;
+//         l++;
       }
 
       __syncthreads();

                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) {
    int k = min(n - bj, nb);
    int l = 0;
    while (k > 0) {

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

//   if (trans == CBlasNoTrans) {
    A += /*(uplo == CBlasUpper) ? */(bj + threadIdx.y) * lda + threadIdx.x/* : (bj + threadIdx.y) * lda + bj + threadIdx.x*/;
    B += /*(uplo == CBlasUpper) ? */bi + ti/* : bj * ldb + bi + ti*/;
//   }
//   else {
//     A += (uplo == CBlasUpper) ? (bj + threadIdx.y) * lda + bj + threadIdx.x : threadIdx.y * lda + bj + threadIdx.x;
//     B += (uplo == CBlasUpper) ? bj * ldb + bi + ti : bi + ti;
//   }
  X += bj * ldx + bi + ti;

  __shared__ float a[kb][nb];
