     B += ldb;
   }
 
+  if (m - bi - ti > 0)
+    dscal(n - bj, alpha, x, X, ldx);
+}
+
+template <CBlasDiag diag,
+          unsigned int mb, unsigned int nb, unsigned int kb,
+          unsigned int bx, unsigned int by>
+__global__ void dtrmmRLN(int m, int n,
+                        double alpha, const double * __restrict__ A, int lda,
+                        const double * __restrict__ B, int ldb,
+                        double * __restrict__ X, int ldx) {
+
+  const int bi = blockIdx.x * mb;       // Starting row of block of X
+  const int bj = blockIdx.y * nb;       // Starting column of block of X
+  int ti = threadIdx.y * bx + threadIdx.x;
+
+  A += (bj + threadIdx.y) * lda + bj + threadIdx.x;
+  B += bj * ldb + bi + ti;
+  X += bj * ldx + bi + ti;
+
+  __shared__ double a[kb][nb];
+
+  double x[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+
+  // For Upper/Trans and Lower/NoTrans process diagonal first
+  int k = min(n - bj, nb);
+  while (k > 0) {
+#pragma unroll
+    for (int j = 0; j < nb; j += by)
+      a[threadIdx.x][j + threadIdx.y] = A[j * lda];
+
+    __syncthreads();
+
+    if (k < kb) break;
+
+    INNER_RIGHT_LOOP_INC( 1); B += ldb;
+    INNER_RIGHT_LOOP_INC( 2); B += ldb;
+    INNER_RIGHT_LOOP_INC( 3); B += ldb;
+    INNER_RIGHT_LOOP_INC( 4); B += ldb;
+    INNER_RIGHT_LOOP_INC( 5); B += ldb;
+    INNER_RIGHT_LOOP_INC( 6); B += ldb;
+    INNER_RIGHT_LOOP_INC( 7); B += ldb;
+    INNER_RIGHT_LOOP_INC( 8); B += ldb;
+
+    __syncthreads();
+
+    A += kb;
+    k -= kb;
+  }
+
+  if (k >  0) { INNER_RIGHT_LOOP_INC( 1); B += ldb;
+  if (k >  1) { INNER_RIGHT_LOOP_INC( 2); B += ldb;
+  if (k >  2) { INNER_RIGHT_LOOP_INC( 3); B += ldb;
+  if (k >  3) { INNER_RIGHT_LOOP_INC( 4); B += ldb;
+  if (k >  4) { INNER_RIGHT_LOOP_INC( 5); B += ldb;
+  if (k >  5) { INNER_RIGHT_LOOP_INC( 6); B += ldb;
+  if (k >  6) { INNER_RIGHT_LOOP_INC( 7); B += ldb;
+  if (k >  7) { INNER_RIGHT_LOOP_INC( 8); }}}}}}}}
+
+  // Process non-diagonal blocks as for DGEMM
+  k = n - bj - nb;
+  while (k > 0) {
+#pragma unroll
+    for (int j = 0; j < nb; j += by)
+      a[threadIdx.x][j + threadIdx.y] = A[j * lda];
+
+    __syncthreads();
+
+    if (k < kb) break;
+
+#pragma unroll
+    for (int l = 0; l < kb; l++) {
+      daxpy(B[0], a[l], x);
+      B += ldb;
+    }
+
+    __syncthreads();
+
+    A += kb;
+    k -= kb;
+  }
+
+  for (int l = 0; l < k; l++) {
+    daxpy(B[0], a[l], x);
+    B += ldb;
+  }
+
+  if (m - bi - ti > 0)
+    dscal(n - bj, alpha, x, X, ldx);
+}
+
+template <CBlasDiag diag,
+          unsigned int mb, unsigned int nb, unsigned int kb,
+          unsigned int bx, unsigned int by>
+__global__ void dtrmmRLT(int m, int n,
+                         double alpha, const double * __restrict__ A, int lda,
+                         const double * __restrict__ B, int ldb,
+                         double * __restrict__ X, int ldx) {
+
+  const int bi = blockIdx.x * mb;       // Starting row of block of X
+  const int bj = blockIdx.y * nb;       // Starting column of block of X
+  int ti = threadIdx.y * bx + threadIdx.x;
+
+  A += threadIdx.y * lda + bj + threadIdx.x;
+  B += bi + ti;
+  X += bj * ldx + bi + ti;
+
+  __shared__ double a[kb][nb];
+
+  double x[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+
+  // Process non-diagonal blocks as for DGEMM
+  int k = bj;
+  while (k > 0) {
+#pragma unroll
+    for (int l = 0; l < kb; l += by)
+      a[l + threadIdx.y][threadIdx.x] = A[l * lda];
+
+    __syncthreads();
+
+#pragma unroll
+    for (int l = 0; l < kb; l++) {
+      daxpy(B[0], a[l], x);
+      B += ldb;
+    }
+
+    __syncthreads();
+
+    A += kb * lda;
+    k -= kb;
+  }
+
   // For Upper/NoTrans and Lower/Trans process diagonal last
-  if (uplo == CBlasUpper && trans == CBlasNoTrans ||
-      uplo == CBlasLower && trans != CBlasNoTrans) {
-
-    __syncthreads();
-
-    int k = min(n - bj, nb);
-    while (k > 0) {
-      if (trans == CBlasNoTrans) {
-#pragma unroll
-        for (int j = 0; j < nb; j += by)
-          a[threadIdx.x][j + threadIdx.y] =
-            (diag != CBlasNonUnit && threadIdx.x == j + threadIdx.y) ? 1.0 : A[j * lda];
-      }
-      else {
-#pragma unroll
-        for (int l = 0; l < kb; l += by)
-          a[l + threadIdx.y][threadIdx.x] =
-            (diag != CBlasNonUnit && threadIdx.x == l + threadIdx.y) ? 1.0 : A[l * lda];
-      }
-
-      __syncthreads();
-
-      if (k < kb) break;
-
-// #pragma unroll
-//       for (int ll = 0; ll < kb; ll++) {
-//         daxpy(nb - ll, B[0], &a[ll][ll], &x[ll]);
-        daxpy(8, B[0], &a[0][0], &x[0]); B += ldb;
-        daxpy(7, B[0], &a[1][1], &x[1]); B += ldb;
-        daxpy(6, B[0], &a[2][2], &x[2]); B += ldb;
-        daxpy(5, B[0], &a[3][3], &x[3]); B += ldb;
-        daxpy(4, B[0], &a[4][4], &x[4]); B += ldb;
-        daxpy(3, B[0], &a[5][5], &x[5]); B += ldb;
-        daxpy(2, B[0], &a[6][6], &x[6]); B += ldb;
-        daxpy(1, B[0], &a[7][7], &x[7]); B += ldb;
-//         B += ldb;
-//       }
-
-      __syncthreads();
-
-      A += (trans == CBlasNoTrans) ? kb : kb * lda;
-      k -= kb;
-    }
-
-//     for (int ll = 0; ll < k; ll++) {
-//       daxpy(nb - ll, B[0], &a[ll][ll], &x[ll]);
-//       B += ldb;
-//     }
-    if (k > 0) { daxpy(8, B[0], &a[0][0], &x[0]); B += ldb;
-    if (k > 1) { daxpy(7, B[0], &a[1][1], &x[1]); B += ldb;
-    if (k > 2) { daxpy(6, B[0], &a[2][2], &x[2]); B += ldb;
-    if (k > 3) { daxpy(5, B[0], &a[3][3], &x[3]); B += ldb;
-    if (k > 4) { daxpy(4, B[0], &a[4][4], &x[4]); B += ldb;
-    if (k > 5) { daxpy(3, B[0], &a[5][5], &x[5]); B += ldb;
-    if (k > 6) { daxpy(2, B[0], &a[6][6], &x[6]); B += ldb;
-    if (k > 7) { daxpy(1, B[0], &a[7][7], &x[7]); B += ldb; }}}}}}}}
-  }
-
-  n -= bj;
-  m -= bi + ti;
-  if (n <= 0 || m <= 0) return;
-  X[0] = alpha * x[ 0]; if ( 1 >= n) return; X += ldx;
-  X[0] = alpha * x[ 1]; if ( 2 >= n) return; X += ldx;
-  X[0] = alpha * x[ 2]; if ( 3 >= n) return; X += ldx;
-  X[0] = alpha * x[ 3]; if ( 4 >= n) return; X += ldx;
-  X[0] = alpha * x[ 4]; if ( 5 >= n) return; X += ldx;
-  X[0] = alpha * x[ 5]; if ( 6 >= n) return; X += ldx;
-  X[0] = alpha * x[ 6]; if ( 7 >= n) return; X += ldx;
-  X[0] = alpha * x[ 7];
+  k = min(n - bj, nb);
+  while (k > 0) {
+#pragma unroll
+    for (int l = 0; l < kb; l += by)
+      a[l + threadIdx.y][threadIdx.x] = A[l * lda];
+
+    __syncthreads();
+
+    if (k < kb) break;
+
+    INNER_RIGHT_LOOP_DEC( 0); B += ldb; INNER_RIGHT_LOOP_DEC( 1); B += ldb;
+    INNER_RIGHT_LOOP_DEC( 2); B += ldb; INNER_RIGHT_LOOP_DEC( 3); B += ldb;
+    INNER_RIGHT_LOOP_DEC( 4); B += ldb; INNER_RIGHT_LOOP_DEC( 5); B += ldb;
+    INNER_RIGHT_LOOP_DEC( 6); B += ldb; INNER_RIGHT_LOOP_DEC( 7); B += ldb;
+
+    __syncthreads();
+
+    A += kb * lda;
+    k -= kb;
+  }
+
+  if (k > 0) { INNER_RIGHT_LOOP_DEC( 0); B += ldb;
+  if (k > 1) { INNER_RIGHT_LOOP_DEC( 1); B += ldb;
+  if (k > 2) { INNER_RIGHT_LOOP_DEC( 2); B += ldb;
+  if (k > 3) { INNER_RIGHT_LOOP_DEC( 3); B += ldb;
+  if (k > 4) { INNER_RIGHT_LOOP_DEC( 4); B += ldb;
+  if (k > 5) { INNER_RIGHT_LOOP_DEC( 5); B += ldb;
+  if (k > 6) { INNER_RIGHT_LOOP_DEC( 6); B += ldb;
+  if (k > 7) { INNER_RIGHT_LOOP_DEC( 7); }}}}}}}}
+
+  if (m - bi - ti > 0)
+    dscal(n - bj, alpha, x, X, ldx);
 }
 
 #endif
 
-template void dtrmm2L<CBlasUpper, CBlasNoTrans, CBlasUnit,    64,  8, 16, 16,  4>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2L<CBlasUpper, CBlasNoTrans, CBlasNonUnit, 64,  8, 16, 16,  4>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2L<CBlasUpper, CBlasTrans,   CBlasUnit,    32, 16,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2L<CBlasUpper, CBlasTrans,   CBlasNonUnit, 32, 16,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2L<CBlasLower, CBlasNoTrans, CBlasUnit,    64,  8, 16, 16,  4>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2L<CBlasLower, CBlasNoTrans, CBlasNonUnit, 64,  8, 16, 16,  4>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2L<CBlasLower, CBlasTrans,   CBlasUnit,    32, 16,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2L<CBlasLower, CBlasTrans,   CBlasNonUnit, 32, 16,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-
-template void dtrmm2R<CBlasUpper, CBlasNoTrans, CBlasUnit,    64,  8,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2R<CBlasUpper, CBlasNoTrans, CBlasNonUnit, 64,  8,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2R<CBlasUpper, CBlasTrans,   CBlasUnit,    64,  8,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2R<CBlasUpper, CBlasTrans,   CBlasNonUnit, 64,  8,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2R<CBlasLower, CBlasNoTrans, CBlasUnit,    64,  8,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2R<CBlasLower, CBlasNoTrans, CBlasNonUnit, 64,  8,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2R<CBlasLower, CBlasTrans,   CBlasUnit,    64,  8,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
-template void dtrmm2R<CBlasLower, CBlasTrans,   CBlasNonUnit, 64,  8,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmLUN<CBlasUnit,    64,  8, 16, 16,  4>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmLUN<CBlasNonUnit, 64,  8, 16, 16,  4>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmLUT<CBlasUnit,    32, 16,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmLUT<CBlasNonUnit, 32, 16,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmLLN<CBlasUnit,    64,  8, 16, 16,  4>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmLLN<CBlasNonUnit, 64,  8, 16, 16,  4>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmLLT<CBlasUnit,    32, 16,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmLLT<CBlasNonUnit, 32, 16,  8,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+
+template void dtrmmRUN<CBlasUnit,    64,  8, 16, 16,  4>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmRUN<CBlasNonUnit, 64,  8, 16, 16,  4>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmRUT<CBlasUnit,    64,  8, 16,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmRUT<CBlasNonUnit, 64,  8, 16,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmRLN<CBlasUnit,    64,  8, 16, 16,  4>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmRLN<CBlasNonUnit, 64,  8, 16, 16,  4>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmRLT<CBlasUnit,    64,  8, 16,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);
+template void dtrmmRLT<CBlasNonUnit, 64,  8, 16,  8,  8>(int, int, double, const double * __restrict__, int, const double * __restrict__, int, double * __restrict__, int);


    __syncthreads();

    if (k < kb) break;

    if (trans == CBlasNoTrans) {
#pragma unroll
      for (int l = 0; l < kb; l++) {
        daxpy(A[0], b_hi[l], b_lo[l], x);
        A += lda;
      }
    }
    else {
#pragma unroll
      for (int l = 0; l < kb; l++)
        daxpy(__hiloint2double(a_hi[ti][l], a_lo[ti][l]), &b_hi[l][tj], &b_lo[l][tj], x);
    }

    __syncthreads();

    B += kb;
    k -= kb;
  }

  if (trans == CBlasNoTrans) {
    for (int l = 0; l < k; l++) {
      daxpy(A[0], b_hi[l], b_lo[l], x);
      A += lda;
    }
  }
  else {
    for (int l = 0; l < k; l++)
      daxpy(__hiloint2double(a_hi[ti][l], a_lo[ti][l]), &b_hi[l][tj], &b_lo[l][tj], x);
  }

  // For Upper/Trans and Lower/NoTrans process diagonal last
  if (uplo == CBlasUpper && trans != CBlasNoTrans ||
      uplo == CBlasLower && trans == CBlasNoTrans) {

    __syncthreads();

    int k = min(m - bi, mb);
    int l = 0;
    while (k > 0) {
      if (trans != CBlasNoTrans) {
#pragma unroll
        for (int i = 0; i < mb; i += by) {
          a_hi[i + threadIdx.y][threadIdx.x] = __double2hiint(A[i * lda]);
          a_lo[i + threadIdx.y][threadIdx.x] = __double2loint(A[i * lda]);
        }
        A += kb;
      }

#pragma unroll
      for (int j = 0; j < nb; j += by) {
        b_hi[threadIdx.x][j + threadIdx.y] = __double2hiint(B[j * ldb]);
        b_lo[threadIdx.x][j + threadIdx.y] = __double2loint(B[j * ldb]);
      }

      __syncthreads();

      if (k < kb) break;

      if (diag == CBlasNonUnit) {
#pragma unroll
        for (int ll = 0; ll < kb; ll++) {
          if (ti >= l++)
            daxpy((trans == CBlasNoTrans) ? A[0]
                                          : __hiloint2double(a_hi[ti][ll], a_lo[ti][ll]),
                  (trans == CBlasNoTrans) ? b_hi[ll] : &b_hi[ll][tj],
                  (trans == CBlasNoTrans) ? b_lo[ll] : &b_lo[ll][tj], x);
          if (trans == CBlasNoTrans)
            A += lda;
        }
      }
      else {
#pragma unroll
        for (int ll = 0; ll < kb; ll++) {
          if (ti == l)
            daxpy(1.0, (trans == CBlasNoTrans) ? b_hi[ll] : &b_hi[ll][tj],
                       (trans == CBlasNoTrans) ? b_lo[ll] : &b_lo[ll][tj], x);
          else if (ti > l)
            daxpy((trans == CBlasNoTrans) ? A[0]
                                          : __hiloint2double(a_hi[ti][ll], a_lo[ti][ll]),
                  (trans == CBlasNoTrans) ? b_hi[ll] : &b_hi[ll][tj],
                  (trans == CBlasNoTrans) ? b_lo[ll] : &b_lo[ll][tj], x);
          if (trans == CBlasNoTrans)
            A += lda;
          l++;
        }
      }

      __syncthreads();

      B += kb;
      k -= kb;
    }

    if (diag == CBlasNonUnit) {
      for (int ll = 0; ll < k; ll++) {
        if (ti >= l++)
          daxpy((trans == CBlasNoTrans) ? A[0]
                                          : __hiloint2double(a_hi[ti][ll], a_lo[ti][ll]),
                (trans == CBlasNoTrans) ? b_hi[ll] : &b_hi[ll][tj],
                (trans == CBlasNoTrans) ? b_lo[ll] : &b_lo[ll][tj], x);
        if (trans == CBlasNoTrans)
          A += lda;
      }
    }
    else {
      for (int ll = 0; ll < k; ll++) {
        if (ti == l)
          daxpy(1.0, (trans == CBlasNoTrans) ? b_hi[ll] : &b_hi[ll][tj],
                     (trans == CBlasNoTrans) ? b_lo[ll] : &b_lo[ll][tj], x);
        else if (ti > l)
          daxpy((trans == CBlasNoTrans) ? A[0]
                                          : __hiloint2double(a_hi[ti][ll], a_lo[ti][ll]),
                (trans == CBlasNoTrans) ? b_hi[ll] : &b_hi[ll][tj],
                (trans == CBlasNoTrans) ? b_lo[ll] : &b_lo[ll][tj], x);
        if (trans == CBlasNoTrans)
          A += lda;
        l++;
      }
    }
  }

  n -= bj + tj;
  m -= bi + ti;
  if (n <= 0 || m <= 0) return;
  X[0] = alpha * x[ 0]; if ( 1 >= n) return; X += ldx;
  X[0] = alpha * x[ 1]; if ( 2 >= n) return; X += ldx;
  X[0] = alpha * x[ 2]; if ( 3 >= n) return; X += ldx;
  X[0] = alpha * x[ 3]; if ( 4 >= n) return; X += ldx;
  X[0] = alpha * x[ 4]; if ( 5 >= n) return; X += ldx;
  X[0] = alpha * x[ 5]; if ( 6 >= n) return; X += ldx;
  X[0] = alpha * x[ 6]; if ( 7 >= n) return; X += ldx;
  X[0] = alpha * x[ 7];
}

template <CBlasUplo uplo, CBlasTranspose trans, CBlasDiag diag,
          unsigned int mb, unsigned int nb, unsigned int kb,
          unsigned int bx, unsigned int by>
__global__ void dtrmm2R(int m, int n,
                        double alpha, const double * __restrict__ A, int lda, const double * __restrict__ B, int ldb,
                        double * __restrict__ X, int ldx) {

  const int bi = blockIdx.x * mb;       // Starting row of block of X
  const int bj = blockIdx.y * nb;       // Starting column of block of X
  int ti = threadIdx.y * bx + threadIdx.x;

  if (trans == CBlasNoTrans) {
    A += (uplo == CBlasUpper) ? (bj + threadIdx.y) * lda + threadIdx.x
                              : (bj + threadIdx.y) * lda + bj + threadIdx.x;
    B += (uplo == CBlasUpper) ? bi + ti : bj * ldb + bi + ti;
  }
  else {
    A += (uplo == CBlasUpper) ? (bj + threadIdx.y) * lda + bj + threadIdx.x
                              : threadIdx.y * lda + bj + threadIdx.x;
    B += (uplo == CBlasUpper) ? bj * ldb + bi + ti : bi + ti;
  }
  X += bj * ldx + bi + ti;

  __shared__ int a_hi[kb][nb];


    __syncthreads();

    A += kb;
    k -= kb;
  }

  // For Upper/NoTrans and Lower/Trans process diagonal last
  k = min(n - bj, nb);
  while (k > 0) {
#pragma unroll
    for (int j = 0; j < nb; j += by)
      a[threadIdx.x][j + threadIdx.y] = A[j * lda];

    __syncthreads();

    if (k < kb) break;

    INNER_RIGHT_LOOP_DEC( 0); B += ldb; INNER_RIGHT_LOOP_DEC( 1); B += ldb;
    INNER_RIGHT_LOOP_DEC( 2); B += ldb; INNER_RIGHT_LOOP_DEC( 3); B += ldb;
    INNER_RIGHT_LOOP_DEC( 4); B += ldb; INNER_RIGHT_LOOP_DEC( 5); B += ldb;
    INNER_RIGHT_LOOP_DEC( 6); B += ldb; INNER_RIGHT_LOOP_DEC( 7); B += ldb;

    __syncthreads();

    A += kb;
    k -= kb;
  }

  if (k > 0) { INNER_RIGHT_LOOP_DEC( 0); B += ldb;
  if (k > 1) { INNER_RIGHT_LOOP_DEC( 1); B += ldb;
  if (k > 2) { INNER_RIGHT_LOOP_DEC( 2); B += ldb;
  if (k > 3) { INNER_RIGHT_LOOP_DEC( 3); B += ldb;
  if (k > 4) { INNER_RIGHT_LOOP_DEC( 4); B += ldb;
  if (k > 5) { INNER_RIGHT_LOOP_DEC( 5); B += ldb;
  if (k > 6) { INNER_RIGHT_LOOP_DEC( 6); B += ldb;
  if (k > 7) { INNER_RIGHT_LOOP_DEC( 7); }}}}}}}}

  if (m - bi - ti > 0)
    dscal(n - bj, alpha, x, X, ldx);
}

template <CBlasDiag diag,
          unsigned int mb, unsigned int nb, unsigned int kb,
          unsigned int bx, unsigned int by>
__global__ void dtrmmRUT(int m, int n,
                         double alpha, const double * __restrict__ A, int lda,
                         const double * __restrict__ B, int ldb,
                         double * __restrict__ X, int ldx) {

  const int bi = blockIdx.x * mb;       // Starting row of block of X
  const int bj = blockIdx.y * nb;       // Starting column of block of X
  int ti = threadIdx.y * bx + threadIdx.x;


  A += (bj + threadIdx.y) * lda + bj + threadIdx.x;
  B += bj * ldb + bi + ti;
  X += bj * ldx + bi + ti;

  __shared__ double a[kb][nb];

  double x[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };

  // For Upper/Trans and Lower/NoTrans process diagonal first
  int k = min(n - bj, nb);
  while (k > 0) {
#pragma unroll
    for (int l = 0; l < kb; l += by)
      a[l + threadIdx.y][threadIdx.x] = A[l * lda];

    __syncthreads();

    if (k < kb) break;

    INNER_RIGHT_LOOP_INC( 1); B += ldb;
    INNER_RIGHT_LOOP_INC( 2); B += ldb;
    INNER_RIGHT_LOOP_INC( 3); B += ldb;
    INNER_RIGHT_LOOP_INC( 4); B += ldb;
    INNER_RIGHT_LOOP_INC( 5); B += ldb;
    INNER_RIGHT_LOOP_INC( 6); B += ldb;
    INNER_RIGHT_LOOP_INC( 7); B += ldb;
    INNER_RIGHT_LOOP_INC( 8); B += ldb;

    __syncthreads();

    A += kb * lda;
    k -= kb;
  }

  if (k >  0) { INNER_RIGHT_LOOP_INC( 1); B += ldb;
  if (k >  1) { INNER_RIGHT_LOOP_INC( 2); B += ldb;
  if (k >  2) { INNER_RIGHT_LOOP_INC( 3); B += ldb;
  if (k >  3) { INNER_RIGHT_LOOP_INC( 4); B += ldb;
  if (k >  4) { INNER_RIGHT_LOOP_INC( 5); B += ldb;
  if (k >  5) { INNER_RIGHT_LOOP_INC( 6); B += ldb;
  if (k >  6) { INNER_RIGHT_LOOP_INC( 7); B += ldb;
  if (k >  7) { INNER_RIGHT_LOOP_INC( 8); }}}}}}}}

  // Process non-diagonal blocks as for DGEMM
  k = n - bj - nb;
  while (k > 0) {
#pragma unroll
    for (int l = 0; l < kb; l += by)
      a[l + threadIdx.y][threadIdx.x] = A[l * lda];

    __syncthreads();

    if (k < kb) break;

#pragma unroll
    for (int l = 0; l < kb; l++) {
      daxpy(B[0], a[l], x);
      B += ldb;
    }

    __syncthreads();

    A += kb * lda;
    k -= kb;
  }

