From f2ccc504551026c9430e4eeb4149e111341110f2 Mon Sep 17 00:00:00 2001 From: Abdulmalik Alquwayfili Date: Sun, 8 Nov 2026 08:22:38 +0300 Subject: [PATCH] Thread s/d omatcopy for large matrices omatcopy always ran on one thread. Split A into blocks of rows (row-major) or columns (column-major) with blas_level1_thread and run the kernel on each block. For the transposed copies BLAS_TRANSB_T moves B by the block width. Copies below 65536 * GEMM_MULTITHREAD_THRESHOLD elements stay on one thread. The transposed ones give each thread at least 32 rows or columns, since narrower bands make the threads write into the same cache lines of B. Add domatcopy tests that are large enough to use the threaded path. --- interface/omatcopy.c | 78 +++++++++++++++++++++ utest/test_extensions/test_domatcopy.c | 94 +++++++++++++++++++++++++- 2 files changed, 171 insertions(+), 1 deletion(-) 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