       }
     }
     else {
-//       typedef char x[(kb % bx == 0) ? 1 : -1];  // bx must be a multiple of kb
-//       typedef char y[(mb % by == 0) ? 1 : -1];  // by must be a multiple of mb
-//       typedef char z[(nb % by == 0) ? 1 : -1];  // by must be a multiple of nb
+      // 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


  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

  }

  if (trans == CBlasNoTrans) {
    // Read A from global memory
    for (int l = 0; l < k; l++) {
      saxpy(A[0], b[l], c);
      A += lda;
    }
  }
  else {
    // Read A' from shared memory
    for (int l = 0; l < k; l++)
      saxpy(a[ti % mb][l], &b[l][16 * (ti / mb)], c);
  }

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