 
   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
+      // C = aAA' + bC so read B into shared memory and transpose leaving A
+      // untransposed in global memory
 #pragma unroll
       for (int l = 0; l < kb; l += by) {
 #pragma unroll

  }

  if (trans == CBlasNoTrans) {
    for (int l = 0; l < k; l++) {
      saxpy(A[0], b[l], c);
      A += lda;
    }
  }
  else {
    for (int l = 0; l < k; l++)
      saxpy(a[(bx * by == mb) ? ti : ti % mb][l],
            &b[l][(bx * by == mb) ? 0 : 16 * (ti / mb)], c);
  }

  const unsigned int i = (bx * by == mb) ? bi + ti : bi + ti % mb;
  const unsigned int j = (bx * by == mb) ? bj : bj + 16 * (ti / mb);
  if (i < n) {
    n -= j;
    if (n <= 0) return;
    if (beta == 0.0f) {
      if ((uplo == CBlasUpper && i <= j +  0) || (uplo == CBlasLower && i >= j +  0)) C[0] = alpha * c[ 0]; if ( 1 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  1) || (uplo == CBlasLower && i >= j +  1)) C[0] = alpha * c[ 1]; if ( 2 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  2) || (uplo == CBlasLower && i >= j +  2)) C[0] = alpha * c[ 2]; if ( 3 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  3) || (uplo == CBlasLower && i >= j +  3)) C[0] = alpha * c[ 3]; if ( 4 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  4) || (uplo == CBlasLower && i >= j +  4)) C[0] = alpha * c[ 4]; if ( 5 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  5) || (uplo == CBlasLower && i >= j +  5)) C[0] = alpha * c[ 5]; if ( 6 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  6) || (uplo == CBlasLower && i >= j +  6)) C[0] = alpha * c[ 6]; if ( 7 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  7) || (uplo == CBlasLower && i >= j +  7)) C[0] = alpha * c[ 7]; if ( 8 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  8) || (uplo == CBlasLower && i >= j +  8)) C[0] = alpha * c[ 8]; if ( 9 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  9) || (uplo == CBlasLower && i >= j +  9)) C[0] = alpha * c[ 9]; if (10 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j + 10) || (uplo == CBlasLower && i >= j + 10)) C[0] = alpha * c[10]; if (11 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j + 11) || (uplo == CBlasLower && i >= j + 11)) C[0] = alpha * c[11]; if (12 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j + 12) || (uplo == CBlasLower && i >= j + 12)) C[0] = alpha * c[12]; if (13 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j + 13) || (uplo == CBlasLower && i >= j + 13)) C[0] = alpha * c[13]; if (14 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j + 14) || (uplo == CBlasLower && i >= j + 14)) C[0] = alpha * c[14]; if (15 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j + 15) || (uplo == CBlasLower && i >= j + 15)) C[0] = alpha * c[15];
    }
    else {
      if ((uplo == CBlasUpper && i <= j +  0) || (uplo == CBlasLower && i >= j +  0)) C[0] = alpha * c[ 0] + beta * C[0]; if ( 1 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  1) || (uplo == CBlasLower && i >= j +  1)) C[0] = alpha * c[ 1] + beta * C[0]; if ( 2 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  2) || (uplo == CBlasLower && i >= j +  2)) C[0] = alpha * c[ 2] + beta * C[0]; if ( 3 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  3) || (uplo == CBlasLower && i >= j +  3)) C[0] = alpha * c[ 3] + beta * C[0]; if ( 4 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  4) || (uplo == CBlasLower && i >= j +  4)) C[0] = alpha * c[ 4] + beta * C[0]; if ( 5 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  5) || (uplo == CBlasLower && i >= j +  5)) C[0] = alpha * c[ 5] + beta * C[0]; if ( 6 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  6) || (uplo == CBlasLower && i >= j +  6)) C[0] = alpha * c[ 6] + beta * C[0]; if ( 7 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  7) || (uplo == CBlasLower && i >= j +  7)) C[0] = alpha * c[ 7] + beta * C[0]; if ( 8 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  8) || (uplo == CBlasLower && i >= j +  8)) C[0] = alpha * c[ 8] + beta * C[0]; if ( 9 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j +  9) || (uplo == CBlasLower && i >= j +  9)) C[0] = alpha * c[ 9] + beta * C[0]; if (10 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j + 10) || (uplo == CBlasLower && i >= j + 10)) C[0] = alpha * c[10] + beta * C[0]; if (11 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j + 11) || (uplo == CBlasLower && i >= j + 11)) C[0] = alpha * c[11] + beta * C[0]; if (12 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j + 12) || (uplo == CBlasLower && i >= j + 12)) C[0] = alpha * c[12] + beta * C[0]; if (13 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j + 13) || (uplo == CBlasLower && i >= j + 13)) C[0] = alpha * c[13] + beta * C[0]; if (14 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j + 14) || (uplo == CBlasLower && i >= j + 14)) C[0] = alpha * c[14] + beta * C[0]; if (15 >= n) return; C += ldc;
      if ((uplo == CBlasUpper && i <= j + 15) || (uplo == CBlasLower && i >= j + 15)) C[0] = alpha * c[15] + beta * C[0];
    }
  }
}

    if (k < kb) break;

    if (trans == CBlasNoTrans) {
      // Read A from global memory
#pragma unroll
      for (int l = 0; l < kb; l++) {
        saxpy(A[0], b[l], c);
