From b439d63b872cc2a362380bdb26279aee56059f88 Mon Sep 17 00:00:00 2001 From: Abdulmalik Alquwayfili Date: Sun, 18 Oct 2026 19:42:26 +0300 Subject: [PATCH] Add an AVX DOMATCOPY_RT/CT kernel for x86_64 The double precision transposed omatcopy used the generic C kernels from kernel/arm on x86_64. This adds a blocked kernel that transposes 4x4 blocks in AVX registers, with a plain C fallback for targets without AVX. When the leading dimensions allow it, the first rows and columns are peeled so the blocked part starts on 32-byte boundaries, which keeps 16-byte aligned (malloc'd) buffers from running slower than before. --- kernel/x86_64/KERNEL | 3 + kernel/x86_64/domatcopy_t_avx.c | 136 ++++++++++++++++++++++++++++++++ 2 files changed, 139 insertions(+) create mode 100644 kernel/x86_64/domatcopy_t_avx.c diff --git a/kernel/x86_64/KERNEL b/kernel/x86_64/KERNEL index a8ba70a318..1c67d56c48 100644 --- a/kernel/x86_64/KERNEL +++ b/kernel/x86_64/KERNEL @@ -506,3 +506,6 @@ ZSUMKERNEL = zsum_sse2.S SOMATCOPY_RT = omatcopy_rt.c SOMATCOPY_CT = omatcopy_ct.c + +DOMATCOPY_RT = domatcopy_t_avx.c +DOMATCOPY_CT = domatcopy_t_avx.c diff --git a/kernel/x86_64/domatcopy_t_avx.c b/kernel/x86_64/domatcopy_t_avx.c new file mode 100644 index 0000000000..1d3e86db5b --- /dev/null +++ b/kernel/x86_64/domatcopy_t_avx.c @@ -0,0 +1,136 @@ +/*************************************************************************** +Copyright (c) 2026, The OpenBLAS Project +All rights reserved. +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: +1. Redistributions of source code must retain the above copyright +notice, this list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright +notice, this list of conditions and the following disclaimer in +the documentation and/or other materials provided with the +distribution. +3. Neither the name of the OpenBLAS project nor the names of +its contributors may be used to endorse or promote products +derived from this software without specific prior written permission. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE OPENBLAS PROJECT OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE +USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*****************************************************************************/ + +#include "common.h" +#include +#if defined(__AVX__) +#include +#endif + +/* + * Double precision transposed copy B := alpha * A^T for the OMATCOPY_RT + * (row-major, ROWM) and OMATCOPY_CT (column-major) kernels. + * + * Both orders do the same memory operation, b[j*ldb + i] = alpha * a[i*lda + j], + * the column-major one with rows and cols swapped. The matrix is walked in + * 32x32 tiles of 4x4 blocks, transposed in AVX registers. + */ + +#define TILE 32 + +static void dtrans(BLASLONG m, BLASLONG n, FLOAT alpha, FLOAT *a, BLASLONG lda, FLOAT *b, BLASLONG ldb) +{ + BLASLONG i, j, i0, j0, i1, j1; + BLASLONG m4 = m & ~3, n4 = n & ~3; +#if defined(__AVX__) + __m256d va = _mm256_set1_pd(alpha); +#endif + + for (i0 = 0; i0 < m4; i0 += TILE) { + i1 = (i0 + TILE < m4) ? i0 + TILE : m4; + for (j0 = 0; j0 < n4; j0 += TILE) { + j1 = (j0 + TILE < n4) ? j0 + TILE : n4; + for (i = i0; i < i1; i += 4) { + FLOAT *ap = a + i * lda; + FLOAT *bp = b + i; + for (j = j0; j < j1; j += 4) { +#if defined(__AVX__) + __m256d r0 = _mm256_mul_pd(_mm256_loadu_pd(ap + j), va); + __m256d r1 = _mm256_mul_pd(_mm256_loadu_pd(ap + lda + j), va); + __m256d r2 = _mm256_mul_pd(_mm256_loadu_pd(ap + 2 * lda + j), va); + __m256d r3 = _mm256_mul_pd(_mm256_loadu_pd(ap + 3 * lda + j), va); + __m256d t0 = _mm256_unpacklo_pd(r0, r1); + __m256d t1 = _mm256_unpackhi_pd(r0, r1); + __m256d t2 = _mm256_unpacklo_pd(r2, r3); + __m256d t3 = _mm256_unpackhi_pd(r2, r3); + _mm256_storeu_pd(bp + (j + 0) * ldb, _mm256_permute2f128_pd(t0, t2, 0x20)); + _mm256_storeu_pd(bp + (j + 1) * ldb, _mm256_permute2f128_pd(t1, t3, 0x20)); + _mm256_storeu_pd(bp + (j + 2) * ldb, _mm256_permute2f128_pd(t0, t2, 0x31)); + _mm256_storeu_pd(bp + (j + 3) * ldb, _mm256_permute2f128_pd(t1, t3, 0x31)); +#else + BLASLONG ii, jj; + for (ii = 0; ii < 4; ii++) + for (jj = 0; jj < 4; jj++) + bp[(j + jj) * ldb + ii] = alpha * ap[ii * lda + j + jj]; +#endif + } + } + } + } + for (i = 0; i < m4; i++) + for (j = n4; j < n; j++) + b[j * ldb + i] = alpha * a[i * lda + j]; + for (i = m4; i < m; i++) + for (j = 0; j < n; j++) + b[j * ldb + i] = alpha * a[i * lda + j]; +} + +/* + * Buffers from malloc are often only 16-byte aligned, and with them the + * blocked loop above runs slower than the generic kernel at some sizes. + * When a leading dimension is a multiple of 4, every row then has the same + * alignment, so do the first few rows and columns here and start the + * blocked part on 32-byte boundaries. + */ +static void dtrans_aligned(BLASLONG m, BLASLONG n, FLOAT alpha, FLOAT *a, BLASLONG lda, FLOAT *b, BLASLONG ldb) +{ + BLASLONG i, j, pi = 0, pj = 0; + + if ((ldb & 3) == 0) pi = ((32 - ((uintptr_t)b & 31)) & 31) / sizeof(FLOAT); + if ((lda & 3) == 0) pj = ((32 - ((uintptr_t)a & 31)) & 31) / sizeof(FLOAT); + if (pi > m) pi = m; + if (pj > n) pj = n; + + for (i = 0; i < pi; i++) + for (j = 0; j < n; j++) + b[j * ldb + i] = alpha * a[i * lda + j]; + for (i = pi; i < m; i++) + for (j = 0; j < pj; j++) + b[j * ldb + i] = alpha * a[i * lda + j]; + + dtrans(m - pi, n - pj, alpha, a + pi * lda + pj, lda, b + pj * ldb + pi, ldb); +} + +int CNAME(BLASLONG rows, BLASLONG cols, FLOAT alpha, FLOAT *a, BLASLONG lda, FLOAT *b, BLASLONG ldb) +{ + if (rows <= 0) return 0; + if (cols <= 0) return 0; +#ifdef ROWM + dtrans_aligned(rows, cols, alpha, a, lda, b, ldb); +#else + if (alpha == 0.0) { + /* same as the generic CT kernel: write zeros without reading A */ + BLASLONG i, j; + for (i = 0; i < cols; i++) + for (j = 0; j < rows; j++) + b[j * ldb + i] = 0.0; + return 0; + } + dtrans_aligned(cols, rows, alpha, a, lda, b, ldb); +#endif + return 0; +}