   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];
 }
 
 /**

 * 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, 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.
    if (diag == CBlasNonUnit) {
#pragma unroll
      for (int l = 0; l < kb; l++) {
        if (k + l >= ti)
          saxpy(A[0], b[l], x);
        A += lda;
      }
    }
    else {
#pragma unroll
      for (int l = 0; l < kb; l++) {
        if (k + l == ti)
          scopy(b[l], x);
        if (k + l > ti)
          saxpy(A[0], b[l], x);
        A += lda;
      }
    }

    __syncthreads();

    B += kb;
    k += kb;
  }

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

  if (n <= 0) return;
