     k -= kb;
   }
 
-  if (diag == CBlasNonUnit) {
-    for (int ll = 0; ll < k; ll++) {
-      if (ti <= l++)
-        saxpy(A[0], b[ll], x);
-      A += lda;
-    }
-  }
-  else {
-    for (int ll = 0; ll < k; ll++) {
-      if (ti == l)
-        saxpy(1.0f, b[ll], x);
-      else if (ti < l)
-        saxpy(A[0], b[ll], x);
-      A += lda;
-      l++;
+  for (int l = 0; l < k; l++) {
+    saxpy(A[0], b[l], x);
+    A += lda;
+  }
+
+  // For CBlasLower process diagonal last
+  if (uplo == CBlasLower) {
+    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) {
+      for (int ll = 0; ll < k; ll++) {
+        if (ti >= l++)
+          saxpy(A[0], b[ll], x);
+        A += lda;
+      }
+    }
+    else {
+      for (int ll = 0; ll < k; ll++) {
+        if (ti == l)
+          saxpy(1.0f, b[ll], x);
+        else if (ti > l)
+          saxpy(A[0], b[ll], x);
+        A += lda;
+        l++;
+      }
     }
   }
 

template void strmm2<CBlasLeft,  CBlasUpper, CBlasNoTrans, CBlasNonUnit, 64, 16, 16, 16,  4>(int, int, float, const float * __restrict__, int, const float * __restrict__, int, float * __restrict__, int);
// template void strmm2<CBlasLeft,  CBlasUpper, CBlasTrans,   CBlasUnit,    32, 32,  8,  8,  8>(int, int, float, const float * __restrict__, int, const float * __restrict__, int, float * __restrict__, int);
// template void strmm2<CBlasLeft,  CBlasUpper, CBlasTrans,   CBlasNonUnit, 32, 32,  8,  8,  8>(int, int, float, const float * __restrict__, int, const float * __restrict__, int, float * __restrict__, int);
// template void strmm2<CBlasLeft,  CBlasLower, CBlasNoTrans, CBlasUnit,    64, 16, 16, 16,  4>(int, int, float, const float * __restrict__, int, const float * __restrict__, int, float * __restrict__, int);
// template void strmm2<CBlasLeft,  CBlasLower, CBlasNoTrans, CBlasNonUnit, 64, 16, 16, 16,  4>(int, int, float, const float * __restrict__, int, const float * __restrict__, int, float * __restrict__, int);
// template void strmm2<CBlasLeft,  CBlasLower, CBlasTrans,   CBlasUnit,    32, 32,  8,  8,  8>(int, int, float, const float * __restrict__, int, const float * __restrict__, int, float * __restrict__, int);
// template void strmm2<CBlasLeft,  CBlasLower, CBlasTrans,   CBlasNonUnit, 32, 32,  8,  8,  8>(int, int, float, const float * __restrict__, int, const float * __restrict__, int, float * __restrict__, int);
// template void strmm2<CBlasRight, CBlasUpper, CBlasNoTrans, CBlasUnit,    64, 16, 16, 16,  4>(int, int, float, const float * __restrict__, int, const float * __restrict__, int, float * __restrict__, int);

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

  A += (uplo == CBlasUpper) ? bi * lda + bi + ti : bi + ti;
  B += (uplo == CBlasUpper) ? (bj + threadIdx.y) * ldb + bi + threadIdx.x : (bj + threadIdx.y) * ldb + threadIdx.x;
  X += bj * ldx + bi + ti;

  __shared__ float b[kb][nb];
