 //   bi *= mb;
 //   bj *= nb;
 
-  const int bi = blockIdx.x * mb;
-  const int bj = blockIdx.y * nb;
-
+  // Using an SGEMM kernel, SSYRK is:
+  // C = alpha * A * B + beta * C
+  // with A = A and B = A' when trans == CBlasNoTrans, and
+  // with A = A' and B = A when trans == CBlasTrans
+  const float * __restrict__ B = A;
+
+  const int bi = blockIdx.x * mb;       // Starting row of block of C
+  const int bj = blockIdx.y * nb;       // Starting column of block of C
   const int ti = threadIdx.y * bx + threadIdx.x;        // Unwrapped thread index [0, bx * by]
 
+  if (uplo == CBlasUpper) {
+    if (bj + nb - 1 < bi)
+      return;
+  }
+  else if (uplo == CBlasLower) {
+    if (bi + mb - 1 < bj)
+      return;
+  }
+
   /*
-   * Cause blocks that are entirely above or below the diagonal to exit now.
-   */
-//   if (mb == nb && blockIdx.y < blockIdx.x) return;
-  if ((uplo == CBlasUpper && bj + nb - 1 < bi) ||
-      (uplo == CBlasLower && bi + mb - 1 < bj))
-    return;
-
-  /*
-   * If trans == CBlasNoTrans, A is A and B is A'. If trans != CBlasNoTrans, A
-   * is A' and B is A.
-   */
-  const float * __restrict__ B = A;
-
-  /*
-   * Compute our starting points in A, "B" and C.
+   * Compute our starting points in A, B and C.
    *
-   * For trans != CBlasNoTrans A is cached in shared memory so the unwrapped
+   * For trans == CBlasTrans A is cached in shared memory so the unwrapped
    * thread index can be re-wrapped around mb when calculating C.
    *
    * If trans == CBlasNoTrans then bx * by == mb (checked later on) so there
    * doesn't need to be a separate check for trans == CBlasNoTrans in
    * calculating the start of C here.
    */
-  A += (trans == CBlasNoTrans) ? bi + ti : (bi + threadIdx.y) * lda + threadIdx.x;
-  B += (trans == CBlasNoTrans) ? threadIdx.y * lda + bj + threadIdx.x : (bj + threadIdx.y) * lda + threadIdx.x;
-  C += (bx * by == mb) ? bj * ldc + bi + ti : (bj + 16 * (ti / mb)) * ldc + bi + ti % mb;
+  if (trans == CBlasNoTrans) {
+    A += bi + ti;
+    B += threadIdx.y * lda + bj + threadIdx.x;
+    C += bj * ldc + bi + ti;
+  }
+  else {
+    A += (bi + threadIdx.y) * lda + threadIdx.x;
+    B += (bj + threadIdx.y) * lda + threadIdx.x;
+    C += (bj + 16 * (ti / mb)) * ldc + bi + ti % mb;
+  }
 
   /*
-   * Blocks of A and "B" in shared memory and C in registers.
+   * Blocks of A and B in shared memory and C in registers.
    */
   __shared__ float a[mb][kb + 1];       // Optimised away when trans == CBlasNoTrans
   __shared__ float b[kb][(trans == CBlasNoTrans) ? nb : nb + 1];


  while (k > 0) {
    if (trans == CBlasNoTrans) {
//       typedef char x[(nb % bx == 0) ? 1 : -1];  // bx must be a multiple of nb
//       typedef char y[(kb % by == 0) ? 1 : -1];  // by must be a multiple of kb
#pragma unroll
      for (int l = 0; l < kb; l += by) {
#pragma unroll

      }
    }
    else {
      // C = aA'A + bC so read A into shared memory and transpose before reading
      // B into shared memory untransposed
#pragma unroll
      for (int l = 0; l < kb; l += bx) {
#pragma unroll
