 
   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
-  if (uplo == CBlasUpper && trans != CBlasNoTrans ||
-      uplo == CBlasLower && trans == CBlasNoTrans) {
-    int k = min(n - bj, nb);
-    while (k > 0) {
-      if (trans == CBlasNoTrans) {
-#pragma unroll
-        for (int j = 0; j < nb; j += by) {
-          a_hi[threadIdx.x][j + threadIdx.y] =
-            (diag != CBlasNonUnit && threadIdx.x == j + threadIdx.y) ? __double2hiint(1.0) : __double2hiint(A[j * lda]);
-          a_lo[threadIdx.x][j + threadIdx.y] =
-            (diag != CBlasNonUnit && threadIdx.x == j + threadIdx.y) ? __double2loint(1.0) : __double2loint(A[j * lda]);
-        }
-      }
-      else {
-#pragma unroll
-        for (int l = 0; l < kb; l += by) {
-          a_hi[l + threadIdx.y][threadIdx.x] =
-            (diag != CBlasNonUnit && threadIdx.x == l + threadIdx.y) ? __double2hiint(1.0) : __double2hiint(A[l * lda]);
-          a_lo[l + threadIdx.y][threadIdx.x] =
-            (diag != CBlasNonUnit && threadIdx.x == l + threadIdx.y) ? __double2loint(1.0) : __double2loint(A[l * lda]);
-        }
-      }
-
-      __syncthreads();
-
-      if (k < kb) break;
-
-// #pragma unroll
-//       for (int ll = 0; ll < kb; ll++) {
-//         daxpy(ll + 1, B[0], a[ll], x);
-        daxpy( 1, B[0], a_hi[0], a_lo[0], x); B += ldb;
-        daxpy( 2, B[0], a_hi[1], a_lo[1], x); B += ldb;
-        daxpy( 3, B[0], a_hi[2], a_lo[2], x); B += ldb;
-        daxpy( 4, B[0], a_hi[3], a_lo[3], x); B += ldb;
-        daxpy( 5, B[0], a_hi[4], a_lo[4], x); B += ldb;
-        daxpy( 6, B[0], a_hi[5], a_lo[5], x); B += ldb;
-        daxpy( 7, B[0], a_hi[6], a_lo[6], x); B += ldb;
-        daxpy( 8, B[0], a_hi[7], a_lo[7], x); B += ldb;
-//         B += ldb;
-//       }
-
-      __syncthreads();
-
-      A += (trans == CBlasNoTrans) ? kb : kb * lda;
-      k -= kb;
-    }
-
-    for (int ll = 0; ll < k; ll++) {
-      daxpy(ll + 1, B[0], a_hi[ll], a_lo[ll], x);
-      B += ldb;
-    }
-
-    __syncthreads();
-  }
-
   // Process non-diagonal blocks as for DGEMM
-  int k = (trans == CBlasNoTrans) ? ((uplo == CBlasUpper) ? bj : n - bj - nb)
-                                  : ((uplo == CBlasUpper) ? n - bj - nb : bj);
-  while (k > 0) {
-    if (trans == CBlasNoTrans) {
-#pragma unroll
-      for (int j = 0; j < nb; j += by) {
-        a_hi[threadIdx.x][j + threadIdx.y] = __double2hiint(A[j * lda]);
-        a_lo[threadIdx.x][j + threadIdx.y] = __double2loint(A[j * lda]);
-      }
-    }
-    else {
-#pragma unroll
-      for (int l = 0; l < kb; l += by) {
-        a_hi[l + threadIdx.y][threadIdx.x] = __double2hiint(A[l * lda]);
-        a_lo[l + threadIdx.y][threadIdx.x] = __double2loint(A[l * lda]);
-      }
-    }
-
-    __syncthreads();
-
-    if (k < kb) break;
+  int k = bj;
+  while (k > 0) {
+#pragma unroll
+    for (int j = 0; j < nb; j += by) {
+      a_hi[threadIdx.x][j + threadIdx.y] = __double2hiint(A[j * lda]);
+      a_lo[threadIdx.x][j + threadIdx.y] = __double2loint(A[j * lda]);
+    }
+
+    __syncthreads();
 
 #pragma unroll
     for (int l = 0; l < kb; l++) {


    __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];


  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 = bi;
  while (k > 0) {
#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) {
