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

//
//     __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 (trans == CBlasNoTrans) {
#pragma unroll

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

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