  * per multiprocessor is limited by the register usage.
  */
 template void strmm<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 strmm<CBlasLeft,  CBlasUpper, CBlasNoTrans, CBlasUnit,    64, 16, 16, 16,  4>(int, int, float, const float * __restrict__, int, const float * __restrict__, int, float * __restrict__, int);
+template void strmm<CBlasLeft,  CBlasUpper, CBlasNoTrans, CBlasUnit,    64, 16, 16, 16,  4>(int, int, float, const float * __restrict__, int, const float * __restrict__, int, float * __restrict__, int);
 // template void strmm<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 strmm<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 strmm<CBlasLeft,  CBlasLower, CBlasNoTrans, CBlasNonUnit, 64, 16, 16, 16,  4>(int, int, float, const float * __restrict__, int, const float * __restrict__, int, float * __restrict__, int);

    if (k + kb > m) break;

    // Read A straight from global memory.
#pragma unroll
    for (int l = 0; l < kb; l++) {
      if (k + l >= ti)
        saxpy(A[0], b[l], x);
      A += lda;
    }

    __syncthreads();

    B += kb;
    k += kb;
  }

  int kk = m - k;
  for (int l = 0; l < kk; l++) {
    if (k + l >= ti)
      saxpy(A[0], b[l], x);
    A += lda;
  }

  if (n <= 0) return;

  y[ 4] += alpha * x[ 4]; y[ 5] += alpha * x[ 5]; y[ 6] += alpha * x[ 6]; y[ 7] += alpha * x[ 7];
  y[ 8] += alpha * x[ 8]; y[ 9] += alpha * x[ 9]; y[10] += alpha * x[10]; y[11] += alpha * x[11];
  y[12] += alpha * x[12]; y[13] += alpha * x[13]; y[14] += alpha * x[14]; y[15] += alpha * x[15];
}

// y(1:16) = x(1:16)
__device__ void scopy(const float * x, float * y) {
  y[ 0] = x[ 0]; y[ 1] = x[ 1]; y[ 2] = x[ 2]; y[ 3] = x[ 3];
  y[ 4] = x[ 4]; y[ 5] = x[ 5]; y[ 6] = x[ 6]; y[ 7] = x[ 7];
  y[ 8] = x[ 8]; y[ 9] = x[ 9]; y[10] = x[10]; y[11] = x[11];
  y[12] = x[12]; y[13] = x[13]; y[14] = x[14]; y[15] = x[15];
}

/**
