diff --git a/kernel/power/KERNEL.POWER10 b/kernel/power/KERNEL.POWER10 index c009e33cf4..c4b0529858 100644 --- a/kernel/power/KERNEL.POWER10 +++ b/kernel/power/KERNEL.POWER10 @@ -3,6 +3,22 @@ #CGEMM_BETA = ../generic/zgemm_beta.c #ZGEMM_BETA = ../generic/zgemm_beta.c +# BGEMM (BFloat16-in, BFloat16-out) for POWER10 using plain MMA (xvbf16ger2pp). +# The kernel source is shared with SBGEMM; -DBGEMM is injected by Makefile.L3 +# to switch the SAVE_ACC macros to convert float32 accumulators back to BF16. +# UNROLL_M=16 (A-panel, packed by ncopy_16), UNROLL_N=8 (B-panel, packed by +# ncopy_8) — matches the 16x8 primary tile in sbgemm_kernel_power10.c. +BGEMM_BETA = ../generic/gemm_beta.c +BGEMMKERNEL = sbgemm_kernel_power10.c +BGEMMINCOPY = sbgemm_ncopy_16_power10.c +BGEMMITCOPY = sbgemm_tcopy_16_power10.c +BGEMMONCOPY = sbgemm_ncopy_8_power10.c +BGEMMOTCOPY = sbgemm_tcopy_8_power10.c +BGEMMINCOPYOBJ = bgemm_incopy$(TSUFFIX).$(SUFFIX) +BGEMMITCOPYOBJ = bgemm_itcopy$(TSUFFIX).$(SUFFIX) +BGEMMONCOPYOBJ = bgemm_oncopy$(TSUFFIX).$(SUFFIX) +BGEMMOTCOPYOBJ = bgemm_otcopy$(TSUFFIX).$(SUFFIX) + SBGEMM_BETA = ../generic/gemm_beta.c SBGEMMKERNEL = sbgemm_kernel_power10.c SBGEMMINCOPY = sbgemm_ncopy_16_power10.c diff --git a/kernel/power/sbgemm_kernel_power10.c b/kernel/power/sbgemm_kernel_power10.c index 33d744abdd..98794d7259 100644 --- a/kernel/power/sbgemm_kernel_power10.c +++ b/kernel/power/sbgemm_kernel_power10.c @@ -27,17 +27,15 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include "common.h" #include #if defined(BFLOAT16) && defined(BFLOAT16CONVERSION) -static float -bfloat16tof32 (bfloat16 f16) -{ - float result = 0; - unsigned short *q = (unsigned short *) (&result); +static float bfloat16tof32 (bfloat16 f16) { + float result = 0; + unsigned short *q = (unsigned short *) (&result); #if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ - q[0] = f16; + q[0] = f16; #else - q[1] = f16; + q[1] = f16; #endif - return result; + return result; } #define BF16TOF32(x) (bfloat16tof32(x)) @@ -46,8 +44,72 @@ bfloat16tof32 (bfloat16 f16) #endif typedef __vector unsigned char vec_t; +/* Under BGEMM, FLOAT resolves to bfloat16 (see common.h), but the MMA + * accumulators always produce float32. Declare internal float vector types + * explicitly so arithmetic stays in float32 regardless of FLOAT. */ +#ifdef BGEMM +typedef float v4sf_t __attribute__ ((vector_size (16))); +typedef float v2sf_t __attribute__ ((vector_size (8))); +/* Scalar float alpha derived from the bfloat16 alpha parameter. */ +#else typedef FLOAT v4sf_t __attribute__ ((vector_size (16))); typedef FLOAT v2sf_t __attribute__ ((vector_size (8))); +#endif + +#ifdef BGEMM +/* + * STORE4_BF16: C[0..3] += alpha * acc for a BF16 output matrix. + * + * All arithmetic in float32. One VSX conversion instruction for the store: + * xvcvspbf16 — convert float32 sum to BF16 (4 lanes, 1 instruction) + * + * The read of existing C is done element-wise with BF16TOF32 (scalar widening). + * This is the same approach used throughout the rest of the kernel for C access. + * + * Steps: + * 1. Read 4 existing BF16 values from C, widen to float32 (BF16TOF32) + * 2. Scale acc by float32 alpha, add to existing C — all in float32 + * 3. Convert float32 sum → BF16 via xvcvspbf16 (1 VSX instruction, 4 lanes) + * 4. Extract BF16 results via uint16 element indices 0,2,4,6 and store + */ +#define STORE4_BF16(ptr, acc_vec, alpha_f32) \ + do { \ + v4sf_t _c = { \ + BF16TOF32 ((ptr)[0]), BF16TOF32 ((ptr)[1]), \ + BF16TOF32 ((ptr)[2]), BF16TOF32 ((ptr)[3]) \ + }; \ + v4sf_t _sum = _c + (acc_vec) * (alpha_f32); \ + vec_t _conv = __builtin_vsx_xvcvspbf16 ((vec_t)_sum); \ + typedef __vector unsigned short vec_u16_t; \ + vec_u16_t _cv = (vec_u16_t) _conv; \ + (ptr)[0] = _cv[0]; \ + (ptr)[1] = _cv[2]; \ + (ptr)[2] = _cv[4]; \ + (ptr)[3] = _cv[6]; \ + } while (0) + +/* Same as STORE4_BF16 but for 2 BF16 lanes (n&2 remainder paths). */ +#define STORE2_BF16(ptr, acc_vec, alpha_f32) \ + do { \ + v4sf_t _c = { BF16TOF32 ((ptr)[0]), BF16TOF32 ((ptr)[1]), 0, 0 }; \ + v4sf_t _sum = _c + (acc_vec) * (alpha_f32); \ + vec_t _conv = __builtin_vsx_xvcvspbf16 ((vec_t)_sum); \ + typedef __vector unsigned short vec_u16_t; \ + vec_u16_t _cv = (vec_u16_t) _conv; \ + (ptr)[0] = _cv[0]; \ + (ptr)[1] = _cv[2]; \ + } while (0) + +/* Scalar float32 → BF16 using xvcvspbf16 on a 1-element vector. + * Used for single-element remainder paths (m&1, n&1 tails). */ +static inline bfloat16 f32_to_bf16_scalar (float f) +{ + v4sf_t v = { f, 0, 0, 0 }; + vec_t conv = __builtin_vsx_xvcvspbf16 ((vec_t)v); + /* Extract BF16 from high 16b of first slot via vector element access */ + return (((__vector unsigned short)conv)[0]); +} +#endif /* BGEMM */ /* * BFLOAT16 xvbf16ger2pp instruction needs 4×2 matrix of @@ -57,95 +119,163 @@ typedef FLOAT v2sf_t __attribute__ ((vector_size (8))); #define MERGE_HIGH(x, y) (vec_t) vec_mergeh ((vector short)x, (vector short)y) #define MERGE_LOW(x, y) (vec_t) vec_mergel ((vector short)x, (vector short)y) -#define SAVE_ACC(ACC, J) \ - __builtin_mma_disassemble_acc ((void *)result, ACC); \ - rowC = (v4sf_t *) &CO[0* ldc+J]; \ - rowC[0] += result[0] * alpha; \ - rowC = (v4sf_t *) &CO[1*ldc+J]; \ - rowC[0] += result[1] * alpha; \ - rowC = (v4sf_t *) &CO[2*ldc+J]; \ - rowC[0] += result[2] * alpha; \ - rowC = (v4sf_t *) &CO[3*ldc+J]; \ - rowC[0] += result[3] * alpha; -#define SAVE_ACC1(ACC, J) \ - __builtin_mma_disassemble_acc ((void *)result, ACC); \ - rowC = (v4sf_t *) &CO[4* ldc+J]; \ - rowC[0] += result[0] * alpha; \ - rowC = (v4sf_t *) &CO[5*ldc+J]; \ - rowC[0] += result[1] * alpha; \ - rowC = (v4sf_t *) &CO[6*ldc+J]; \ - rowC[0] += result[2] * alpha; \ - rowC = (v4sf_t *) &CO[7*ldc+J]; \ - rowC[0] += result[3] * alpha; -#define SAVE4x2_ACC(ACC, J) \ - __builtin_mma_disassemble_acc ((void *)result, ACC); \ - rowC = (v2sf_t *) &CO[0* ldc+J]; \ - rowC[0] += result[0] * alpha; \ - rowC = (v2sf_t *) &CO[1* ldc+J]; \ - rowC[0] += result[2] * alpha; \ - rowC = (v2sf_t *) &CO[2* ldc+J]; \ - rowC[0] += result[4] * alpha; \ - rowC = (v2sf_t *) &CO[3* ldc+J]; \ - rowC[0] += result[6] * alpha; -#define SAVE4x2_ACC1(ACC, J) \ - __builtin_mma_disassemble_acc ((void *)result, ACC); \ - rowC = (v2sf_t *) &CO[4* ldc+J]; \ - rowC[0] += result[0] * alpha; \ - rowC = (v2sf_t *) &CO[5* ldc+J]; \ - rowC[0] += result[2] * alpha; \ - rowC = (v2sf_t *) &CO[6* ldc+J]; \ - rowC[0] += result[4] * alpha; \ - rowC = (v2sf_t *) &CO[7* ldc+J]; \ - rowC[0] += result[6] * alpha; +#ifndef BGEMM +/* ---- SBGEMM: accumulator result is float32; C is float* ---- */ +#define SAVE_ACC(ACC, J) \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + rowC = (v4sf_t *) &CO[0* ldc+J]; \ + rowC[0] += result[0] * alpha; \ + rowC = (v4sf_t *) &CO[1*ldc+J]; \ + rowC[0] += result[1] * alpha; \ + rowC = (v4sf_t *) &CO[2*ldc+J]; \ + rowC[0] += result[2] * alpha; \ + rowC = (v4sf_t *) &CO[3*ldc+J]; \ + rowC[0] += result[3] * alpha; + +#define SAVE_ACC1(ACC, J) \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + rowC = (v4sf_t *) &CO[4* ldc+J]; \ + rowC[0] += result[0] * alpha; \ + rowC = (v4sf_t *) &CO[5*ldc+J]; \ + rowC[0] += result[1] * alpha; \ + rowC = (v4sf_t *) &CO[6*ldc+J]; \ + rowC[0] += result[2] * alpha; \ + rowC = (v4sf_t *) &CO[7*ldc+J]; \ + rowC[0] += result[3] * alpha; + +#define SAVE4x2_ACC(ACC, J) \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + rowC = (v2sf_t *) &CO[0* ldc+J]; \ + rowC[0] += result[0] * alpha; \ + rowC = (v2sf_t *) &CO[1* ldc+J]; \ + rowC[0] += result[2] * alpha; \ + rowC = (v2sf_t *) &CO[2* ldc+J]; \ + rowC[0] += result[4] * alpha; \ + rowC = (v2sf_t *) &CO[3* ldc+J]; \ + rowC[0] += result[6] * alpha; + +#define SAVE4x2_ACC1(ACC, J) \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + rowC = (v2sf_t *) &CO[4* ldc+J]; \ + rowC[0] += result[0] * alpha; \ + rowC = (v2sf_t *) &CO[5* ldc+J]; \ + rowC[0] += result[2] * alpha; \ + rowC = (v2sf_t *) &CO[6* ldc+J]; \ + rowC[0] += result[4] * alpha; \ + rowC = (v2sf_t *) &CO[7* ldc+J]; \ + rowC[0] += result[6] * alpha; - #define SAVE4x2_ACC_SCALAR(ACC) { \ - __builtin_mma_disassemble_acc ((void *)result, ACC); \ - res[0] = result[0] * alpha; \ - res[1] = result[1] * alpha; \ - res[2] = result[2] * alpha; \ - res[3] = result[3] * alpha; \ - CO[0 * ldc] += res[0][0]; \ - CO[1 * ldc] += res[1][0]; \ - CO[2 * ldc] += res[2][0]; \ - CO[3 * ldc] += res[3][0]; \ +#define SAVE4x2_ACC_SCALAR(ACC) { \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + res[0] = result[0] * alpha; \ + res[1] = result[1] * alpha; \ + res[2] = result[2] * alpha; \ + res[3] = result[3] * alpha; \ + CO[0 * ldc] += res[0][0]; \ + CO[1 * ldc] += res[1][0]; \ + CO[2 * ldc] += res[2][0]; \ + CO[3 * ldc] += res[3][0]; \ } - #define SAVE4x2_ACC1_SCALAR(ACC) { \ - __builtin_mma_disassemble_acc ((void *)result, ACC); \ - res[0] = result[0] * alpha; \ - res[1] = result[1] * alpha; \ - res[2] = result[2] * alpha; \ - res[3] = result[3] * alpha; \ - CO[4 * ldc] += res[0][0]; \ - CO[5 * ldc] += res[1][0]; \ - CO[6 * ldc] += res[2][0]; \ - CO[7 * ldc] += res[3][0]; \ +#define SAVE4x2_ACC1_SCALAR(ACC) { \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + res[0] = result[0] * alpha; \ + res[1] = result[1] * alpha; \ + res[2] = result[2] * alpha; \ + res[3] = result[3] * alpha; \ + CO[4 * ldc] += res[0][0]; \ + CO[5 * ldc] += res[1][0]; \ + CO[6 * ldc] += res[2][0]; \ + CO[7 * ldc] += res[3][0]; \ } -#define MMA __builtin_mma_xvbf16ger2pp +#define SAVE2x4_ACC(ACC, J) \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + rowC = (v4sf_t *) &CO[0* ldc+J]; \ + rowC[0] += result[0] * alpha; \ + rowC = (v4sf_t *) &CO[1* ldc+J]; \ + rowC[0] += result[1] * alpha; + +#else /* BGEMM: accumulator result is float32; C is bfloat16* */ + +/* Disassemble, scale by float alpha, convert each float32 lane to BF16 and + * store to the BF16 output row. CO is bfloat16*, J is column offset. */ +#define SAVE_ACC(ACC, J) \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + STORE4_BF16 (&CO[0 * ldc + (J)], result[0], falpha); \ + STORE4_BF16 (&CO[1 * ldc + (J)], result[1], falpha); \ + STORE4_BF16 (&CO[2 * ldc + (J)], result[2], falpha); \ + STORE4_BF16 (&CO[3 * ldc + (J)], result[3], falpha); \ + +#define SAVE_ACC1(ACC, J) \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + STORE4_BF16 (&CO[4 * ldc + (J)], result[0], falpha); \ + STORE4_BF16 (&CO[5 * ldc + (J)], result[1], falpha); \ + STORE4_BF16 (&CO[6 * ldc + (J)], result[2], falpha); \ + STORE4_BF16 (&CO[7 * ldc + (J)], result[3], falpha); \ -#define SAVE2x4_ACC(ACC, J) \ - __builtin_mma_disassemble_acc ((void *)result, ACC); \ - rowC = (v4sf_t *) &CO[0* ldc+J]; \ - rowC[0] += result[0] * alpha; \ - rowC = (v4sf_t *) &CO[1* ldc+J]; \ - rowC[0] += result[1] * alpha; +/* SAVE4x2_ACC: 2-wide B side — accumulator row i maps to C row i, cols J..J+1 */ +#define SAVE4x2_ACC(ACC, J) \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + STORE2_BF16 (&CO[0 * ldc + (J)], result[0], falpha); \ + STORE2_BF16 (&CO[1 * ldc + (J)], result[1], falpha); \ + STORE2_BF16 (&CO[2 * ldc + (J)], result[2], falpha); \ + STORE2_BF16 (&CO[3 * ldc + (J)], result[3], falpha); \ + +#define SAVE4x2_ACC1(ACC, J) \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + STORE2_BF16 (&CO[4 * ldc + (J)], result[0], falpha); \ + STORE2_BF16 (&CO[5 * ldc + (J)], result[1], falpha); \ + STORE2_BF16 (&CO[6 * ldc + (J)], result[2], falpha); \ + STORE2_BF16 (&CO[7 * ldc + (J)], result[3], falpha); \ + +#define SAVE4x2_ACC_SCALAR(ACC) \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + res[0] = result[0] * falpha; \ + res[1] = result[1] * falpha; \ + res[2] = result[2] * falpha; \ + res[3] = result[3] * falpha; \ + CO[0 * ldc] = f32_to_bf16_scalar (BF16TOF32 (CO[0 * ldc]) + res[0][0]); \ + CO[1 * ldc] = f32_to_bf16_scalar (BF16TOF32 (CO[1 * ldc]) + res[1][0]); \ + CO[2 * ldc] = f32_to_bf16_scalar (BF16TOF32 (CO[2 * ldc]) + res[2][0]); \ + CO[3 * ldc] = f32_to_bf16_scalar (BF16TOF32 (CO[3 * ldc]) + res[3][0]); + +#define SAVE4x2_ACC1_SCALAR(ACC) \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + res[0] = result[0] * falpha; \ + res[1] = result[1] * falpha; \ + res[2] = result[2] * falpha; \ + res[3] = result[3] * falpha; \ + CO[4 * ldc] = f32_to_bf16_scalar (BF16TOF32 (CO[4 * ldc]) + res[0][0]); \ + CO[5 * ldc] = f32_to_bf16_scalar (BF16TOF32 (CO[5 * ldc]) + res[1][0]); \ + CO[6 * ldc] = f32_to_bf16_scalar (BF16TOF32 (CO[6 * ldc]) + res[2][0]); \ + CO[7 * ldc] = f32_to_bf16_scalar (BF16TOF32 (CO[7 * ldc]) + res[3][0]); + +#define SAVE2x4_ACC(ACC, J) \ + __builtin_mma_disassemble_acc ((void *)result, ACC); \ + STORE4_BF16 (&CO[0 * ldc + (J)], result[0], falpha); \ + STORE4_BF16 (&CO[1 * ldc + (J)], result[1], falpha); \ + +#endif /* BGEMM */ + +/* MMA instruction is identical for both SBGEMM and BGEMM */ +#define MMA __builtin_mma_xvbf16ger2pp #define SET_ACC_ZERO4() \ - __builtin_mma_xxsetaccz (&acc0); \ - __builtin_mma_xxsetaccz (&acc1); \ - __builtin_mma_xxsetaccz (&acc2); \ - __builtin_mma_xxsetaccz (&acc3); + __builtin_mma_xxsetaccz (&acc0); \ + __builtin_mma_xxsetaccz (&acc1); \ + __builtin_mma_xxsetaccz (&acc2); \ + __builtin_mma_xxsetaccz (&acc3); #define SET_ACC_ZERO8() \ - __builtin_mma_xxsetaccz (&acc0); \ - __builtin_mma_xxsetaccz (&acc1); \ - __builtin_mma_xxsetaccz (&acc2); \ - __builtin_mma_xxsetaccz (&acc3); \ - __builtin_mma_xxsetaccz (&acc4); \ - __builtin_mma_xxsetaccz (&acc5); \ - __builtin_mma_xxsetaccz (&acc6); \ - __builtin_mma_xxsetaccz (&acc7); + __builtin_mma_xxsetaccz (&acc0); \ + __builtin_mma_xxsetaccz (&acc1); \ + __builtin_mma_xxsetaccz (&acc2); \ + __builtin_mma_xxsetaccz (&acc3); \ + __builtin_mma_xxsetaccz (&acc4); \ + __builtin_mma_xxsetaccz (&acc5); \ + __builtin_mma_xxsetaccz (&acc6); \ + __builtin_mma_xxsetaccz (&acc7); #define PREFETCH1(x, y) asm volatile ("dcbt %0, %1" : : "b" (x), "r" (y) : "memory"); /************************************************************************************* @@ -155,832 +285,805 @@ int CNAME (BLASLONG m, BLASLONG n, BLASLONG k, FLOAT alpha, IFLOAT * A, IFLOAT * B, FLOAT * C, BLASLONG ldc) { - BLASLONG i1; - v4sf_t valpha = { alpha, alpha, alpha, alpha }; - vector short vzero = { 0, 0, 0, 0, 0, 0, 0, 0 }; - /* Loop for n >= 8. */ - for (i1 = 0; i1 < (n >> 3); i1++) - { - BLASLONG j; - FLOAT *CO; - IFLOAT *AO; - CO = C; - C += ldc << 3; - AO = A; - PREFETCH1 (A, 128); - PREFETCH1 (A, 256); - /* Loop for m >= 16. */ - for (j = 0; j < (m >> 4); j++) - { - IFLOAT *BO = B; - v4sf_t *rowC; - v4sf_t result[4]; - __vector_quad acc0, acc1, acc2, acc3, acc4, acc5, acc6, acc7; - SET_ACC_ZERO8 (); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vec_t *rowA = (vec_t *) & (AO[l << 5]); - vec_t *rowB = (vec_t *) & (BO[l << 4]); - MMA (&acc0, rowB[0], rowA[0]); - MMA (&acc1, rowB[1], rowA[0]); - MMA (&acc2, rowB[0], rowA[1]); - MMA (&acc3, rowB[1], rowA[1]); - MMA (&acc4, rowB[0], rowA[2]); - MMA (&acc5, rowB[1], rowA[2]); - MMA (&acc6, rowB[0], rowA[3]); - MMA (&acc7, rowB[1], rowA[3]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 4; - vec_t *rowA = (vec_t *) & (AO[l << 1]); - vec_t *rowB = (vec_t *) & (BO[l]); - vec_t rowB_h = MERGE_HIGH (rowB[0], vzero); - vec_t rowB_l = MERGE_LOW (rowB[0], vzero); - vec_t rowA_h = MERGE_HIGH (rowA[0], vzero); - vec_t rowA_l = MERGE_LOW (rowA[0], vzero); - vec_t rowA2_h = MERGE_HIGH (rowA[1], vzero); - vec_t rowA2_l = MERGE_LOW (rowA[1], vzero); - MMA (&acc0, rowB_h, rowA_h); - MMA (&acc1, rowB_l, rowA_h); - MMA (&acc2, rowB_h, rowA_l); - MMA (&acc3, rowB_l, rowA_l); - MMA (&acc4, rowB_h, rowA2_h); - MMA (&acc5, rowB_l, rowA2_h); - MMA (&acc6, rowB_h, rowA2_l); - MMA (&acc7, rowB_l, rowA2_l); - } - SAVE_ACC (&acc0, 0); - SAVE_ACC (&acc2, 4); - SAVE_ACC1 (&acc1, 0); - SAVE_ACC1 (&acc3, 4); - SAVE_ACC (&acc4, 8); - SAVE_ACC (&acc6, 12); - SAVE_ACC1 (&acc5, 8); - SAVE_ACC1 (&acc7, 12); - CO += 16; - - AO += (k << 4); - BO += (k << 3); - } - if (m & 8) - { - IFLOAT *BO = B; - v4sf_t *rowC; - v4sf_t result[4]; - __vector_quad acc0, acc1, acc2, acc3; - SET_ACC_ZERO4 (); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vec_t *rowA = (vec_t *) & (AO[l << 4]); - vec_t *rowB = (vec_t *) & (BO[l << 4]); + BLASLONG i1; +#ifdef BGEMM + /* alpha is bfloat16 under BGEMM; convert once to float for all arithmetic. */ + float falpha = BF16TOF32 (alpha); + v4sf_t valpha = { falpha, falpha, falpha, falpha }; +#else + v4sf_t valpha = { alpha, alpha, alpha, alpha }; +#endif + vector short vzero = { 0, 0, 0, 0, 0, 0, 0, 0 }; + /* Loop for n >= 8. */ + for (i1 = 0; i1 < (n >> 3); i1++) { + BLASLONG j; + FLOAT *CO; + IFLOAT *AO; + CO = C; + C += ldc << 3; + AO = A; + PREFETCH1 (A, 128); + PREFETCH1 (A, 256); + /* Loop for m >= 16. */ + for (j = 0; j < (m >> 4); j++) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + __vector_quad acc0, acc1, acc2, acc3, acc4, acc5, acc6, acc7; + SET_ACC_ZERO8 (); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) + { + vec_t *rowA = (vec_t *) & (AO[l << 5]); + vec_t *rowB = (vec_t *) & (BO[l << 4]); + MMA (&acc0, rowB[0], rowA[0]); + MMA (&acc1, rowB[1], rowA[0]); + MMA (&acc2, rowB[0], rowA[1]); + MMA (&acc3, rowB[1], rowA[1]); + MMA (&acc4, rowB[0], rowA[2]); + MMA (&acc5, rowB[1], rowA[2]); + MMA (&acc6, rowB[0], rowA[3]); + MMA (&acc7, rowB[1], rowA[3]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 4; + vec_t *rowA = (vec_t *) & (AO[l << 1]); + vec_t *rowB = (vec_t *) & (BO[l]); + vec_t rowB_h = MERGE_HIGH (rowB[0], vzero); + vec_t rowB_l = MERGE_LOW (rowB[0], vzero); + vec_t rowA_h = MERGE_HIGH (rowA[0], vzero); + vec_t rowA_l = MERGE_LOW (rowA[0], vzero); + vec_t rowA2_h = MERGE_HIGH (rowA[1], vzero); + vec_t rowA2_l = MERGE_LOW (rowA[1], vzero); + MMA (&acc0, rowB_h, rowA_h); + MMA (&acc1, rowB_l, rowA_h); + MMA (&acc2, rowB_h, rowA_l); + MMA (&acc3, rowB_l, rowA_l); + MMA (&acc4, rowB_h, rowA2_h); + MMA (&acc5, rowB_l, rowA2_h); + MMA (&acc6, rowB_h, rowA2_l); + MMA (&acc7, rowB_l, rowA2_l); + } + SAVE_ACC (&acc0, 0); + SAVE_ACC (&acc2, 4); + SAVE_ACC1 (&acc1, 0); + SAVE_ACC1 (&acc3, 4); + SAVE_ACC (&acc4, 8); + SAVE_ACC (&acc6, 12); + SAVE_ACC1 (&acc5, 8); + SAVE_ACC1 (&acc7, 12); + CO += 16; + AO += (k << 4); + BO += (k << 3); + } + if (m & 8) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + __vector_quad acc0, acc1, acc2, acc3; + SET_ACC_ZERO4 (); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vec_t *rowA = (vec_t *) & (AO[l << 4]); + vec_t *rowB = (vec_t *) & (BO[l << 4]); - MMA (&acc0, rowB[0], rowA[0]); - MMA (&acc1, rowB[1], rowA[0]); - MMA (&acc2, rowB[0], rowA[1]); - MMA (&acc3, rowB[1], rowA[1]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 4; - vec_t *rowA = (vec_t *) & (AO[l]); - vec_t *rowB = (vec_t *) & (BO[l]); - vec_t rowB_h = MERGE_HIGH (rowB[0], vzero); - vec_t rowB_l = MERGE_LOW (rowB[0], vzero); - vec_t rowA_h = MERGE_HIGH (rowA[0], vzero); - vec_t rowA_l = MERGE_LOW (rowA[0], vzero); - MMA (&acc0, rowB_h, rowA_h); - MMA (&acc1, rowB_l, rowA_h); - MMA (&acc2, rowB_h, rowA_l); - MMA (&acc3, rowB_l, rowA_l); - } - SAVE_ACC (&acc0, 0); - SAVE_ACC (&acc2, 4); - SAVE_ACC1 (&acc1, 0); - SAVE_ACC1 (&acc3, 4); - CO += 8; - AO += (k << 3); - BO += (k << 3); - } - if (m & 4) - { - IFLOAT *BO = B; - v4sf_t *rowC; - v4sf_t result[4]; - __vector_quad acc0, acc1; - __builtin_mma_xxsetaccz (&acc0); - __builtin_mma_xxsetaccz (&acc1); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vec_t *rowA = (vec_t *) & (AO[l << 3]); - vec_t *rowB = (vec_t *) & (BO[l << 4]); - MMA (&acc0, rowB[0], rowA[0]); - MMA (&acc1, rowB[1], rowA[0]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 3; - vector short rowA = - { AO[l + 0], 0, AO[l + 1], 0, AO[l + 2], 0, AO[l + 3], 0 }; - vec_t *rowB = (vec_t *) & (BO[l << 1]); - MMA (&acc0, MERGE_HIGH (rowB[0], vzero), (vec_t) rowA); - MMA (&acc1, MERGE_LOW (rowB[0], vzero), (vec_t) rowA); - } - SAVE_ACC (&acc0, 0); - SAVE_ACC1 (&acc1, 0); - CO += 4; - AO += (k << 2); - BO += (k << 3); - } - if (m & 2) - { - IFLOAT *BO = B; - v2sf_t *rowC; - v2sf_t result[8]; - __vector_quad acc0, acc1; - __builtin_mma_xxsetaccz (&acc0); - __builtin_mma_xxsetaccz (&acc1); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vector short rowA = - { AO[(l << 2) + 0], AO[(l << 2) + 2], AO[(l << 2) + 1], - AO[(l << 2) + 3], - 0, 0, 0, 0 - }; - vec_t *rowB = (vec_t *) & (BO[l << 4]); - MMA (&acc0, rowB[0], (vec_t) rowA); - MMA (&acc1, rowB[1], (vec_t) rowA); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 2; - vector short rowA = { AO[l + 0], 0, AO[l + 1], 0, 0, 0, 0, 0 }; - vec_t *rowB = (vec_t *) & (BO[(l << 2)]); - MMA (&acc0, MERGE_HIGH (rowB[0], vzero), (vec_t) rowA); - MMA (&acc1, MERGE_LOW (rowB[0], vzero), (vec_t) rowA); - } - SAVE4x2_ACC (&acc0, 0); - SAVE4x2_ACC1 (&acc1, 0); - CO += 2; - AO += (k << 1); - BO += (k << 3); - } - if (m & 1) - { - IFLOAT *BO = B; - v4sf_t result[4], res[4]; - __vector_quad acc0, acc1; - __builtin_mma_xxsetaccz (&acc0); - __builtin_mma_xxsetaccz (&acc1); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vector short rowA = - { AO[(l << 1) + 0], AO[(l << 1) + 1], 0, 0, 0, 0, 0, 0}; - vec_t *rowB = (vec_t *) & (BO[l << 4]); - MMA (&acc0, rowB[0], (vec_t) rowA); - MMA (&acc1, rowB[1], (vec_t) rowA); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 1; - vector short rowA = { AO[l], 0, 0, 0, 0, 0, 0, 0 }; - vec_t *rowB = (vec_t *) & (BO[(l << 3)]); - MMA (&acc0, MERGE_HIGH (rowB[0], vzero), (vec_t) rowA); - MMA (&acc1, MERGE_LOW (rowB[0], vzero), (vec_t) rowA); - } - SAVE4x2_ACC_SCALAR (&acc0); - SAVE4x2_ACC1_SCALAR (&acc1); - CO += 1; - AO += k; - BO += (k << 3); - } - B += k << 3; + MMA (&acc0, rowB[0], rowA[0]); + MMA (&acc1, rowB[1], rowA[0]); + MMA (&acc2, rowB[0], rowA[1]); + MMA (&acc3, rowB[1], rowA[1]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 4; + vec_t *rowA = (vec_t *) & (AO[l]); + vec_t *rowB = (vec_t *) & (BO[l]); + vec_t rowB_h = MERGE_HIGH (rowB[0], vzero); + vec_t rowB_l = MERGE_LOW (rowB[0], vzero); + vec_t rowA_h = MERGE_HIGH (rowA[0], vzero); + vec_t rowA_l = MERGE_LOW (rowA[0], vzero); + MMA (&acc0, rowB_h, rowA_h); + MMA (&acc1, rowB_l, rowA_h); + MMA (&acc2, rowB_h, rowA_l); + MMA (&acc3, rowB_l, rowA_l); + } + SAVE_ACC (&acc0, 0); + SAVE_ACC (&acc2, 4); + SAVE_ACC1 (&acc1, 0); + SAVE_ACC1 (&acc3, 4); + CO += 8; + AO += (k << 3); + BO += (k << 3); + } + if (m & 4) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + __vector_quad acc0, acc1; + __builtin_mma_xxsetaccz (&acc0); + __builtin_mma_xxsetaccz (&acc1); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vec_t *rowA = (vec_t *) & (AO[l << 3]); + vec_t *rowB = (vec_t *) & (BO[l << 4]); + MMA (&acc0, rowB[0], rowA[0]); + MMA (&acc1, rowB[1], rowA[0]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 3; + vector short rowA = { AO[l + 0], 0, AO[l + 1], 0, AO[l + 2], 0, AO[l + 3], 0 }; + vec_t *rowB = (vec_t *) & (BO[l << 1]); + MMA (&acc0, MERGE_HIGH (rowB[0], vzero), (vec_t) rowA); + MMA (&acc1, MERGE_LOW (rowB[0], vzero), (vec_t) rowA); + } + SAVE_ACC (&acc0, 0); + SAVE_ACC1 (&acc1, 0); + CO += 4; + AO += (k << 2); + BO += (k << 3); + } + if (m & 2) { + IFLOAT *BO = B; +#ifndef BGEMM + v2sf_t *rowC; + v2sf_t result[8]; +#else + v4sf_t result[4]; +#endif + __vector_quad acc0, acc1; + __builtin_mma_xxsetaccz (&acc0); + __builtin_mma_xxsetaccz (&acc1); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vector short rowA = { AO[(l << 2) + 0], AO[(l << 2) + 2], AO[(l << 2) + 1], + AO[(l << 2) + 3], 0, 0, 0, 0 }; + vec_t *rowB = (vec_t *) & (BO[l << 4]); + MMA (&acc0, rowB[0], (vec_t) rowA); + MMA (&acc1, rowB[1], (vec_t) rowA); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 2; + vector short rowA = { AO[l + 0], 0, AO[l + 1], 0, 0, 0, 0, 0 }; + vec_t *rowB = (vec_t *) & (BO[(l << 2)]); + MMA (&acc0, MERGE_HIGH (rowB[0], vzero), (vec_t) rowA); + MMA (&acc1, MERGE_LOW (rowB[0], vzero), (vec_t) rowA); + } + SAVE4x2_ACC (&acc0, 0); + SAVE4x2_ACC1 (&acc1, 0); + CO += 2; + AO += (k << 1); + BO += (k << 3); + } + if (m & 1) { + IFLOAT *BO = B; + v4sf_t result[4], res[4]; + __vector_quad acc0, acc1; + __builtin_mma_xxsetaccz (&acc0); + __builtin_mma_xxsetaccz (&acc1); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vector short rowA = { AO[(l << 1) + 0], AO[(l << 1) + 1], 0, 0, 0, 0, 0, 0}; + vec_t *rowB = (vec_t *) & (BO[l << 4]); + MMA (&acc0, rowB[0], (vec_t) rowA); + MMA (&acc1, rowB[1], (vec_t) rowA); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 1; + vector short rowA = { AO[l], 0, 0, 0, 0, 0, 0, 0 }; + vec_t *rowB = (vec_t *) & (BO[(l << 3)]); + MMA (&acc0, MERGE_HIGH (rowB[0], vzero), (vec_t) rowA); + MMA (&acc1, MERGE_LOW (rowB[0], vzero), (vec_t) rowA); + } + SAVE4x2_ACC_SCALAR (&acc0); + SAVE4x2_ACC1_SCALAR (&acc1); + CO += 1; + AO += k; + BO += (k << 3); + } + B += k << 3; } - if (n & 4) - { - BLASLONG j; - FLOAT *CO; - IFLOAT *AO; - CO = C; - C += ldc << 2; - AO = A; - /* Loop for m >= 32. */ - for (j = 0; j < (m >> 5); j++) - { - IFLOAT *BO = B; - IFLOAT *A1 = AO + (16 * k); - v4sf_t *rowC; - v4sf_t result[4]; - __vector_quad acc0, acc1, acc2, acc3, acc4, acc5, acc6, acc7; - SET_ACC_ZERO8 (); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vec_t *rowA = (vec_t *) & (AO[l << 5]); - vec_t *rowA1 = (vec_t *) & (A1[l << 5]); - vec_t *rowB = (vec_t *) & (BO[l << 3]); - MMA (&acc0, rowB[0], rowA[0]); - MMA (&acc1, rowB[0], rowA[1]); - MMA (&acc2, rowB[0], rowA[2]); - MMA (&acc3, rowB[0], rowA[3]); - MMA (&acc4, rowB[0], rowA1[0]); - MMA (&acc5, rowB[0], rowA1[1]); - MMA (&acc6, rowB[0], rowA1[2]); - MMA (&acc7, rowB[0], rowA1[3]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 3; - vec_t *rowA = (vec_t *) & (AO[(l << 2)]); - vec_t *rowA1 = (vec_t *) & (A1[(l << 2)]); - vector short rowB_mrg = - { BO[l], 0, BO[l + 1], 0, BO[l + 2], 0, BO[l + 3], 0 }; - MMA (&acc0, (vec_t)rowB_mrg, MERGE_HIGH (rowA[0], vzero)); - MMA (&acc1, (vec_t)rowB_mrg, MERGE_LOW (rowA[0], vzero)); - MMA (&acc2, (vec_t)rowB_mrg, MERGE_HIGH (rowA[1], vzero)); - MMA (&acc3, (vec_t)rowB_mrg, MERGE_LOW (rowA[1], vzero)); - MMA (&acc4, (vec_t)rowB_mrg, MERGE_HIGH (rowA1[0], vzero)); - MMA (&acc5, (vec_t)rowB_mrg, MERGE_LOW (rowA1[0], vzero)); - MMA (&acc6, (vec_t)rowB_mrg, MERGE_HIGH (rowA1[1], vzero)); - MMA (&acc7, (vec_t)rowB_mrg, MERGE_LOW (rowA1[1], vzero)); - } - - SAVE_ACC (&acc0, 0); - SAVE_ACC (&acc1, 4); - CO += 8; - SAVE_ACC (&acc2, 0); - SAVE_ACC (&acc3, 4); - CO += 8; - SAVE_ACC (&acc4, 0); - SAVE_ACC (&acc5, 4); - CO += 8; - SAVE_ACC (&acc6, 0); - SAVE_ACC (&acc7, 4); - CO += 8; - AO += k << 5; - BO += k << 2; - } - if (m & 16) - { - IFLOAT *BO = B; - v4sf_t *rowC; - v4sf_t result[4]; - __vector_quad acc0, acc1, acc2, acc3; - SET_ACC_ZERO4 (); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vec_t *rowA = (vec_t *) & (AO[l << 5]); - vec_t *rowB = (vec_t *) & (BO[l << 3]); - MMA (&acc0, rowB[0], rowA[0]); - MMA (&acc1, rowB[0], rowA[1]); - MMA (&acc2, rowB[0], rowA[2]); - MMA (&acc3, rowB[0], rowA[3]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 3; - vec_t *rowA = (vec_t *) & (AO[(l << 2)]); - vector short rowB_mrg = - { BO[l], 0, BO[l + 1], 0, BO[l + 2], 0, BO[l + 3], 0 }; - MMA (&acc0, (vec_t)rowB_mrg, MERGE_HIGH (rowA[0], vzero)); - MMA (&acc1, (vec_t)rowB_mrg, MERGE_LOW (rowA[0], vzero)); - MMA (&acc2, (vec_t)rowB_mrg, MERGE_HIGH (rowA[1], vzero)); - MMA (&acc3, (vec_t)rowB_mrg, MERGE_LOW (rowA[1], vzero)); - } - - SAVE_ACC (&acc0, 0); - SAVE_ACC (&acc1, 4); - CO += 8; - SAVE_ACC (&acc2, 0); - SAVE_ACC (&acc3, 4); - CO += 8; - AO += k << 4; - BO += k << 2; - } - if (m & 8) - { - IFLOAT *BO = B; - v4sf_t *rowC; - v4sf_t result[4]; - __vector_quad acc0, acc1; - __builtin_mma_xxsetaccz (&acc0); - __builtin_mma_xxsetaccz (&acc1); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vec_t *rowA = (vec_t *) & (AO[l << 4]); - vec_t *rowB = (vec_t *) & (BO[l << 3]); - MMA (&acc0, rowB[0], rowA[0]); - MMA (&acc1, rowB[0], rowA[1]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 3; - vec_t *rowA = (vec_t *) & (AO[l << 1]); - vector short rowB_mrg = - { BO[l], 0, BO[l + 1], 0, BO[l + 2], 0, BO[l + 3], 0 }; - MMA (&acc0, (vec_t)rowB_mrg, MERGE_HIGH (rowA[0], vzero)); - MMA (&acc1, (vec_t)rowB_mrg, MERGE_LOW (rowA[0], vzero)); - } - SAVE_ACC (&acc0, 0); - SAVE_ACC (&acc1, 4); - CO += 8; - AO += k << 3; - BO += k << 2; - } - if (m & 4) - { - IFLOAT *BO = B; - v4sf_t *rowC; - __vector_quad acc0; - v4sf_t result[4]; - BLASLONG l = 0; - __builtin_mma_xxsetaccz (&acc0); - for (l = 0; l < k / 2; l++) - { - vec_t *rowA = (vec_t *) & (AO[l << 3]); - vec_t *rowB = (vec_t *) & (BO[l << 3]); - MMA (&acc0, rowB[0], rowA[0]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 3; - vector short rowA = - { AO[l], 0, AO[l + 1], 0, AO[l + 2], 0, AO[l + 3], 0 }; - vector short rowB_mrg = - { BO[l], 0, BO[l + 1], 0, BO[l + 2], 0, BO[l + 3], 0 }; - MMA (&acc0, (vec_t)(rowB_mrg), (vec_t) rowA); - } - SAVE_ACC (&acc0, 0); - CO += 4; - AO += k << 2; - BO += k << 2; - } - if (m & 2) - { - IFLOAT *BO = B; - v2sf_t *rowC; - v2sf_t result[8]; - __vector_quad acc0; - BLASLONG l = 0; - __builtin_mma_xxsetaccz (&acc0); - for (l = 0; l < k / 2; l++) - { - vector short rowA = - { AO[(l << 2) + 0], AO[(l << 2) + 2], AO[(l << 2) + 1], - AO[(l << 2) + 3], - 0, 0, 0, 0 - }; - vec_t *rowB = (vec_t *) & (BO[l << 3]); - MMA (&acc0, rowB[0], (vec_t) rowA); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 2; - vector short rowA = { AO[l], 0, AO[l + 1], 0, 0, 0, 0, 0 }; - vector short rowB_mrg = - { BO[(l<<1)], 0, BO[(l<<1) + 1], 0, BO[(l<<1) + 2], 0, - BO[(l<<1) + 3], 0 - }; - MMA (&acc0, (vec_t)(rowB_mrg), (vec_t) rowA); - } - SAVE4x2_ACC (&acc0, 0); - CO += 2; - AO += k << 1; - BO += k << 2; - } - if (m & 1) - { - IFLOAT *BO = B; - v4sf_t result[4], res[4]; - __vector_quad acc0; - BLASLONG l = 0; - __builtin_mma_xxsetaccz (&acc0); - for (l = 0; l < k / 2; l++) - { - vector short rowA = - { AO[(l << 1) + 0], AO[(l << 1) + 1], 0, - 0, 0, 0, 0 - }; - vec_t *rowB = (vec_t *) & (BO[l << 3]); - MMA (&acc0, rowB[0], (vec_t) rowA); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 1; - vector short rowA = { AO[l], 0, 0, 0, 0, 0, 0, 0 }; - vector short rowB_mrg = - { BO[(l<<2) + 0], 0, BO[(l<<2) + 1], 0, BO[(l <<2) + 2], 0, - BO[(l<<2) + 3], 0 - }; - MMA (&acc0, (vec_t)(rowB_mrg), (vec_t) rowA); - } - SAVE4x2_ACC_SCALAR (&acc0); - AO += k; - BO += (k << 2); - CO += 1; - } + if (n & 4) { + BLASLONG j; + FLOAT *CO; + IFLOAT *AO; + CO = C; + C += ldc << 2; + AO = A; + /* Loop for m >= 32. */ + for (j = 0; j < (m >> 5); j++) { + IFLOAT *BO = B; + IFLOAT *A1 = AO + (16 * k); +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + __vector_quad acc0, acc1, acc2, acc3, acc4, acc5, acc6, acc7; + SET_ACC_ZERO8 (); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vec_t *rowA = (vec_t *) & (AO[l << 5]); + vec_t *rowA1 = (vec_t *) & (A1[l << 5]); + vec_t *rowB = (vec_t *) & (BO[l << 3]); + MMA (&acc0, rowB[0], rowA[0]); + MMA (&acc1, rowB[0], rowA[1]); + MMA (&acc2, rowB[0], rowA[2]); + MMA (&acc3, rowB[0], rowA[3]); + MMA (&acc4, rowB[0], rowA1[0]); + MMA (&acc5, rowB[0], rowA1[1]); + MMA (&acc6, rowB[0], rowA1[2]); + MMA (&acc7, rowB[0], rowA1[3]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 3; + vec_t *rowA = (vec_t *) & (AO[(l << 2)]); + vec_t *rowA1 = (vec_t *) & (A1[(l << 2)]); + vector short rowB_mrg = { BO[l], 0, BO[l + 1], 0, BO[l + 2], 0, BO[l + 3], 0 }; + MMA (&acc0, (vec_t)rowB_mrg, MERGE_HIGH (rowA[0], vzero)); + MMA (&acc1, (vec_t)rowB_mrg, MERGE_LOW (rowA[0], vzero)); + MMA (&acc2, (vec_t)rowB_mrg, MERGE_HIGH (rowA[1], vzero)); + MMA (&acc3, (vec_t)rowB_mrg, MERGE_LOW (rowA[1], vzero)); + MMA (&acc4, (vec_t)rowB_mrg, MERGE_HIGH (rowA1[0], vzero)); + MMA (&acc5, (vec_t)rowB_mrg, MERGE_LOW (rowA1[0], vzero)); + MMA (&acc6, (vec_t)rowB_mrg, MERGE_HIGH (rowA1[1], vzero)); + MMA (&acc7, (vec_t)rowB_mrg, MERGE_LOW (rowA1[1], vzero)); + } + SAVE_ACC (&acc0, 0); + SAVE_ACC (&acc1, 4); + CO += 8; + SAVE_ACC (&acc2, 0); + SAVE_ACC (&acc3, 4); + CO += 8; + SAVE_ACC (&acc4, 0); + SAVE_ACC (&acc5, 4); + CO += 8; + SAVE_ACC (&acc6, 0); + SAVE_ACC (&acc7, 4); + CO += 8; + AO += k << 5; + BO += k << 2; + } + if (m & 16) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + __vector_quad acc0, acc1, acc2, acc3; + SET_ACC_ZERO4 (); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vec_t *rowA = (vec_t *) & (AO[l << 5]); + vec_t *rowB = (vec_t *) & (BO[l << 3]); + MMA (&acc0, rowB[0], rowA[0]); + MMA (&acc1, rowB[0], rowA[1]); + MMA (&acc2, rowB[0], rowA[2]); + MMA (&acc3, rowB[0], rowA[3]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 3; + vec_t *rowA = (vec_t *) & (AO[(l << 2)]); + vector short rowB_mrg = { BO[l], 0, BO[l + 1], 0, BO[l + 2], 0, BO[l + 3], 0 }; + MMA (&acc0, (vec_t)rowB_mrg, MERGE_HIGH (rowA[0], vzero)); + MMA (&acc1, (vec_t)rowB_mrg, MERGE_LOW (rowA[0], vzero)); + MMA (&acc2, (vec_t)rowB_mrg, MERGE_HIGH (rowA[1], vzero)); + MMA (&acc3, (vec_t)rowB_mrg, MERGE_LOW (rowA[1], vzero)); + } + SAVE_ACC (&acc0, 0); + SAVE_ACC (&acc1, 4); + CO += 8; + SAVE_ACC (&acc2, 0); + SAVE_ACC (&acc3, 4); + CO += 8; + AO += k << 4; + BO += k << 2; + } + if (m & 8) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + __vector_quad acc0, acc1; + __builtin_mma_xxsetaccz (&acc0); + __builtin_mma_xxsetaccz (&acc1); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vec_t *rowA = (vec_t *) & (AO[l << 4]); + vec_t *rowB = (vec_t *) & (BO[l << 3]); + MMA (&acc0, rowB[0], rowA[0]); + MMA (&acc1, rowB[0], rowA[1]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 3; + vec_t *rowA = (vec_t *) & (AO[l << 1]); + vector short rowB_mrg = { BO[l], 0, BO[l + 1], 0, BO[l + 2], 0, BO[l + 3], 0 }; + MMA (&acc0, (vec_t)rowB_mrg, MERGE_HIGH (rowA[0], vzero)); + MMA (&acc1, (vec_t)rowB_mrg, MERGE_LOW (rowA[0], vzero)); + } + SAVE_ACC (&acc0, 0); + SAVE_ACC (&acc1, 4); + CO += 8; + AO += k << 3; + BO += k << 2; + } + if (m & 4) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + __vector_quad acc0; + v4sf_t result[4]; + BLASLONG l = 0; + __builtin_mma_xxsetaccz (&acc0); + for (l = 0; l < k / 2; l++) { + vec_t *rowA = (vec_t *) & (AO[l << 3]); + vec_t *rowB = (vec_t *) & (BO[l << 3]); + MMA (&acc0, rowB[0], rowA[0]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 3; + vector short rowA = { AO[l], 0, AO[l + 1], 0, AO[l + 2], 0, AO[l + 3], 0 }; + vector short rowB_mrg = { BO[l], 0, BO[l + 1], 0, BO[l + 2], 0, BO[l + 3], 0 }; + MMA (&acc0, (vec_t)(rowB_mrg), (vec_t) rowA); + } + SAVE_ACC (&acc0, 0); + CO += 4; + AO += k << 2; + BO += k << 2; + } + if (m & 2) { + IFLOAT *BO = B; +#ifndef BGEMM + v2sf_t *rowC; + v2sf_t result[8]; +#else + v4sf_t result[4]; +#endif + __vector_quad acc0; + BLASLONG l = 0; + __builtin_mma_xxsetaccz (&acc0); + for (l = 0; l < k / 2; l++) { + vector short rowA = { AO[(l << 2) + 0], AO[(l << 2) + 2], AO[(l << 2) + 1], + AO[(l << 2) + 3], 0, 0, 0, 0 }; + vec_t *rowB = (vec_t *) & (BO[l << 3]); + MMA (&acc0, rowB[0], (vec_t) rowA); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 2; + vector short rowA = { AO[l], 0, AO[l + 1], 0, 0, 0, 0, 0 }; + vector short rowB_mrg = { BO[(l<<1)], 0, BO[(l<<1) + 1], 0, BO[(l<<1) + 2], 0, + BO[(l<<1) + 3], 0 }; + MMA (&acc0, (vec_t)(rowB_mrg), (vec_t) rowA); + } + SAVE4x2_ACC (&acc0, 0); + CO += 2; + AO += k << 1; + BO += k << 2; + } + if (m & 1) { + IFLOAT *BO = B; + v4sf_t result[4], res[4]; + __vector_quad acc0; + BLASLONG l = 0; + __builtin_mma_xxsetaccz (&acc0); + for (l = 0; l < k / 2; l++) { + vector short rowA = { AO[(l << 1) + 0], AO[(l << 1) + 1], 0, 0, 0, 0, 0 }; + vec_t *rowB = (vec_t *) & (BO[l << 3]); + MMA (&acc0, rowB[0], (vec_t) rowA); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 1; + vector short rowA = { AO[l], 0, 0, 0, 0, 0, 0, 0 }; + vector short rowB_mrg = { BO[(l<<2) + 0], 0, BO[(l<<2) + 1], 0, BO[(l <<2) + 2], 0, + BO[(l<<2) + 3], 0 }; + MMA (&acc0, (vec_t)(rowB_mrg), (vec_t) rowA); + } + SAVE4x2_ACC_SCALAR (&acc0); + AO += k; + BO += (k << 2); + CO += 1; + } B += k << 2; } - if (n & 2) - { - BLASLONG j; - FLOAT *CO; - IFLOAT *AO; - CO = C; - C += ldc << 1; - AO = A; - /* Loop for m >= 32. */ - for (j = 0; j < (m >> 5); j++) - { - IFLOAT *BO = B; - v4sf_t *rowC; - v4sf_t result[4]; - IFLOAT *A1 = AO + (16 * k); - __vector_quad acc0, acc1, acc2, acc3, acc4, acc5, acc6, acc7; - SET_ACC_ZERO8 (); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vector short rowB = - { BO[(l << 2) + 0], BO[(l << 2) + 2], BO[(l << 2) + 1], - BO[(l << 2) + 3], - 0, 0, 0, 0 - }; - vec_t *rowA = (vec_t *) & (AO[l << 5]); - vec_t *rowA1 = (vec_t *) & (A1[l << 5]); - MMA (&acc0, (vec_t) rowB, rowA[0]); - MMA (&acc1, (vec_t) rowB, rowA[1]); - MMA (&acc2, (vec_t) rowB, rowA[2]); - MMA (&acc3, (vec_t) rowB, rowA[3]); - MMA (&acc4, (vec_t) rowB, rowA1[0]); - MMA (&acc5, (vec_t) rowB, rowA1[1]); - MMA (&acc6, (vec_t) rowB, rowA1[2]); - MMA (&acc7, (vec_t) rowB, rowA1[3]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 2; - vector short rowB = { BO[l + 0], 0, BO[l + 1], 0, 0, 0, 0, 0 }; - vec_t *rowA = (vec_t *) & (AO[l << 3]); - vec_t *rowA1 = (vec_t *) & (A1[l << 3]); - MMA (&acc0, (vec_t) rowB, MERGE_HIGH (rowA[0], vzero)); - MMA (&acc1, (vec_t) rowB, MERGE_LOW (rowA[0], vzero)); - MMA (&acc2, (vec_t) rowB, MERGE_HIGH (rowA[1], vzero)); - MMA (&acc3, (vec_t) rowB, MERGE_LOW (rowA[1], vzero)); - MMA (&acc4, (vec_t) rowB, MERGE_HIGH (rowA1[0], vzero)); - MMA (&acc5, (vec_t) rowB, MERGE_LOW (rowA1[0], vzero)); - MMA (&acc6, (vec_t) rowB, MERGE_HIGH (rowA1[1], vzero)); - MMA (&acc7, (vec_t) rowB, MERGE_LOW (rowA1[1], vzero)); - } - SAVE2x4_ACC (&acc0, 0); - SAVE2x4_ACC (&acc1, 4); - SAVE2x4_ACC (&acc2, 8); - SAVE2x4_ACC (&acc3, 12); - CO += 16; - SAVE2x4_ACC (&acc4, 0); - SAVE2x4_ACC (&acc5, 4); - SAVE2x4_ACC (&acc6, 8); - SAVE2x4_ACC (&acc7, 12); - CO += 16; - AO += k << 5; - BO += k << 1; - } - if (m & 16) - { - IFLOAT *BO = B; - v4sf_t *rowC; - v4sf_t result[4]; - __vector_quad acc0, acc1, acc2, acc3; - SET_ACC_ZERO4 (); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vector short rowB = - { BO[(l << 2) + 0], BO[(l << 2) + 2], BO[(l << 2) + 1], - BO[(l << 2) + 3], - 0, 0, 0, 0 - }; - vec_t *rowA = (vec_t *) & (AO[l << 5]); - MMA (&acc0, (vec_t) rowB, rowA[0]); - MMA (&acc1, (vec_t) rowB, rowA[1]); - MMA (&acc2, (vec_t) rowB, rowA[2]); - MMA (&acc3, (vec_t) rowB, rowA[3]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 2; - vector short rowB = { BO[l + 0], 0, BO[l + 1], 0, 0, 0, 0, 0 }; - vec_t *rowA = (vec_t *) & (AO[l << 3]); - MMA (&acc0, (vec_t) rowB, MERGE_HIGH (rowA[0], vzero )); - MMA (&acc1, (vec_t) rowB, MERGE_LOW (rowA[0], vzero)); - MMA (&acc2, (vec_t) rowB, MERGE_HIGH (rowA[1], vzero)); - MMA (&acc3, (vec_t) rowB, MERGE_LOW (rowA[1], vzero)); - } - SAVE2x4_ACC (&acc0, 0); - SAVE2x4_ACC (&acc1, 4); - SAVE2x4_ACC (&acc2, 8); - SAVE2x4_ACC (&acc3, 12); - CO += 16; - AO += k << 4; - BO += k << 1; - } - if (m & 8) - { - IFLOAT *BO = B; - v4sf_t *rowC; - v4sf_t result[4]; - __vector_quad acc0, acc1; - __builtin_mma_xxsetaccz (&acc0); - __builtin_mma_xxsetaccz (&acc1); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vector short rowB = - { BO[(l << 2) + 0], BO[(l << 2) + 2], BO[(l << 2) + 1], - BO[(l << 2) + 3], - 0, 0, 0, 0 - }; - vec_t *rowA = (vec_t *) & (AO[l << 4]); - MMA (&acc0, (vec_t) rowB, rowA[0]); - MMA (&acc1, (vec_t) rowB, rowA[1]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 2; - vector short rowB = { BO[l + 0], 0, BO[l + 1], 0, 0, 0, 0, 0 }; - vec_t *rowA = (vec_t *) & (AO[(l << 2)]); - MMA (&acc0, (vec_t) rowB, MERGE_HIGH (rowA[0], vzero)); - MMA (&acc1, (vec_t) rowB, MERGE_LOW (rowA[0], vzero)); - } - SAVE2x4_ACC (&acc0, 0); - SAVE2x4_ACC (&acc1, 4); - CO += 8; - AO += k << 3; - BO += k << 1; - } - if (m & 4) - { - IFLOAT *BO = B; - v4sf_t *rowC; - v4sf_t result[4]; - __vector_quad acc0; - __builtin_mma_xxsetaccz (&acc0); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vector short rowB = - { BO[(l << 2) + 0], BO[(l << 2) + 2], BO[(l << 2) + 1], - BO[(l << 2) + 3], - 0, 0, 0, 0 - }; - vec_t *rowA = (vec_t *) & (AO[l << 3]); - MMA (&acc0, (vec_t) rowB, rowA[0]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 2; - vector short rowB = { BO[l + 0], 0, BO[l + 1], 0, 0, 0, 0, 0 }; - vector short rowA = - { AO[(l << 1)], 0, AO[(l << 1) + 1] , 0 , AO[(l<<1) + 2], - 0, AO[(l << 1) + 3], 0 }; - MMA (&acc0, (vec_t) rowB, (vec_t)(rowA)); - } - SAVE2x4_ACC (&acc0, 0); - CO += 4; - AO += k << 2; - BO += k << 1; - } - if (m & 2) - { - IFLOAT *BO = B; - BLASLONG l = 0; - v4sf_t t = { 0, 0, 0, 0 }; - for (l = 0; l < (k << 1); l += 2) - { - v4sf_t rowA = - { BF16TOF32 (AO[l]), BF16TOF32 (AO[l]), BF16TOF32 (AO[l + 1]), - BF16TOF32 (AO[l + 1]) - }; - v4sf_t rowB = - { BF16TOF32 (BO[l]), BF16TOF32 (BO[l + 1]), BF16TOF32 (BO[l]), - BF16TOF32 (BO[l + 1]) - }; - t += rowA * rowB; - } - t = t * valpha; - CO[0 * ldc] += t[0]; - CO[1 * ldc] += t[1]; - CO[0 * ldc + 1] += t[2]; - CO[1 * ldc + 1] += t[3]; - CO += 2; - AO += k << 1; - BO += k << 1; - } - if (m & 1) - { - IFLOAT *BO = B; - BLASLONG l = 0; - v4sf_t t = { 0, 0, 0, 0 }; - for (l = 0; l < k; l++) - { - v4sf_t rowA = { BF16TOF32 (AO[l]), BF16TOF32 (AO[l]), 0, 0 }; - v4sf_t rowB = - { BF16TOF32 (BO[l << 1]), BF16TOF32 (BO[(l << 1) + 1]), 0, - 0 - }; - t += rowA * rowB; - } - CO[0 * ldc] += t[0] * alpha; - CO[1 * ldc] += t[1] * alpha; - CO += 1; - AO += k; - BO += k << 1; - } - B += k << 1; + + if (n & 2) { + BLASLONG j; + FLOAT *CO; + IFLOAT *AO; + CO = C; + C += ldc << 1; + AO = A; + /* Loop for m >= 32. */ + for (j = 0; j < (m >> 5); j++) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + IFLOAT *A1 = AO + (16 * k); + __vector_quad acc0, acc1, acc2, acc3, acc4, acc5, acc6, acc7; + SET_ACC_ZERO8 (); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vector short rowB = { BO[(l << 2) + 0], BO[(l << 2) + 2], BO[(l << 2) + 1], + BO[(l << 2) + 3], 0, 0, 0, 0 }; + vec_t *rowA = (vec_t *) & (AO[l << 5]); + vec_t *rowA1 = (vec_t *) & (A1[l << 5]); + MMA (&acc0, (vec_t) rowB, rowA[0]); + MMA (&acc1, (vec_t) rowB, rowA[1]); + MMA (&acc2, (vec_t) rowB, rowA[2]); + MMA (&acc3, (vec_t) rowB, rowA[3]); + MMA (&acc4, (vec_t) rowB, rowA1[0]); + MMA (&acc5, (vec_t) rowB, rowA1[1]); + MMA (&acc6, (vec_t) rowB, rowA1[2]); + MMA (&acc7, (vec_t) rowB, rowA1[3]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 2; + vector short rowB = { BO[l + 0], 0, BO[l + 1], 0, 0, 0, 0, 0 }; + vec_t *rowA = (vec_t *) & (AO[l << 3]); + vec_t *rowA1 = (vec_t *) & (A1[l << 3]); + MMA (&acc0, (vec_t) rowB, MERGE_HIGH (rowA[0], vzero)); + MMA (&acc1, (vec_t) rowB, MERGE_LOW (rowA[0], vzero)); + MMA (&acc2, (vec_t) rowB, MERGE_HIGH (rowA[1], vzero)); + MMA (&acc3, (vec_t) rowB, MERGE_LOW (rowA[1], vzero)); + MMA (&acc4, (vec_t) rowB, MERGE_HIGH (rowA1[0], vzero)); + MMA (&acc5, (vec_t) rowB, MERGE_LOW (rowA1[0], vzero)); + MMA (&acc6, (vec_t) rowB, MERGE_HIGH (rowA1[1], vzero)); + MMA (&acc7, (vec_t) rowB, MERGE_LOW (rowA1[1], vzero)); + } + SAVE2x4_ACC (&acc0, 0); + SAVE2x4_ACC (&acc1, 4); + SAVE2x4_ACC (&acc2, 8); + SAVE2x4_ACC (&acc3, 12); + CO += 16; + SAVE2x4_ACC (&acc4, 0); + SAVE2x4_ACC (&acc5, 4); + SAVE2x4_ACC (&acc6, 8); + SAVE2x4_ACC (&acc7, 12); + CO += 16; + AO += k << 5; + BO += k << 1; + } + if (m & 16) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + __vector_quad acc0, acc1, acc2, acc3; + SET_ACC_ZERO4 (); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vector short rowB = { BO[(l << 2) + 0], BO[(l << 2) + 2], BO[(l << 2) + 1], + BO[(l << 2) + 3], 0, 0, 0, 0 }; + vec_t *rowA = (vec_t *) & (AO[l << 5]); + MMA (&acc0, (vec_t) rowB, rowA[0]); + MMA (&acc1, (vec_t) rowB, rowA[1]); + MMA (&acc2, (vec_t) rowB, rowA[2]); + MMA (&acc3, (vec_t) rowB, rowA[3]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 2; + vector short rowB = { BO[l + 0], 0, BO[l + 1], 0, 0, 0, 0, 0 }; + vec_t *rowA = (vec_t *) & (AO[l << 3]); + MMA (&acc0, (vec_t) rowB, MERGE_HIGH (rowA[0], vzero )); + MMA (&acc1, (vec_t) rowB, MERGE_LOW (rowA[0], vzero)); + MMA (&acc2, (vec_t) rowB, MERGE_HIGH (rowA[1], vzero)); + MMA (&acc3, (vec_t) rowB, MERGE_LOW (rowA[1], vzero)); + } + SAVE2x4_ACC (&acc0, 0); + SAVE2x4_ACC (&acc1, 4); + SAVE2x4_ACC (&acc2, 8); + SAVE2x4_ACC (&acc3, 12); + CO += 16; + AO += k << 4; + BO += k << 1; + } + if (m & 8) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + __vector_quad acc0, acc1; + __builtin_mma_xxsetaccz (&acc0); + __builtin_mma_xxsetaccz (&acc1); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vector short rowB = { BO[(l << 2) + 0], BO[(l << 2) + 2], BO[(l << 2) + 1], + BO[(l << 2) + 3], 0, 0, 0, 0 }; + vec_t *rowA = (vec_t *) & (AO[l << 4]); + MMA (&acc0, (vec_t) rowB, rowA[0]); + MMA (&acc1, (vec_t) rowB, rowA[1]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 2; + vector short rowB = { BO[l + 0], 0, BO[l + 1], 0, 0, 0, 0, 0 }; + vec_t *rowA = (vec_t *) & (AO[(l << 2)]); + MMA (&acc0, (vec_t) rowB, MERGE_HIGH (rowA[0], vzero)); + MMA (&acc1, (vec_t) rowB, MERGE_LOW (rowA[0], vzero)); + } + SAVE2x4_ACC (&acc0, 0); + SAVE2x4_ACC (&acc1, 4); + CO += 8; + AO += k << 3; + BO += k << 1; + } + if (m & 4) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + __vector_quad acc0; + __builtin_mma_xxsetaccz (&acc0); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vector short rowB = { BO[(l << 2) + 0], BO[(l << 2) + 2], BO[(l << 2) + 1], + BO[(l << 2) + 3], 0, 0, 0, 0 }; + vec_t *rowA = (vec_t *) & (AO[l << 3]); + MMA (&acc0, (vec_t) rowB, rowA[0]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 2; + vector short rowB = { BO[l + 0], 0, BO[l + 1], 0, 0, 0, 0, 0 }; + vector short rowA = + { AO[(l << 1)], 0, AO[(l << 1) + 1] , 0 , AO[(l<<1) + 2], + 0, AO[(l << 1) + 3], 0 }; + MMA (&acc0, (vec_t) rowB, (vec_t)(rowA)); + } + SAVE2x4_ACC (&acc0, 0); + CO += 4; + AO += k << 2; + BO += k << 1; + } + if (m & 2) { + IFLOAT *BO = B; + BLASLONG l = 0; + v4sf_t t = { 0, 0, 0, 0 }; + for (l = 0; l < (k << 1); l += 2) { + v4sf_t rowA = { BF16TOF32 (AO[l]), BF16TOF32 (AO[l]), BF16TOF32 (AO[l + 1]), + BF16TOF32 (AO[l + 1]) }; + v4sf_t rowB = { BF16TOF32 (BO[l]), BF16TOF32 (BO[l + 1]), BF16TOF32 (BO[l]), + BF16TOF32 (BO[l + 1]) }; + t += rowA * rowB; + } + t = t * valpha; +#ifdef BGEMM + CO[0 * ldc] = f32_to_bf16_scalar (BF16TOF32 (CO[0 * ldc]) + t[0]); + CO[1 * ldc] = f32_to_bf16_scalar (BF16TOF32 (CO[1 * ldc]) + t[1]); + CO[0 * ldc + 1] = f32_to_bf16_scalar (BF16TOF32 (CO[0 * ldc + 1]) + t[2]); + CO[1 * ldc + 1] = f32_to_bf16_scalar (BF16TOF32 (CO[1 * ldc + 1]) + t[3]); +#else + CO[0 * ldc] += t[0]; + CO[1 * ldc] += t[1]; + CO[0 * ldc + 1] += t[2]; + CO[1 * ldc + 1] += t[3]; +#endif + CO += 2; + AO += k << 1; + BO += k << 1; + } + if (m & 1) { + IFLOAT *BO = B; + BLASLONG l = 0; + v4sf_t t = { 0, 0, 0, 0 }; + for (l = 0; l < k; l++) { + v4sf_t rowA = { BF16TOF32 (AO[l]), BF16TOF32 (AO[l]), 0, 0 }; + v4sf_t rowB = { BF16TOF32 (BO[l << 1]), BF16TOF32 (BO[(l << 1) + 1]), 0, 0 }; + t += rowA * rowB; + } +#ifdef BGEMM + CO[0 * ldc] = f32_to_bf16_scalar (BF16TOF32 (CO[0 * ldc]) + t[0] * falpha); + CO[1 * ldc] = f32_to_bf16_scalar (BF16TOF32 (CO[1 * ldc]) + t[1] * falpha); +#else + CO[0 * ldc] += t[0] * alpha; + CO[1 * ldc] += t[1] * alpha; +#endif + CO += 1; + AO += k; + BO += k << 1; + } + B += k << 1; } - if (n & 1) - { - BLASLONG j; - FLOAT *CO; - IFLOAT *AO; - CO = C; - C += ldc; - AO = A; - /* Loop for m >= 16. */ - for (j = 0; j < (m >> 4); j++) - { - IFLOAT *BO = B; - v4sf_t *rowC; - v4sf_t result[4]; - __vector_quad acc0, acc1, acc2, acc3; - SET_ACC_ZERO4 (); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vector short rowB = - { BO[l << 1], BO[(l << 1) + 1], 0, 0, 0, 0, 0, 0}; - vec_t *rowA = (vec_t *) & (AO[l << 5]); - MMA (&acc0, (vec_t) rowB, rowA[0]); - MMA (&acc1, (vec_t) rowB, rowA[1]); - MMA (&acc2, (vec_t) rowB, rowA[2]); - MMA (&acc3, (vec_t) rowB, rowA[3]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 1; - vector short rowB = { BO[l], 0, 0, 0, 0, 0, 0, 0 }; - vec_t *rowA = (vec_t *) & (AO[(l << 4)]); - MMA (&acc0, (vec_t) rowB, MERGE_HIGH (rowA[0], vzero)); - MMA (&acc1, (vec_t) rowB, MERGE_LOW (rowA[0], vzero)); - MMA (&acc2, (vec_t) rowB, MERGE_HIGH (rowA[1], vzero)); - MMA (&acc3, (vec_t) rowB, MERGE_LOW (rowA[1], vzero)); - } - rowC = (v4sf_t *) &CO[0]; - __builtin_mma_disassemble_acc ((void *)result, &acc0); - rowC[0] += result[0] * alpha; - __builtin_mma_disassemble_acc ((void *)result, &acc1); - rowC[1] += result[0] * alpha; - __builtin_mma_disassemble_acc ((void *)result, &acc2); - rowC[2] += result[0] * alpha; - __builtin_mma_disassemble_acc ((void *)result, &acc3); - rowC[3] += result[0] * alpha; - AO += k << 4; - BO += k; - CO += 16; - } - /* Loop for m >= 8. */ - if (m & 8) - { - IFLOAT *BO = B; - v4sf_t *rowC; - v4sf_t result[4]; - __vector_quad acc0, acc1; - __builtin_mma_xxsetaccz (&acc0); - __builtin_mma_xxsetaccz (&acc1); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vector short rowB = - { BO[l << 1], BO[(l << 1) + 1], 0, 0, 0, 0, 0, 0}; - vec_t *rowA = (vec_t *) & (AO[l << 4]); - MMA (&acc0, (vec_t) rowB, rowA[0]); - MMA (&acc1, (vec_t) rowB, rowA[1]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 1; - vector short rowB = { BO[l], 0, 0, 0, 0, 0, 0, 0 }; - vec_t *rowA = (vec_t *) & (AO[(l << 3)]); - MMA (&acc0, (vec_t) rowB, MERGE_HIGH (rowA[0], vzero)); - MMA (&acc1, (vec_t) rowB, MERGE_LOW (rowA[0], vzero)); - } - rowC = (v4sf_t *) &CO[0]; - __builtin_mma_disassemble_acc ((void *)result, &acc0); - rowC[0] += result[0] * alpha; - __builtin_mma_disassemble_acc ((void *)result, &acc1); - rowC[1] += result[0] * alpha; - AO += k << 3; - BO += k; - CO += 8; - } - /* Loop for m >= 4. */ - if (m & 4) - { - IFLOAT *BO = B; - v4sf_t *rowC; - v4sf_t result[4]; - __vector_quad acc0; - __builtin_mma_xxsetaccz (&acc0); - BLASLONG l = 0; - for (l = 0; l < k / 2; l++) - { - vector short rowB = - { BO[l << 1], BO[(l << 1) + 1], 0, 0, 0, 0, 0, 0}; - vec_t *rowA = (vec_t *) & (AO[l << 3]); - MMA (&acc0, (vec_t) rowB, rowA[0]); - } - if (k % 2 == 1) - { - if (k > 1) - l = (k / 2) << 1; - vector short rowB = { BO[l], 0, 0, 0, 0, 0, 0, 0 }; - vector short rowA = - { AO[(l << 2)], 0, AO[(l << 2) + 1] , 0 , - AO[(l << 2) + 2], 0, AO[(l << 2) + 3], 0 }; - MMA (&acc0, (vec_t) rowB, (vec_t)(rowA)); - } - rowC = (v4sf_t *) &CO[0]; - __builtin_mma_disassemble_acc ((void *)result, &acc0); - rowC[0] += result[0] * alpha; - AO += k << 2; - BO += k; - CO += 4; - } - /* Loop for m >= 2. */ - if (m & 2) - { - IFLOAT *BO = B; - BLASLONG l = 0; - v4sf_t t = { 0, 0, 0, 0 }; - for (l = 0; l < k; l++) - { - v4sf_t rowB = { BF16TOF32 (BO[l]), BF16TOF32 (BO[l]), 0, 0 }; - v4sf_t rowA = - { BF16TOF32 (AO[l << 1]), BF16TOF32 (AO[(l << 1) + 1]), 0, - 0 - }; - t += rowA * rowB; - } - t = t * valpha; - CO[0] += t[0]; - CO[1] += t[1]; - AO += k << 1; - BO += k; - CO += 2; - } - /* Loop for m = 1. */ - if (m & 1) - { - IFLOAT *BO = B; - BLASLONG l = 0; - FLOAT t = 0; - for (l = 0; l < k; l++) - { - t += BF16TOF32 (AO[l]) * BF16TOF32 (BO[l]); - } - AO += k; - BO += k; - CO[0] += t * alpha; - CO += 1; - } - B += k; + if (n & 1) { + BLASLONG j; + FLOAT *CO; + IFLOAT *AO; + CO = C; + C += ldc; + AO = A; + /* Loop for m >= 16. */ + for (j = 0; j < (m >> 4); j++) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + __vector_quad acc0, acc1, acc2, acc3; + SET_ACC_ZERO4 (); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vector short rowB = { BO[l << 1], BO[(l << 1) + 1], 0, 0, 0, 0, 0, 0}; + vec_t *rowA = (vec_t *) & (AO[l << 5]); + MMA (&acc0, (vec_t) rowB, rowA[0]); + MMA (&acc1, (vec_t) rowB, rowA[1]); + MMA (&acc2, (vec_t) rowB, rowA[2]); + MMA (&acc3, (vec_t) rowB, rowA[3]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 1; + vector short rowB = { BO[l], 0, 0, 0, 0, 0, 0, 0 }; + vec_t *rowA = (vec_t *) & (AO[(l << 4)]); + MMA (&acc0, (vec_t) rowB, MERGE_HIGH (rowA[0], vzero)); + MMA (&acc1, (vec_t) rowB, MERGE_LOW (rowA[0], vzero)); + MMA (&acc2, (vec_t) rowB, MERGE_HIGH (rowA[1], vzero)); + MMA (&acc3, (vec_t) rowB, MERGE_LOW (rowA[1], vzero)); + } +#ifdef BGEMM + __builtin_mma_disassemble_acc ((void *)result, &acc0); + STORE4_BF16 (&CO[0], result[0], falpha); + __builtin_mma_disassemble_acc ((void *)result, &acc1); + STORE4_BF16 (&CO[4], result[0], falpha); + __builtin_mma_disassemble_acc ((void *)result, &acc2); + STORE4_BF16 (&CO[8], result[0], falpha); + __builtin_mma_disassemble_acc ((void *)result, &acc3); + STORE4_BF16 (&CO[12], result[0], falpha); +#else + rowC = (v4sf_t *) &CO[0]; + __builtin_mma_disassemble_acc ((void *)result, &acc0); + rowC[0] += result[0] * alpha; + __builtin_mma_disassemble_acc ((void *)result, &acc1); + rowC[1] += result[0] * alpha; + __builtin_mma_disassemble_acc ((void *)result, &acc2); + rowC[2] += result[0] * alpha; + __builtin_mma_disassemble_acc ((void *)result, &acc3); + rowC[3] += result[0] * alpha; +#endif + AO += k << 4; + BO += k; + CO += 16; + } + /* Loop for m >= 8. */ + if (m & 8) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + __vector_quad acc0, acc1; + __builtin_mma_xxsetaccz (&acc0); + __builtin_mma_xxsetaccz (&acc1); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vector short rowB = { BO[l << 1], BO[(l << 1) + 1], 0, 0, 0, 0, 0, 0}; + vec_t *rowA = (vec_t *) & (AO[l << 4]); + MMA (&acc0, (vec_t) rowB, rowA[0]); + MMA (&acc1, (vec_t) rowB, rowA[1]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 1; + vector short rowB = { BO[l], 0, 0, 0, 0, 0, 0, 0 }; + vec_t *rowA = (vec_t *) & (AO[(l << 3)]); + MMA (&acc0, (vec_t) rowB, MERGE_HIGH (rowA[0], vzero)); + MMA (&acc1, (vec_t) rowB, MERGE_LOW (rowA[0], vzero)); + } +#ifdef BGEMM + __builtin_mma_disassemble_acc ((void *)result, &acc0); + STORE4_BF16 (&CO[0], result[0], falpha); + __builtin_mma_disassemble_acc ((void *)result, &acc1); + STORE4_BF16 (&CO[4], result[0], falpha); +#else + rowC = (v4sf_t *) &CO[0]; + __builtin_mma_disassemble_acc ((void *)result, &acc0); + rowC[0] += result[0] * alpha; + __builtin_mma_disassemble_acc ((void *)result, &acc1); + rowC[1] += result[0] * alpha; +#endif + AO += k << 3; + BO += k; + CO += 8; + } + /* Loop for m >= 4. */ + if (m & 4) { + IFLOAT *BO = B; +#ifndef BGEMM + v4sf_t *rowC; +#endif + v4sf_t result[4]; + __vector_quad acc0; + __builtin_mma_xxsetaccz (&acc0); + BLASLONG l = 0; + for (l = 0; l < k / 2; l++) { + vector short rowB = { BO[l << 1], BO[(l << 1) + 1], 0, 0, 0, 0, 0, 0}; + vec_t *rowA = (vec_t *) & (AO[l << 3]); + MMA (&acc0, (vec_t) rowB, rowA[0]); + } + if (k % 2 == 1) { + if (k > 1) + l = (k / 2) << 1; + vector short rowB = { BO[l], 0, 0, 0, 0, 0, 0, 0 }; + vector short rowA = { AO[(l << 2)], 0, AO[(l << 2) + 1] , 0 , + AO[(l << 2) + 2], 0, AO[(l << 2) + 3], 0 }; + MMA (&acc0, (vec_t) rowB, (vec_t)(rowA)); + } +#ifdef BGEMM + __builtin_mma_disassemble_acc ((void *)result, &acc0); + STORE4_BF16 (&CO[0], result[0], falpha); +#else + rowC = (v4sf_t *) &CO[0]; + __builtin_mma_disassemble_acc ((void *)result, &acc0); + rowC[0] += result[0] * alpha; +#endif + AO += k << 2; + BO += k; + CO += 4; + } + /* Loop for m >= 2. */ + if (m & 2) { + IFLOAT *BO = B; + BLASLONG l = 0; + v4sf_t t = { 0, 0, 0, 0 }; + for (l = 0; l < k; l++) { + v4sf_t rowB = { BF16TOF32 (BO[l]), BF16TOF32 (BO[l]), 0, 0 }; + v4sf_t rowA = { BF16TOF32 (AO[l << 1]), BF16TOF32 (AO[(l << 1) + 1]), 0, 0 }; + t += rowA * rowB; + } + t = t * valpha; +#ifdef BGEMM + CO[0] = f32_to_bf16_scalar (BF16TOF32 (CO[0]) + t[0]); + CO[1] = f32_to_bf16_scalar (BF16TOF32 (CO[1]) + t[1]); +#else + CO[0] += t[0]; + CO[1] += t[1]; +#endif + AO += k << 1; + BO += k; + CO += 2; + } + /* Loop for m = 1. */ + if (m & 1) { + IFLOAT *BO = B; + BLASLONG l = 0; + float t = 0; /* float, not FLOAT: accumulate in fp32 regardless of BGEMM */ + for (l = 0; l < k; l++) { + t += BF16TOF32 (AO[l]) * BF16TOF32 (BO[l]); + } + AO += k; + BO += k; +#ifdef BGEMM + CO[0] = f32_to_bf16_scalar (BF16TOF32 (CO[0]) + (float)t * falpha); +#else + CO[0] += t * alpha; +#endif + CO += 1; + } + B += k; } - return 0; + return 0; } diff --git a/param.h b/param.h index 8e7ab29f24..bc179416f3 100644 --- a/param.h +++ b/param.h @@ -2654,6 +2654,20 @@ USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #define SBGEMM_DEFAULT_P 512 #define SBGEMM_DEFAULT_Q 1024 #define SBGEMM_DEFAULT_R 4096 + +/* BGEMM (BF16->BF16) on POWER10: tile matches the 16x8 primary tile in + * sbgemm_kernel_power10.c. A-panel packed by sbgemm_ncopy_16, B-panel by + * sbgemm_ncopy_8 (identical BF16 pair interleaving, reused for BGEMM). */ +#undef BGEMM_DEFAULT_UNROLL_M +#undef BGEMM_DEFAULT_UNROLL_N +#undef BGEMM_DEFAULT_P +#undef BGEMM_DEFAULT_Q +#undef BGEMM_DEFAULT_R +#define BGEMM_DEFAULT_UNROLL_M 16 +#define BGEMM_DEFAULT_UNROLL_N 8 +#define BGEMM_DEFAULT_P 512 +#define BGEMM_DEFAULT_Q 1024 +#define BGEMM_DEFAULT_R 4096 #endif #if defined(SPARC) && defined(V7) diff --git a/utest/Makefile b/utest/Makefile index c3bf7e706b..897559304f 100644 --- a/utest/Makefile +++ b/utest/Makefile @@ -30,6 +30,7 @@ OBJS_EXT+=$(DIR_EXT)/test_cgemv_t.o $(DIR_EXT)/test_zgemv_t.o $(DIR_EXT)/test_cg OBJS_EXT+=$(DIR_EXT)/test_sgemmt.o $(DIR_EXT)/test_dgemmt.o $(DIR_EXT)/test_cgemmt.o $(DIR_EXT)/test_zgemmt.o OBJS_EXT+=$(DIR_EXT)/test_ztrmv.o $(DIR_EXT)/test_ctrmv.o $(DIR_EXT)/test_ztrsv.o $(DIR_EXT)/test_ctrsv.o OBJS_EXT+=$(DIR_EXT)/test_zgemm.o $(DIR_EXT)/test_cgemm.o $(DIR_EXT)/test_zgbmv.o $(DIR_EXT)/test_cgbmv.o +OBJS_EXT+=$(DIR_EXT)/test_bgemm.o ifneq ($(NO_LAPACK), 1) OBJS += test_potrs.o diff --git a/utest/test_extensions/test_bgemm.c b/utest/test_extensions/test_bgemm.c new file mode 100644 index 0000000000..9e76b4a246 --- /dev/null +++ b/utest/test_extensions/test_bgemm.c @@ -0,0 +1,372 @@ +/***************************************************************************** +Copyright (c) 2025, 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 COPYRIGHT OWNER 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. +*****************************************************************************/ + +/** + * Unit tests for BGEMM: BFloat16-in, BFloat16-out GEMM. + * + * Strategy: compute the same operation via SBGEMM (BF16-in, float32-out) as + * a trusted reference, then widen the BGEMM BF16 output back to float32 and + * compare element-wise. Because both kernels use the same MMA accumulation + * the only difference is the store path; relative tolerance is set to allow + * for the one extra BF16 rounding on the output. + * + * BF16 helper: bfloat16_to_float() converts a stored bfloat16 value to + * float32 for comparison without depending on the BFLOAT16CONVERSION macro. + */ + +#include "utest/openblas_utest.h" +#include "common.h" +#include +#include +#include + +#ifdef BUILD_BFLOAT16 + +#define DATASIZE 100 + +/* ----------------------------------------------------------------------- + * BF16 helpers + * --------------------------------------------------------------------- */ + +/* Widen a bfloat16 bit-pattern to float32 by zero-extending the low 16 bits. + * This matches the OpenBLAS BF16TOF32 macro on both endiannesses because the + * bfloat16 typedef stores the *high* 16 bits of a float32 regardless of host + * endianness. */ +static inline float bf16_to_float(bfloat16 b) +{ + float f = 0.0f; + unsigned short *pf = (unsigned short *)&f; +#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ + pf[0] = (unsigned short)b; +#else + pf[1] = (unsigned short)b; +#endif + return f; +} + +/* Truncate float32 to bfloat16 (round-to-nearest, ties-to-even omitted for + * simplicity — matching what sbstobf16_ produces from a float already + * rounded to BF16 representable values in the test data). */ +static inline bfloat16 float_to_bf16(float f) +{ + unsigned short *p = (unsigned short *)&f; +#if __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ + return (bfloat16)p[0]; +#else + return (bfloat16)p[1]; +#endif +} + +/* ----------------------------------------------------------------------- + * Test data (static to avoid stack pressure for DATASIZE=100) + * --------------------------------------------------------------------- */ + +struct DATA_BGEMM { + bfloat16 a_bf16[DATASIZE * DATASIZE]; /* BF16 input A */ + bfloat16 b_bf16[DATASIZE * DATASIZE]; /* BF16 input B */ + bfloat16 c_bgemm[DATASIZE * DATASIZE]; /* BGEMM BF16 output */ + float c_sbgemm[DATASIZE * DATASIZE]; /* SBGEMM float32 reference */ + bfloat16 c_init[DATASIZE * DATASIZE]; /* saved initial C (BF16) */ + float c_init_f[DATASIZE * DATASIZE]; /* saved initial C (float32) */ +}; + +static struct DATA_BGEMM data_bgemm; + +/* ----------------------------------------------------------------------- + * Core check helper + * + * 1. Fill A, B, C with random BF16 values. + * 2. Run SBGEMM (trusted, float32 C) with float32 alpha/beta. + * 3. Run BGEMM (tested, BF16 C) with BF16 alpha/beta. + * 4. Widen BGEMM output to float32 and diff against SBGEMM reference. + * Returns the max absolute difference (normalised by |ref| + 1). + * --------------------------------------------------------------------- */ +static float check_bgemm(char transa, char transb, + blasint m, blasint n, blasint k, + float alpha_f, float beta_f, + blasint lda, blasint ldb, blasint ldc) +{ + blasint i, j; + float max_rdiff = 0.0f; + + /* Convert scalar parameters to BF16 for BGEMM */ + bfloat16 alpha_bf16 = float_to_bf16(alpha_f); + bfloat16 beta_bf16 = float_to_bf16(beta_f); + + /* Widen back to float so both kernels use the same representable value */ + float alpha_rep = bf16_to_float(alpha_bf16); + float beta_rep = bf16_to_float(beta_bf16); + + /* Determine A and B sizes based on transpose flags */ + blasint a_rows = (transa == 'N' || transa == 'n') ? m : k; + blasint a_cols = (transa == 'N' || transa == 'n') ? k : m; + blasint b_rows = (transb == 'N' || transb == 'n') ? k : n; + blasint b_cols = (transb == 'N' || transb == 'n') ? n : k; + (void)b_rows; /* unused, kept for clarity */ + (void)a_rows; + + /* Fill A and B with random BF16 values in [0.5, 1.5] */ + for (i = 0; i < a_cols * lda; i++) + data_bgemm.a_bf16[i] = float_to_bf16(0.5f + (float)rand() / (float)RAND_MAX); + for (i = 0; i < b_cols * ldb; i++) + data_bgemm.b_bf16[i] = float_to_bf16(0.5f + (float)rand() / (float)RAND_MAX); + + /* Fill C with random BF16 values, saving copies for both kernels */ + for (i = 0; i < n * ldc; i++) { + data_bgemm.c_init[i] = float_to_bf16(0.5f + (float)rand() / (float)RAND_MAX); + data_bgemm.c_init_f[i] = bf16_to_float(data_bgemm.c_init[i]); + } + + /* Copy initial C into both output buffers */ + memcpy(data_bgemm.c_bgemm, data_bgemm.c_init, n * ldc * sizeof(bfloat16)); + memcpy(data_bgemm.c_sbgemm, data_bgemm.c_init_f, n * ldc * sizeof(float)); + + /* Reference: SBGEMM (BF16 A/B, float32 C) */ + BLASFUNC(sbgemm)(&transa, &transb, &m, &n, &k, + &alpha_rep, + data_bgemm.a_bf16, &lda, + data_bgemm.b_bf16, &ldb, + &beta_rep, + data_bgemm.c_sbgemm, &ldc); + + /* Tested: BGEMM (BF16 A/B/C/alpha/beta) */ + BLASFUNC(bgemm)(&transa, &transb, &m, &n, &k, + &alpha_bf16, + data_bgemm.a_bf16, &lda, + data_bgemm.b_bf16, &ldb, + &beta_bf16, + data_bgemm.c_bgemm, &ldc); + + /* Compare: widen BGEMM BF16 output to float32, diff vs SBGEMM float32 */ + for (j = 0; j < n; j++) { + for (i = 0; i < m; i++) { + float ref = data_bgemm.c_sbgemm[j * ldc + i]; + float test = bf16_to_float(data_bgemm.c_bgemm[j * ldc + i]); + float rdiff = fabsf(test - ref) / (fabsf(ref) + 1.0f); + if (rdiff > max_rdiff) + max_rdiff = rdiff; + } + } + return max_rdiff; +} + +/* + * Tolerance: SBGEMM accumulates in float32, stores float32. + * BGEMM does the same accumulation then converts float32 → BF16 on store, + * introducing at most 1 BF16 ULP (≈ 2^-7 ≈ 0.0078). Use 0.01 to match + * the tolerance used in compare_sgemm_sbgemm.c. + */ +#define BGEMM_TOL 0.01f + +/* ======================================================================= + * Test cases: NN, TN, NT, TT × a few representative sizes + edge cases + * ===================================================================== */ + +/* ----- NN ----- */ +CTEST(bgemm, nn_M50_N50_K50) +{ + float rdiff = check_bgemm('N', 'N', 50, 50, 50, + 1.5f, 0.0f, 50, 50, 50); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +CTEST(bgemm, nn_M100_N100_K100) +{ + float rdiff = check_bgemm('N', 'N', 100, 100, 100, + 1.0f, 0.0f, 100, 100, 100); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* Hits the 16x8 primary tile exactly (kernel UNROLL_M=16, UNROLL_N=8) */ +CTEST(bgemm, nn_M16_N8_K16) +{ + float rdiff = check_bgemm('N', 'N', 16, 8, 16, + 1.0f, 0.0f, 16, 16, 16); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* m and n not multiples of tile — exercises all remainder paths */ +CTEST(bgemm, nn_remainder_M17_N9_K7) +{ + float rdiff = check_bgemm('N', 'N', 17, 9, 7, + 1.0f, 0.0f, 17, 7, 17); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* m=1 and n=1 — exercises the scalar (m&1) and (n&1) remainder paths */ +CTEST(bgemm, nn_M1_N1_K8) +{ + float rdiff = check_bgemm('N', 'N', 1, 1, 8, + 1.0f, 0.0f, 1, 8, 1); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* ----- TN ----- */ +CTEST(bgemm, tn_M50_N50_K50) +{ + float rdiff = check_bgemm('T', 'N', 50, 50, 50, + 1.5f, 0.0f, 50, 50, 50); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +CTEST(bgemm, tn_M100_N50_K50) +{ + float rdiff = check_bgemm('T', 'N', 100, 50, 50, + 1.0f, 0.0f, 50, 50, 100); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* ----- NT ----- */ +CTEST(bgemm, nt_M50_N50_K100) +{ + float rdiff = check_bgemm('N', 'T', 50, 50, 100, + 1.0f, 0.0f, 50, 50, 50); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* ----- TT ----- */ +CTEST(bgemm, tt_M50_N50_K50) +{ + float rdiff = check_bgemm('T', 'T', 50, 50, 50, + 1.5f, 2.0f, 50, 50, 50); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* ----- alpha = 0: C should equal beta * C_init ----- */ +CTEST(bgemm, alpha_zero_nn) +{ + float rdiff = check_bgemm('N', 'N', 50, 50, 50, + 0.0f, 2.0f, 50, 50, 50); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* ----- beta = 0: C should be fresh accumulation, no add from old C ----- */ +CTEST(bgemm, beta_zero_nn) +{ + float rdiff = check_bgemm('N', 'N', 50, 50, 50, + 1.0f, 0.0f, 50, 50, 50); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* ----- beta = 1: accumulate into C, no scaling on existing C ----- */ +CTEST(bgemm, beta_one_nn) +{ + float rdiff = check_bgemm('N', 'N', 50, 50, 50, + 2.0f, 1.0f, 50, 50, 50); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* ----- odd k: exercises the k%2==1 remainder path in the kernel ----- */ +CTEST(bgemm, nn_odd_k_M32_N8_K7) +{ + float rdiff = check_bgemm('N', 'N', 32, 8, 7, + 1.0f, 0.0f, 32, 7, 32); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +CTEST(bgemm, nn_odd_k_M16_N8_K1) +{ + float rdiff = check_bgemm('N', 'N', 16, 8, 1, + 1.0f, 0.0f, 16, 1, 16); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* ----- Large n>8 sweep (n&4, n&2, n&1 paths) ----- */ +CTEST(bgemm, nn_n4_path_M32_N4_K16) +{ + float rdiff = check_bgemm('N', 'N', 32, 4, 16, + 1.0f, 0.0f, 32, 16, 32); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +CTEST(bgemm, nn_n2_path_M32_N2_K16) +{ + float rdiff = check_bgemm('N', 'N', 32, 2, 16, + 1.0f, 0.0f, 32, 16, 32); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +CTEST(bgemm, nn_n1_path_M32_N1_K16) +{ + float rdiff = check_bgemm('N', 'N', 32, 1, 16, + 1.0f, 0.0f, 32, 16, 32); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* ----- m remainder paths: m&8, m&4, m&2, m&1 within the n>=8 loop ----- */ +CTEST(bgemm, nn_m8_path_M8_N8_K16) +{ + float rdiff = check_bgemm('N', 'N', 8, 8, 16, + 1.0f, 0.0f, 8, 16, 8); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +CTEST(bgemm, nn_m4_path_M4_N8_K16) +{ + float rdiff = check_bgemm('N', 'N', 4, 8, 16, + 1.0f, 0.0f, 4, 16, 4); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +CTEST(bgemm, nn_m2_path_M2_N8_K16) +{ + float rdiff = check_bgemm('N', 'N', 2, 8, 16, + 1.0f, 0.0f, 2, 16, 2); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +CTEST(bgemm, nn_m1_path_M1_N8_K16) +{ + float rdiff = check_bgemm('N', 'N', 1, 8, 16, + 1.0f, 0.0f, 1, 16, 1); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* ----- m=32 sweep hits the m>=32 fast path (2x16 tile) ----- */ +CTEST(bgemm, nn_m32_N8_K16) +{ + float rdiff = check_bgemm('N', 'N', 32, 8, 16, + 1.0f, 0.0f, 32, 16, 32); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +/* ----- alpha and beta both BF16-representable non-trivial values ----- */ +CTEST(bgemm, nn_alpha1p5_beta2_M50_N50_K50) +{ + /* alpha=1.5, beta=2.0 are exactly BF16-representable */ + float rdiff = check_bgemm('N', 'N', 50, 50, 50, + 1.5f, 2.0f, 50, 50, 50); + ASSERT_DBL_NEAR_TOL(0.0f, rdiff, BGEMM_TOL); +} + +#endif /* BUILD_BFLOAT16 */