diff --git a/interface/omatcopy.c b/interface/omatcopy.c
index f06f40c601..e0fe6e740e 100644
--- a/interface/omatcopy.c
+++ b/interface/omatcopy.c
@@ -47,6 +47,36 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #define BlasNoTrans  0
 #define BlasTrans    1
 
+#ifdef SMP
+#ifndef GEMM_MULTITHREAD_THRESHOLD
+#define GEMM_MULTITHREAD_THRESHOLD 4
+#endif
+
+static int omatcopy_rn(BLASLONG m, BLASLONG n, BLASLONG dummy0, FLOAT alpha, FLOAT *a, BLASLONG lda,
+		       FLOAT *b, BLASLONG ldb, FLOAT *dummy1, BLASLONG dummy2, void *dummy3)
+{
+	return OMATCOPY_K_RN(m, n, alpha, a, lda, b, ldb);
+}
+
+static int omatcopy_rt(BLASLONG m, BLASLONG n, BLASLONG dummy0, FLOAT alpha, FLOAT *a, BLASLONG lda,
+		       FLOAT *b, BLASLONG ldb, FLOAT *dummy1, BLASLONG dummy2, void *dummy3)
+{
+	return OMATCOPY_K_RT(m, n, alpha, a, lda, b, ldb);
+}
+
+static int omatcopy_cn(BLASLONG m, BLASLONG n, BLASLONG dummy0, FLOAT alpha, FLOAT *a, BLASLONG lda,
+		       FLOAT *b, BLASLONG ldb, FLOAT *dummy1, BLASLONG dummy2, void *dummy3)
+{
+	return OMATCOPY_K_CN(n, m, alpha, a, lda, b, ldb);
+}
+
+static int omatcopy_ct(BLASLONG m, BLASLONG n, BLASLONG dummy0, FLOAT alpha, FLOAT *a, BLASLONG lda,
+		       FLOAT *b, BLASLONG ldb, FLOAT *dummy1, BLASLONG dummy2, void *dummy3)
+{
+	return OMATCOPY_K_CT(n, m, alpha, a, lda, b, ldb);
+}
+#endif
+
 OPENBLAS_EXPORT
 
 #ifndef CBLAS 
@@ -56,6 +86,9 @@ void NAME( char* ORDER, char* TRANS, blasint *rows, blasint *cols, FLOAT *alpha,
 	char Order, Trans;
 	int order=-1,trans=-1;
 	blasint info = -1;
+#ifdef SMP
+	int nthreads;
+#endif
 
 	Order = *ORDER;
 	Trans = *TRANS;
@@ -76,6 +109,9 @@ void CNAME(enum CBLAS_ORDER CORDER, enum CBLAS_TRANSPOSE CTRANS, blasint crows,
 	FLOAT   *alpha; 
 	int order=-1,trans=-1;
 	blasint info = -1;
+#ifdef SMP
+	int nthreads;
+#endif
 
 	if ( CORDER == CblasColMajor ) order = BlasColMajor; 
 	if ( CORDER == CblasRowMajor ) order = BlasRowMajor; 
@@ -115,6 +151,48 @@ void CNAME(enum CBLAS_ORDER CORDER, enum CBLAS_TRANSPOSE CTRANS, blasint crows,
 
 	if ((*rows == 0) || (*cols == 0)) return;
 
+#ifdef SMP
+	nthreads = num_cpu_avail(1);
+	/* Small copies stay on one thread */
+	if ((double)*rows * (double)*cols < 65536. * GEMM_MULTITHREAD_THRESHOLD) nthreads = 1;
+	/* In the transposed copies each thread writes a narrow band of every row
+	   (or column) of B. Keep the bands at least 32 elements (256 bytes of double,
+	   128 of float) wide, otherwise the threads keep writing into the same cache
+	   lines. */
+	if (trans == BlasTrans) {
+		blasint split = (order == BlasRowMajor) ? *rows : *cols;
+		if (nthreads > split / 32) nthreads = MAX(1, split / 32);
+	}
+
+	if (nthreads > 1) {
+#ifdef DOUBLE
+		int mode = BLAS_DOUBLE | BLAS_REAL;
+#else
+		int mode = BLAS_SINGLE | BLAS_REAL;
+#endif
+		/* Each thread gets a block of rows (row-major) or columns (column-major)
+		   of A, so the column-major wrappers swap m and n. For the transposed
+		   copies a block of B then starts width elements further along, not width
+		   rows or columns, which is what BLAS_TRANSB_T does. */
+		if (order == BlasRowMajor) {
+			if (trans == BlasNoTrans)
+				blas_level1_thread(mode, *rows, *cols, 0, alpha, a, *lda, b, *ldb, NULL, 0,
+						   (int (*)(void))omatcopy_rn, nthreads);
+			else
+				blas_level1_thread(mode | BLAS_TRANSB_T, *rows, *cols, 0, alpha, a, *lda, b, *ldb, NULL, 0,
+						   (int (*)(void))omatcopy_rt, nthreads);
+		} else {
+			if (trans == BlasNoTrans)
+				blas_level1_thread(mode, *cols, *rows, 0, alpha, a, *lda, b, *ldb, NULL, 0,
+						   (int (*)(void))omatcopy_cn, nthreads);
+			else
+				blas_level1_thread(mode | BLAS_TRANSB_T, *cols, *rows, 0, alpha, a, *lda, b, *ldb, NULL, 0,
+						   (int (*)(void))omatcopy_ct, nthreads);
+		}
+		return;
+	}
+#endif
+
 	if ( order == BlasColMajor )
 	{
 		if ( trans == BlasNoTrans )
diff --git a/utest/test_extensions/test_domatcopy.c b/utest/test_extensions/test_domatcopy.c
index 8869f7b453..bee271c8e8 100644
--- a/utest/test_extensions/test_domatcopy.c
+++ b/utest/test_extensions/test_domatcopy.c
@@ -34,7 +34,7 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include "utest/openblas_utest.h"
 #include "common.h"
 
-#define DATASIZE 100
+#define DATASIZE 600
 
 struct DATA_DOMATCOPY {
     double a_test[DATASIZE * DATASIZE];
@@ -414,6 +414,98 @@ CTEST(domatcopy, rowmajor_notrans_col_100_row_50)
     ASSERT_DBL_NEAR_TOL(0.0, norm, DOUBLE_EPS);
 }
 
+/**
+ * Fortran API specific test
+ * Test domatcopy by comparing it against refernce
+ * with the following options:
+ *
+ * Column Major
+ * Transposition
+ * Large enough to be split between threads
+ * alpha = 2.0
+ */
+CTEST(domatcopy, colmajor_trans_col_600_row_500)
+{
+    blasint m = 500, n = 600;
+    blasint lda = 500, ldb = 600;
+    char order = 'C';
+    char trans = 'T';
+    double alpha = 2.0;
+
+    double norm = check_domatcopy('F', order, trans, m, n, alpha, lda, ldb);
+
+    ASSERT_DBL_NEAR_TOL(0.0, norm, DOUBLE_EPS);
+}
+
+/**
+ * Fortran API specific test
+ * Test domatcopy by comparing it against refernce
+ * with the following options:
+ *
+ * Column Major
+ * Copy only
+ * Large enough to be split between threads
+ * alpha = 2.0
+ */
+CTEST(domatcopy, colmajor_notrans_col_600_row_500)
+{
+    blasint m = 500, n = 600;
+    blasint lda = 500, ldb = 500;
+    char order = 'C';
+    char trans = 'N';
+    double alpha = 2.0;
+
+    double norm = check_domatcopy('F', order, trans, m, n, alpha, lda, ldb);
+
+    ASSERT_DBL_NEAR_TOL(0.0, norm, DOUBLE_EPS);
+}
+
+/**
+ * Fortran API specific test
+ * Test domatcopy by comparing it against refernce
+ * with the following options:
+ *
+ * Row Major
+ * Transposition
+ * Large enough to be split between threads
+ * alpha = 2.0
+ */
+CTEST(domatcopy, rowmajor_trans_col_600_row_500)
+{
+    blasint m = 500, n = 600;
+    blasint lda = 600, ldb = 500;
+    char order = 'R';
+    char trans = 'T';
+    double alpha = 2.0;
+
+    double norm = check_domatcopy('F', order, trans, m, n, alpha, lda, ldb);
+
+    ASSERT_DBL_NEAR_TOL(0.0, norm, DOUBLE_EPS);
+}
+
+/**
+ * Fortran API specific test
+ * Test domatcopy by comparing it against refernce
+ * with the following options:
+ *
+ * Row Major
+ * Copy only
+ * Large enough to be split between threads
+ * alpha = 2.0
+ */
+CTEST(domatcopy, rowmajor_notrans_col_600_row_500)
+{
+    blasint m = 500, n = 600;
+    blasint lda = 600, ldb = 600;
+    char order = 'R';
+    char trans = 'N';
+    double alpha = 2.0;
+
+    double norm = check_domatcopy('F', order, trans, m, n, alpha, lda, ldb);
+
+    ASSERT_DBL_NEAR_TOL(0.0, norm, DOUBLE_EPS);
+}
+
 #ifndef NO_CBLAS
 /**
  * C API specific test