hey! this is a static git repo viewer. commit diffs, file history, older commits and other revisions won't load without javascript.

you are still able to view source files on branch/tag tips, but the experience is degraded.

char-fork/llama.cpp

git clone https://git.t4t.associates/char-fork/llama.cpp

shalinib-ibmggml-cpu: Enable BF16 tiled gemm optimization on PowerPC (#26068)b910200892026-07-27

1// Copyright 2024 Mozilla Foundation
2//
3// Permission is hereby granted, free of charge, to any person obtaining
4// a copy of this software and associated documentation files (the
5// "Software"), to deal in the Software without restriction, including
6// without limitation the rights to use, copy, modify, merge, publish,
7// distribute, sublicense, and/or sell copies of the Software, and to
8// permit persons to whom the Software is furnished to do so, subject to
9// the following conditions:
10//
11// The above copyright notice and this permission notice shall be
12// included in all copies or substantial portions of the Software.
13//
14// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
17// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
18// BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
19// ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
20// CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
21// SOFTWARE.
22
23//
24//                   _   _          ___ _      _   ___
25//                  | |_(_)_ _ _  _| _ ) |    /_\ / __|
26//                  |  _| | ' \ || | _ \ |__ / _ \\__ \.
27//                   \__|_|_||_\_, |___/____/_/ \_\___/
28//                             |__/
29//
30//                    BASIC LINEAR ALGEBRA SUBPROGRAMS
31//
32//
33// This file implements multithreaded CPU matrix multiplication for the
34// common contiguous use case C = Aᵀ * B. These kernels are designed to
35// have excellent performance[1] for matrices that fit in the CPU cache
36// without imposing any overhead such as cache filling or malloc calls.
37//
38// This implementation does not guarantee any upper bound with rounding
39// errors, which grow along with k. Our goal's to maximally exploit the
40// hardware for performance, and then use whatever resources remain for
41// improving numerical accuracy.
42//
43// [1] J. Tunney, ‘LLaMA Now Goes Faster on CPUs’, Mar. 2024. [Online].
44//     Available: https://justine.lol/matmul/. [Accessed: 29-Mar-2024].
45
46#if defined(__GNUC__)
47#pragma GCC diagnostic ignored "-Wpedantic"
48#pragma GCC diagnostic ignored "-Wignored-attributes"
49#endif
50
51#include "sgemm.h"
52#include "ggml-impl.h"
53#include "ggml-cpu-impl.h"
54#include "ggml-quants.h"
55#include "simd-mappings.h"
56
57#include <array>
58#include <type_traits>
59
60#ifdef _MSC_VER
61#define NOINLINE __declspec(noinline)
62#else
63#define NOINLINE __attribute__((__noinline__))
64#endif
65
66#if defined(__ARM_NEON) || defined(__AVX512F__) || defined(__VXE__) || defined(__VXE2__)
67#define VECTOR_REGISTERS 32
68#else
69#define VECTOR_REGISTERS 16
70#endif
71
72#if defined(__riscv_v_intrinsic)
73#define LMUL 4
74#endif
75
76#define MM256_SET_M128I(a, b) _mm256_insertf128_si256(_mm256_castsi128_si256(b), (a), 1)
77
78namespace {
79
80inline float unhalf(ggml_fp16_t d) {
81    return GGML_CPU_FP16_TO_FP32(d);
82}
83
84////////////////////////////////////////////////////////////////////////////////////////////////////
85// VECTORIZED ARITHMETIC OPERATIONS
86
87#if defined(__SSE__) || defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__)
88inline __m128 add(__m128 x, __m128 y) { return _mm_add_ps(x, y); }
89inline __m128 sub(__m128 x, __m128 y) { return _mm_sub_ps(x, y); }
90inline __m128 mul(__m128 x, __m128 y) { return _mm_mul_ps(x, y); }
91#endif  // __SSE__
92
93#if defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__)
94inline __m256 add(__m256 x, __m256 y) { return _mm256_add_ps(x, y); }
95inline __m256 sub(__m256 x, __m256 y) { return _mm256_sub_ps(x, y); }
96inline __m256 mul(__m256 x, __m256 y) { return _mm256_mul_ps(x, y); }
97#endif // __AVX__
98
99#if defined(__AVX512F__)
100inline __m512 add(__m512 x, __m512 y) { return _mm512_add_ps(x, y); }
101inline __m512 sub(__m512 x, __m512 y) { return _mm512_sub_ps(x, y); }
102inline __m512 mul(__m512 x, __m512 y) { return _mm512_mul_ps(x, y); }
103#endif // __AVX512F__
104
105#if defined(__ARM_NEON)
106inline float32x4_t add(float32x4_t x, float32x4_t y) { return vaddq_f32(x, y); }
107inline float32x4_t sub(float32x4_t x, float32x4_t y) { return vsubq_f32(x, y); }
108inline float32x4_t mul(float32x4_t x, float32x4_t y) { return vmulq_f32(x, y); }
109#endif // __ARM_NEON
110
111#if defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC)
112inline float16x8_t add(float16x8_t x, float16x8_t y) { return vaddq_f16(x, y); }
113inline float16x8_t sub(float16x8_t x, float16x8_t y) { return vsubq_f16(x, y); }
114inline float16x8_t mul(float16x8_t x, float16x8_t y) { return vmulq_f16(x, y); }
115#endif // __ARM_FEATURE_FP16_VECTOR_ARITHMETIC
116
117#if defined(__VXE__) || defined(__VXE2__)
118inline float32x4_t add(float32x4_t x, float32x4_t y) { return vec_add(x, y); }
119inline float32x4_t sub(float32x4_t x, float32x4_t y) { return vec_sub(x, y); }
120inline float32x4_t mul(float32x4_t x, float32x4_t y) { return vec_mul(x, y); }
121#endif
122
123#if defined(__MMA__)
124typedef vector unsigned char vec_t;
125typedef __vector_quad acc_t;
126#endif
127////////////////////////////////////////////////////////////////////////////////////////////////////
128// VECTORIZED FUSED MULTIPLY ADD
129
130/**
131 * Computes a * b + c.
132 */
133template <typename T, typename U>
134inline U madd(T a, T b, U c) {
135    return add(mul(a, b), c);
136}
137
138#if defined(__FMA__)
139#if defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__)
140template <>
141inline __m256 madd(__m256 a, __m256 b, __m256 c) {
142    return _mm256_fmadd_ps(a, b, c);
143}
144#endif
145#if defined(__AVX512F__)
146template <>
147inline __m512 madd(__m512 a, __m512 b, __m512 c) {
148    return _mm512_fmadd_ps(a, b, c);
149}
150#endif
151#if defined(__AVX512BF16__)
152template <>
153inline __m512 madd(__m512bh a, __m512bh b, __m512 c) {
154    return _mm512_dpbf16_ps(c, a, b);
155}
156template <>
157inline __m256 madd(__m256bh a, __m256bh b, __m256 c) {
158    return _mm256_dpbf16_ps(c, a, b);
159}
160#endif
161#endif
162
163#if defined(__ARM_FEATURE_FMA)
164template <>
165inline float32x4_t madd(float32x4_t a, float32x4_t b, float32x4_t c) {
166    return vfmaq_f32(c, b, a);
167}
168#if defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC) && !defined(_MSC_VER)
169template <>
170inline float16x8_t madd(float16x8_t a, float16x8_t b, float16x8_t c) {
171    return vfmaq_f16(c, b, a);
172}
173#endif
174#endif
175
176#if defined(__VXE__) || defined(__VXE2__)
177template <>
178inline float32x4_t madd(float32x4_t a, float32x4_t b, float32x4_t c) {
179    return vec_madd(a, b, c);
180}
181#endif
182
183#if defined(__riscv_v_intrinsic)
184template <> inline vfloat32m1_t madd(vfloat32m1_t a, vfloat32m1_t b, vfloat32m1_t c) {
185    return __riscv_vfmacc_vv_f32m1(c, a, b, __riscv_vsetvlmax_e32m1());
186}
187template <> inline vfloat32m2_t madd(vfloat32m2_t a, vfloat32m2_t b, vfloat32m2_t c) {
188    return __riscv_vfmacc_vv_f32m2(c, a, b, __riscv_vsetvlmax_e32m2());
189}
190template <> inline vfloat32m4_t madd(vfloat32m4_t a, vfloat32m4_t b, vfloat32m4_t c) {
191    return __riscv_vfmacc_vv_f32m4(c, a, b, __riscv_vsetvlmax_e32m4());
192}
193template <> inline vfloat32m8_t madd(vfloat32m8_t a, vfloat32m8_t b, vfloat32m8_t c) {
194    return __riscv_vfmacc_vv_f32m8(c, a, b, __riscv_vsetvlmax_e32m8());
195}
196#endif
197
198#if defined(__riscv_zvfh)
199template <> inline vfloat32m1_t madd(vfloat16mf2_t a, vfloat16mf2_t b, vfloat32m1_t c) {
200    return __riscv_vfwmacc_vv_f32m1(c, a, b, __riscv_vsetvlmax_e32m1());
201}
202template <> inline vfloat32m2_t madd(vfloat16m1_t a, vfloat16m1_t b, vfloat32m2_t c) {
203    return __riscv_vfwmacc_vv_f32m2(c, a, b, __riscv_vsetvlmax_e32m2());
204}
205template <> inline vfloat32m4_t madd(vfloat16m2_t a, vfloat16m2_t b, vfloat32m4_t c) {
206    return __riscv_vfwmacc_vv_f32m4(c, a, b, __riscv_vsetvlmax_e32m4());
207}
208template <> inline vfloat32m8_t madd(vfloat16m4_t a, vfloat16m4_t b, vfloat32m8_t c) {
209    return __riscv_vfwmacc_vv_f32m8(c, a, b, __riscv_vsetvlmax_e32m8());
210}
211#endif
212
213#if defined(__riscv_zvfbfwma)
214template <> inline vfloat32m1_t madd(vbfloat16mf2_t a, vbfloat16mf2_t b, vfloat32m1_t c) {
215    return __riscv_vfwmaccbf16_vv_f32m1(c, a, b, __riscv_vsetvlmax_e32m1());
216}
217template <> inline vfloat32m2_t madd(vbfloat16m1_t a, vbfloat16m1_t b, vfloat32m2_t c) {
218    return __riscv_vfwmaccbf16_vv_f32m2(c, a, b, __riscv_vsetvlmax_e32m2());
219}
220template <> inline vfloat32m4_t madd(vbfloat16m2_t a, vbfloat16m2_t b, vfloat32m4_t c) {
221    return __riscv_vfwmaccbf16_vv_f32m4(c, a, b, __riscv_vsetvlmax_e32m4());
222}
223template <> inline vfloat32m8_t madd(vbfloat16m4_t a, vbfloat16m4_t b, vfloat32m8_t c) {
224    return __riscv_vfwmaccbf16_vv_f32m8(c, a, b, __riscv_vsetvlmax_e32m8());
225}
226#endif
227
228////////////////////////////////////////////////////////////////////////////////////////////////////
229// VECTORIZED HORIZONTAL SUM
230
231#if defined(__ARM_NEON)
232inline float hsum(float32x4_t x) {
233    return vaddvq_f32(x);
234}
235#endif // __ARM_NEON
236
237#if defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC) && !defined(_MSC_VER)
238inline float hsum(float16x8_t x) {
239    return vaddvq_f32(vaddq_f32(vcvt_f32_f16(vget_low_f16(x)),
240                                vcvt_f32_f16(vget_high_f16(x))));
241}
242#endif // __ARM_FEATURE_FP16_VECTOR_ARITHMETIC
243
244#if defined(__VXE__) || defined(__VXE2__)
245inline float hsum(float32x4_t x) {
246    float32x4_t tmp = x + vec_reve(x);
247    return tmp[0] + tmp[1];
248}
249#endif
250
251#if defined(__SSE__) || defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__)
252inline float hsum(__m128 x) {
253#if defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__)
254    x = _mm_add_ps(x, _mm_movehl_ps(x, x));
255    x = _mm_add_ss(x, _mm_movehdup_ps(x));
256#else
257    __m128 t;
258    t = _mm_shuffle_ps(x, x, _MM_SHUFFLE(2, 3, 0, 1));
259    x = _mm_add_ps(x, t);
260    t = _mm_movehl_ps(t, x);
261    x = _mm_add_ss(x, t);
262#endif
263    return _mm_cvtss_f32(x);
264}
265#endif
266
267#if defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__)
268inline float hsum(__m256 x) {
269    return hsum(_mm_add_ps(_mm256_extractf128_ps(x, 1),
270                           _mm256_castps256_ps128(x)));
271}
272#endif // __AVX__
273
274#if defined(__AVX512F__)
275inline float hsum(__m512 x) {
276    return _mm512_reduce_add_ps(x);
277}
278#endif // __AVX512F__
279
280#if defined(__riscv_v_intrinsic)
281inline float hsum(vfloat32m1_t x) {
282    return __riscv_vfmv_f_s_f32m1_f32(
283        __riscv_vfredusum_vs_f32m1_f32m1(x, __riscv_vfmv_v_f_f32m1(0, 1), __riscv_vsetvlmax_e32m1()));
284}
285inline float hsum(vfloat32m2_t x) {
286    return __riscv_vfmv_f_s_f32m1_f32(
287        __riscv_vfredusum_vs_f32m2_f32m1(x, __riscv_vfmv_v_f_f32m1(0, 1), __riscv_vsetvlmax_e32m2()));
288}
289inline float hsum(vfloat32m4_t x) {
290    return __riscv_vfmv_f_s_f32m1_f32(
291        __riscv_vfredusum_vs_f32m4_f32m1(x, __riscv_vfmv_v_f_f32m1(0, 1), __riscv_vsetvlmax_e32m4()));
292}
293inline float hsum(vfloat32m8_t x) {
294    return __riscv_vfmv_f_s_f32m1_f32(
295        __riscv_vfredusum_vs_f32m8_f32m1(x, __riscv_vfmv_v_f_f32m1(0, 1), __riscv_vsetvlmax_e32m8()));
296}
297#endif
298
299////////////////////////////////////////////////////////////////////////////////////////////////////
300// VECTORIZED MEMORY LOADING
301
302template <typename T, typename U> T load(const U *);
303
304#if defined(__ARM_NEON)
305template <> inline float32x4_t load(const float *p) {
306    return vld1q_f32(p);
307}
308#if !defined(_MSC_VER)
309// FIXME: this should check for __ARM_FEATURE_FP16_VECTOR_ARITHMETIC
310template <> inline float16x8_t load(const ggml_fp16_t *p) {
311    return vld1q_f16((const float16_t *)p);
312}
313template <> inline float32x4_t load(const ggml_fp16_t *p) {
314    return vcvt_f32_f16(vld1_f16((const float16_t *)p));
315}
316#endif // _MSC_VER
317#endif // __ARM_NEON
318
319#if defined(__VXE__) || defined(__VXE2__)
320template <> inline float32x4_t load(const ggml_fp16_t * p) {
321    float tmp[4];
322
323    for (int i = 0; i < 4; i++) {
324        tmp[i] = GGML_CPU_FP16_TO_FP32(p[i]);
325    }
326
327    return vec_xl(0, (const float *)(tmp));
328}
329template <> inline float32x4_t load(const float * p) {
330    return vec_xl(0, p);
331}
332#endif
333
334#if defined(__SSE__) || defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__)
335template <> inline __m128 load(const float *p) {
336    return _mm_loadu_ps(p);
337}
338#endif  // __SSE__
339
340#if defined(__AVX__) || defined(__AVX2__) || defined(__AVX512F__)
341template <> inline __m256 load(const float *p) {
342    return _mm256_loadu_ps(p);
343}
344#endif // __AVX__
345
346#if defined(__AVX2__) || defined(__AVX512F__)
347template <> inline __m256 load(const ggml_bf16_t *p) {
348    return _mm256_castsi256_ps(
349        _mm256_slli_epi32(_mm256_cvtepu16_epi32(_mm_loadu_si128((const __m128i *)p)), 16));
350}
351#endif // __AVX2__
352
353#if defined(__F16C__)
354template <> inline __m256 load(const ggml_fp16_t *p) {
355    return _mm256_cvtph_ps(_mm_loadu_si128((const __m128i *)p));
356}
357#endif // __F16C__
358
359#if defined(__AVX512F__)
360template <> inline __m512 load(const float *p) {
361    return _mm512_loadu_ps(p);
362}
363template <> inline __m512 load(const ggml_fp16_t *p) {
364    return _mm512_cvtph_ps(_mm256_loadu_si256((const __m256i *)p));
365}
366template <> inline __m512 load(const ggml_bf16_t *p) {
367    return _mm512_castsi512_ps(
368        _mm512_slli_epi32(_mm512_cvtepu16_epi32(_mm256_loadu_si256((const __m256i *)p)), 16));
369}
370#endif // __AVX512F__
371
372#if defined(__AVX512BF16__)
373template <> inline __m512bh load(const ggml_bf16_t *p) {
374    return (__m512bh)_mm512_loadu_ps((const float *)p);
375}
376template <> inline __m256bh load(const ggml_bf16_t *p) {
377    return (__m256bh)_mm256_loadu_ps((const float *)p);
378}
379template <> inline __m512bh load(const float *p) {
380    return _mm512_cvtne2ps_pbh(_mm512_loadu_ps(p + 16), _mm512_loadu_ps(p));
381}
382template <> inline __m256bh load(const float *p) {
383    return _mm512_cvtneps_pbh(_mm512_loadu_ps(p));
384}
385#endif
386
387#if defined(__riscv_v_intrinsic)
388template <> inline vfloat32m1_t load(const float *p) {
389    return __riscv_vle32_v_f32m1(p, __riscv_vsetvlmax_e32m1());
390}
391template <> inline vfloat32m2_t load(const float *p) {
392    return __riscv_vle32_v_f32m2(p, __riscv_vsetvlmax_e32m2());
393}
394template <> inline vfloat32m4_t load(const float *p) {
395    return __riscv_vle32_v_f32m4(p, __riscv_vsetvlmax_e32m4());
396}
397template <> inline vfloat32m8_t load(const float *p) {
398    return __riscv_vle32_v_f32m8(p, __riscv_vsetvlmax_e32m8());
399}
400#endif
401
402#if defined(__riscv_zvfh)
403template <> inline vfloat16mf2_t load(const ggml_fp16_t *p) {
404    return __riscv_vle16_v_f16mf2(reinterpret_cast<const _Float16 *>(p), __riscv_vsetvlmax_e16mf2());
405}
406template <> inline vfloat16m1_t load(const ggml_fp16_t *p) {
407    return __riscv_vle16_v_f16m1(reinterpret_cast<const _Float16 *>(p), __riscv_vsetvlmax_e16m1());
408}
409template <> inline vfloat16m2_t load(const ggml_fp16_t *p) {
410    return __riscv_vle16_v_f16m2(reinterpret_cast<const _Float16 *>(p), __riscv_vsetvlmax_e16m2());
411}
412template <> inline vfloat16m4_t load(const ggml_fp16_t *p) {
413    return __riscv_vle16_v_f16m4(reinterpret_cast<const _Float16 *>(p), __riscv_vsetvlmax_e16m4());
414}
415#endif
416
417#if defined(__riscv_zvfbfwma)
418template <> inline vbfloat16mf2_t load(const ggml_bf16_t *p) {
419    return __riscv_vle16_v_bf16mf2(reinterpret_cast<const __bf16*>(p), __riscv_vsetvlmax_e16mf2());
420}
421template <> inline vbfloat16m1_t load(const ggml_bf16_t *p) {
422    return __riscv_vle16_v_bf16m1(reinterpret_cast<const __bf16*>(p), __riscv_vsetvlmax_e16m1());
423}
424template <> inline vbfloat16m2_t load(const ggml_bf16_t *p) {
425    return __riscv_vle16_v_bf16m2(reinterpret_cast<const __bf16*>(p), __riscv_vsetvlmax_e16m2());
426}
427template <> inline vbfloat16m4_t load(const ggml_bf16_t *p) {
428    return __riscv_vle16_v_bf16m4(reinterpret_cast<const __bf16*>(p), __riscv_vsetvlmax_e16m4());
429}
430#endif
431
432#if defined(__riscv_v_intrinsic)
433template <typename T> T set_zero();
434
435template <> inline vfloat32m1_t set_zero() {
436    return __riscv_vfmv_v_f_f32m1(0.0f, __riscv_vsetvlmax_e32m1());
437}
438template <> inline vfloat32m2_t set_zero() {
439    return __riscv_vfmv_v_f_f32m2(0, __riscv_vsetvlmax_e32m2());
440}
441template <> inline vfloat32m4_t set_zero() {
442    return __riscv_vfmv_v_f_f32m4(0, __riscv_vsetvlmax_e32m4());
443}
444template <> inline vfloat32m8_t set_zero() {
445    return __riscv_vfmv_v_f_f32m8(0, __riscv_vsetvlmax_e32m8());
446}
447#endif
448
449#if defined(__riscv_v_intrinsic)
450template <typename T> size_t vlmax() {
451    if constexpr (std::is_same_v<T, vfloat32m1_t>) { return  __riscv_vsetvlmax_e32m1(); }
452    else if constexpr (std::is_same_v<T, vfloat32m2_t>) { return  __riscv_vsetvlmax_e32m2(); }
453    else if constexpr (std::is_same_v<T, vfloat32m4_t>) { return  __riscv_vsetvlmax_e32m4(); }
454    else if constexpr (std::is_same_v<T, vfloat32m8_t>) { return  __riscv_vsetvlmax_e32m8(); }
455    #if defined (__riscv_zvfh)
456    else if constexpr (std::is_same_v<T, vfloat16mf2_t>) { return  __riscv_vsetvlmax_e16mf2(); }
457    else if constexpr (std::is_same_v<T, vfloat16m1_t>) { return  __riscv_vsetvlmax_e16m1(); }
458    else if constexpr (std::is_same_v<T, vfloat16m2_t>) { return  __riscv_vsetvlmax_e16m2(); }
459    else if constexpr (std::is_same_v<T, vfloat16m4_t>) { return  __riscv_vsetvlmax_e16m4(); }
460    #endif
461    #if defined (__riscv_zvfbfwma)
462    else if constexpr (std::is_same_v<T, vbfloat16mf2_t>) { return  __riscv_vsetvlmax_e16mf2(); }
463    else if constexpr (std::is_same_v<T, vbfloat16m1_t>) { return  __riscv_vsetvlmax_e16m1(); }
464    else if constexpr (std::is_same_v<T, vbfloat16m2_t>) { return  __riscv_vsetvlmax_e16m2(); }
465    else if constexpr (std::is_same_v<T, vbfloat16m4_t>) { return  __riscv_vsetvlmax_e16m4(); }
466    #endif
467    return 0;
468}
469#endif
470
471////////////////////////////////////////////////////////////////////////////////////////////////////
472// FLOATING POINT MATRIX MULTIPLICATION
473
474template <int M>
475static inline int64_t BLOCK_SIZE(size_t m) {
476    const int64_t NB_BLOC_M = (m + M - 1) / M;
477    return (m % NB_BLOC_M == 0) ? m / NB_BLOC_M : (m / NB_BLOC_M) + 1;
478}
479
480static constexpr inline int64_t BLOC_POS(int64_t ib, int64_t ibN, int64_t bloc_size) {
481    return ib < ibN ? ib * bloc_size : ibN * bloc_size + (ib - ibN) * (bloc_size - 1);
482}
483
484template <int KN, typename D, typename V, typename TA, typename TB, typename TC>
485class tinyBLAS {
486  public:
487    tinyBLAS(const ggml_compute_params * params, int64_t k,
488             const TA *A, int64_t lda,
489             const TB *B, int64_t ldb,
490             TC *C, int64_t ldc)
491        : params(params), A(A), B(B), C(C), k(k), lda(lda), ldb(ldb), ldc(ldc) {
492    }
493
494    bool matmul(int64_t m, int64_t n) {
495        if (k % KN != 0)
496            return false;
497        // compute RM for only need tile with size RM&RM-1
498#if VECTOR_REGISTERS == 32
499        if (m % 16 == 0 && (m/16 >= params->nth)) {
500            const int64_t SIZE_N = BLOCK_SIZE<6>(n);
501            mnpack<4, 6, 4>(m, n, SIZE_N, 12);
502            return true;
503        }
504        if (m % 8 == 0 ) {
505            const int64_t SIZE_N = BLOCK_SIZE<6>(n);
506            mnpack<4, 6, 2>(m, n, SIZE_N, 12);
507            return true;
508        }
509        if (m % 4 == 0) {
510            const int64_t SIZE_N = BLOCK_SIZE<6>(n);
511            mnpack<4, 6, 1>(m, n, SIZE_N, 12);
512            return true;
513        }
514#else  // VECTOR_REGISTERS == 16
515        if (m % 16 == 0 && (m/16 >= params->nth)) {
516            const int64_t SIZE_N = BLOCK_SIZE<3>(n);
517            mnpack<4, 3, 4>(m, n, SIZE_N, 24);
518            return true;
519        }
520        if (m % 8 == 0 ) {
521            const int64_t SIZE_N = BLOCK_SIZE<3>(n);
522            mnpack<4, 3, 2>(m, n, SIZE_N, 24);
523            return true;
524        }
525        if (m % 4 == 0) {
526            const int64_t SIZE_N = BLOCK_SIZE<3>(n);
527            mnpack<4, 3, 1>(m, n, SIZE_N, 24);
528            return true;
529        }
530#endif
531        return false;
532    }
533
534  private:
535    template <int RM, int RN, int BM>
536    inline void mnpack(int64_t m, int64_t n, int64_t SIZE_N, int64_t BN) {
537        if (SIZE_N == RN) {
538            return gemm<RM, RN, BM>(m, n, BN);
539        }
540        if constexpr (RN > 1) {
541            return mnpack<RM, RN-1, BM>(m, n, SIZE_N, BN);
542        } else {
543            GGML_LOG_ERROR("mnpack<%d, %d> block size not supported\n", RM, (int)SIZE_N);
544            GGML_ASSERT(false); // we have miss something.
545        }
546    }
547
548    template <int RM, int RN>
549    inline void gemm_bloc(int64_t ii, int64_t jj) {
550        D Cv[RN][RM] = {};
551        for (int64_t l = 0; l < k; l += KN) {
552            // help compiler for op order.
553            if constexpr (RM <= RN) {
554                V Av[RM];
555                for (int64_t i = 0; i < RM; ++i) {
556                    Av[i] = load<V>(A + lda * (ii + i) + l);
557                }
558                for (int64_t j = 0; j < RN; ++j) {
559                    V Bv = load<V>(B + ldb * (jj + j) + l);
560                    for (int64_t i = 0; i < RM; ++i) {
561                        Cv[j][i] = madd(Av[i], Bv, Cv[j][i]);
562                    }
563                }
564            } else {
565                V Bv[RN];
566                for (int64_t j = 0; j < RN; ++j) {
567                    Bv[j] = load<V>(B + ldb * (jj + j) + l);
568                }
569                for (int64_t i = 0; i < RM; ++i) {
570                    V Av = load<V>(A + lda * (ii + i) + l);
571                    for (int64_t j = 0; j < RN; ++j) {
572                        Cv[j][i] = madd(Av, Bv[j], Cv[j][i]);
573                    }
574                }
575            }
576        }
577        for (int64_t j = 0; j < RN; ++j)
578            for (int64_t i = 0; i < RM; ++i)
579                C[ldc * (jj + j) + (ii + i)] = hsum(Cv[j][i]);
580    }
581
582    template <int RM, int RN, int BM>
583    NOINLINE void gemm(int64_t m, int64_t n, int64_t BN) {
584        GGML_ASSERT(m % (RM * BM) == 0);
585        const int64_t ytiles = m / (RM * BM);
586        const int64_t xtiles = (n + RN -1) / RN;
587        const int64_t jj_RN = (xtiles - (xtiles * RN - n));
588
589        // "round" bloc_size to "nearest" BN
590        const int64_t NB_BN = xtiles < BN ? 1 : (xtiles + BN / 2) / BN;
591        const int64_t SIZE_BN = xtiles % NB_BN == 0 ? xtiles / NB_BN : xtiles / NB_BN + 1;
592        const int64_t jj_BN = (NB_BN - (NB_BN * SIZE_BN - xtiles));
593        const int64_t nb_job = ytiles * NB_BN;
594
595        if (params->ith == 0) {
596            GGML_ASSERT( jj_BN * SIZE_BN + (NB_BN - jj_BN) * (SIZE_BN - 1) == xtiles);
597            // Every thread starts at ith, so the first unprocessed chunk is nth.  This save a bit of coordination right at the start.
598            ggml_threadpool_chunk_set(params->threadpool, params->nth);
599        }
600
601        ggml_barrier(params->threadpool);
602
603        int64_t job = params->ith;
604        while (job < nb_job) {
605            const int64_t ii = (job % ytiles) * RM * BM;
606            const int64_t jb =  job / ytiles;
607            const int64_t jr0 = BLOC_POS(jb  , jj_BN, SIZE_BN);
608            const int64_t jrN = BLOC_POS(jb+1, jj_BN, SIZE_BN);
609
610            const int64_t jj0 = BLOC_POS(jr0, jj_RN, RN);
611            const int64_t jj2 = BLOC_POS(jrN, jj_RN, RN);
612            const int64_t jj1 = jj2 < jj_RN * RN ? jj2 : jj_RN * RN;
613
614            for (int64_t bi = 0; bi < BM * RM; bi += RM) {
615                int64_t jj = jj0;
616                for (; jj < jj1; jj += RN) {
617                    gemm_bloc<RM, RN>(ii + bi, jj);
618                }
619                if constexpr (RN > 1) {
620                    for (; jj < jj2; jj += RN - 1) {
621                        gemm_bloc<RM, RN-1>(ii + bi, jj);
622                    }
623                }
624                GGML_ASSERT(jj == jj2);
625            }
626
627            job = ggml_threadpool_chunk_add(params->threadpool, 1);
628        }
629
630        ggml_barrier(params->threadpool);
631        return;
632    }
633
634    const ggml_compute_params * params;
635    const TA *const A;
636    const TB *const B;
637    TC *const C;
638    const int64_t k;
639    const int64_t lda;
640    const int64_t ldb;
641    const int64_t ldc;
642};
643
644#if defined(__riscv_v_intrinsic)
645template <typename D, typename V, typename TA, typename TB, typename TC>
646class tinyBLAS_RVV {
647  public:
648    tinyBLAS_RVV(const ggml_compute_params * params, int64_t k,
649             const TA *A, int64_t lda,
650             const TB *B, int64_t ldb,
651             TC *C, int64_t ldc)
652        : params(params), A(A), B(B), C(C), k(k), lda(lda), ldb(ldb), ldc(ldc) {
653    }
654
655    bool matmul(int64_t m, int64_t n) {
656        if (k % vlmax<V>() != 0) {
657            return false;
658        }
659
660#if LMUL == 1
661        if (m % 16 == 0 && (m/16 >= params->nth)) {
662            const int64_t SIZE_N = BLOCK_SIZE<6>(n);
663            mnpack<4, 6, 4>(m, n, SIZE_N, 12);
664            return true;
665        }
666        if (m % 8 == 0 ) {
667            const int64_t SIZE_N = BLOCK_SIZE<6>(n);
668            mnpack<4, 6, 2>(m, n, SIZE_N, 12);
669            return true;
670        }
671        if (m % 4 == 0) {
672            const int64_t SIZE_N = BLOCK_SIZE<6>(n);
673            mnpack<4, 6, 1>(m, n, SIZE_N, 12);
674            return true;
675        }
676#elif LMUL == 2
677        if (m % 16 == 0 && (m/16 >= params->nth)) {
678            const int64_t SIZE_N = BLOCK_SIZE<3>(n);
679            mnpack<4, 3, 4>(m, n, SIZE_N, 24);
680            return true;
681        }
682        if (m % 8 == 0 ) {
683            const int64_t SIZE_N = BLOCK_SIZE<3>(n);
684            mnpack<4, 3, 2>(m, n, SIZE_N, 24);
685            return true;
686        }
687        if (m % 4 == 0) {
688            const int64_t SIZE_N = BLOCK_SIZE<3>(n);
689            mnpack<4, 3, 1>(m, n, SIZE_N, 24);
690            return true;
691        }
692#else // LMUL = 4
693        if (m % 16 == 0 && (m/16 >= params->nth)) {
694            const int64_t SIZE_N = BLOCK_SIZE<2>(n);
695            mnpack<2, 2, 8>(m, n, SIZE_N, 36);
696            return true;
697        }
698        if (m % 8 == 0 ) {
699            const int64_t SIZE_N = BLOCK_SIZE<2>(n);
700            mnpack<2, 2, 4>(m, n, SIZE_N, 36);
701            return true;
702        }
703        if (m % 4 == 0) {
704            const int64_t SIZE_N = BLOCK_SIZE<2>(n);
705            mnpack<2, 2, 2>(m, n, SIZE_N, 36);
706            return true;
707        }
708#endif
709        return false;
710    }
711
712  private:
713    template<int RM, int RN, int BM>
714    inline void mnpack(int64_t m, int64_t n, int64_t SIZE_N, int64_t BN) {
715        if (SIZE_N == RN) {
716            return gemm<RM, RN, BM>(m, n, BN);
717        }
718        if constexpr (RN > 1) {
719            return mnpack<RM, RN-1, BM>(m, n, SIZE_N, BN);
720        } else {
721            GGML_LOG_ERROR("mnpack<%d, %d> block size not supported\n", RM, (int)SIZE_N);
722            GGML_ASSERT(false); // we have miss something.
723        }
724    }
725
726    inline void gemm_bloc_4x6(int64_t ii, int64_t jj) {
727        size_t vl = vlmax<V>();
728        D Cv00 = set_zero<D>();
729        D Cv01 = set_zero<D>();
730        D Cv02 = set_zero<D>();
731        D Cv03 = set_zero<D>();
732        D Cv10 = set_zero<D>();
733        D Cv11 = set_zero<D>();
734        D Cv12 = set_zero<D>();
735        D Cv13 = set_zero<D>();
736        D Cv20 = set_zero<D>();
737        D Cv21 = set_zero<D>();
738        D Cv22 = set_zero<D>();
739        D Cv23 = set_zero<D>();
740        D Cv30 = set_zero<D>();
741        D Cv31 = set_zero<D>();
742        D Cv32 = set_zero<D>();
743        D Cv33 = set_zero<D>();
744        D Cv40 = set_zero<D>();
745        D Cv41 = set_zero<D>();
746        D Cv42 = set_zero<D>();
747        D Cv43 = set_zero<D>();
748        D Cv50 = set_zero<D>();
749        D Cv51 = set_zero<D>();
750        D Cv52 = set_zero<D>();
751        D Cv53 = set_zero<D>();
752
753        for (int64_t l = 0; l < k; l += vl) {
754            V Bv0 = load<V>(B + ldb * (jj + 0) + l);
755            V Bv1 = load<V>(B + ldb * (jj + 1) + l);
756            V Bv2 = load<V>(B + ldb * (jj + 2) + l);
757            V Bv3 = load<V>(B + ldb * (jj + 3) + l);
758            V Bv4 = load<V>(B + ldb * (jj + 4) + l);
759            V Bv5 = load<V>(B + ldb * (jj + 5) + l);
760
761            V Av0 = load<V>(A + lda * (ii + 0) + l);
762            Cv00 = madd(Av0, Bv0, Cv00);
763            Cv10 = madd(Av0, Bv1, Cv10);
764            Cv20 = madd(Av0, Bv2, Cv20);
765            Cv30 = madd(Av0, Bv3, Cv30);
766            Cv40 = madd(Av0, Bv4, Cv40);
767            Cv50 = madd(Av0, Bv5, Cv50);
768
769            V Av1 = load<V>(A + lda * (ii + 1) + l);
770            Cv01 = madd(Av1, Bv0, Cv01);
771            Cv11 = madd(Av1, Bv1, Cv11);
772            Cv21 = madd(Av1, Bv2, Cv21);
773            Cv31 = madd(Av1, Bv3, Cv31);
774            Cv41 = madd(Av1, Bv4, Cv41);
775            Cv51 = madd(Av1, Bv5, Cv51);
776
777            V Av2 = load<V>(A + lda * (ii + 2) + l);
778            Cv02 = madd(Av2, Bv0, Cv02);
779            Cv12 = madd(Av2, Bv1, Cv12);
780            Cv22 = madd(Av2, Bv2, Cv22);
781            Cv32 = madd(Av2, Bv3, Cv32);
782            Cv42 = madd(Av2, Bv4, Cv42);
783            Cv52 = madd(Av2, Bv5, Cv52);
784
785            V Av3 = load<V>(A + lda * (ii + 3) + l);
786            Cv03 = madd(Av3, Bv0, Cv03);
787            Cv13 = madd(Av3, Bv1, Cv13);
788            Cv23 = madd(Av3, Bv2, Cv23);
789            Cv33 = madd(Av3, Bv3, Cv33);
790            Cv43 = madd(Av3, Bv4, Cv43);
791            Cv53 = madd(Av3, Bv5, Cv53);
792        }
793
794        C[ldc * (jj + 0) + (ii + 0)] = hsum(Cv00);
795        C[ldc * (jj + 0) + (ii + 1)] = hsum(Cv01);
796        C[ldc * (jj + 0) + (ii + 2)] = hsum(Cv02);
797        C[ldc * (jj + 0) + (ii + 3)] = hsum(Cv03);
798        C[ldc * (jj + 1) + (ii + 0)] = hsum(Cv10);
799        C[ldc * (jj + 1) + (ii + 1)] = hsum(Cv11);
800        C[ldc * (jj + 1) + (ii + 2)] = hsum(Cv12);
801        C[ldc * (jj + 1) + (ii + 3)] = hsum(Cv13);
802        C[ldc * (jj + 2) + (ii + 0)] = hsum(Cv20);
803        C[ldc * (jj + 2) + (ii + 1)] = hsum(Cv21);
804        C[ldc * (jj + 2) + (ii + 2)] = hsum(Cv22);
805        C[ldc * (jj + 2) + (ii + 3)] = hsum(Cv23);
806        C[ldc * (jj + 3) + (ii + 0)] = hsum(Cv30);
807        C[ldc * (jj + 3) + (ii + 1)] = hsum(Cv31);
808        C[ldc * (jj + 3) + (ii + 2)] = hsum(Cv32);
809        C[ldc * (jj + 3) + (ii + 3)] = hsum(Cv33);
810        C[ldc * (jj + 4) + (ii + 0)] = hsum(Cv40);
811        C[ldc * (jj + 4) + (ii + 1)] = hsum(Cv41);
812        C[ldc * (jj + 4) + (ii + 2)] = hsum(Cv42);
813        C[ldc * (jj + 4) + (ii + 3)] = hsum(Cv43);
814        C[ldc * (jj + 5) + (ii + 0)] = hsum(Cv50);
815        C[ldc * (jj + 5) + (ii + 1)] = hsum(Cv51);
816        C[ldc * (jj + 5) + (ii + 2)] = hsum(Cv52);
817        C[ldc * (jj + 5) + (ii + 3)] = hsum(Cv53);
818    }
819
820    inline void gemm_bloc_4x5(int64_t ii, int64_t jj) {
821        size_t vl = vlmax<V>();
822        D Cv00 = set_zero<D>();
823        D Cv01 = set_zero<D>();
824        D Cv02 = set_zero<D>();
825        D Cv03 = set_zero<D>();
826        D Cv10 = set_zero<D>();
827        D Cv11 = set_zero<D>();
828        D Cv12 = set_zero<D>();
829        D Cv13 = set_zero<D>();
830        D Cv20 = set_zero<D>();
831        D Cv21 = set_zero<D>();
832        D Cv22 = set_zero<D>();
833        D Cv23 = set_zero<D>();
834        D Cv30 = set_zero<D>();
835        D Cv31 = set_zero<D>();
836        D Cv32 = set_zero<D>();
837        D Cv33 = set_zero<D>();
838        D Cv40 = set_zero<D>();
839        D Cv41 = set_zero<D>();
840        D Cv42 = set_zero<D>();
841        D Cv43 = set_zero<D>();
842
843        for (int64_t l = 0; l < k; l += vl) {
844            V Bv0 = load<V>(B + ldb * (jj + 0) + l);
845            V Bv1 = load<V>(B + ldb * (jj + 1) + l);
846            V Bv2 = load<V>(B + ldb * (jj + 2) + l);
847            V Bv3 = load<V>(B + ldb * (jj + 3) + l);
848            V Bv4 = load<V>(B + ldb * (jj + 4) + l);
849
850            V Av0 = load<V>(A + lda * (ii + 0) + l);
851            Cv00 = madd(Av0, Bv0, Cv00);
852            Cv10 = madd(Av0, Bv1, Cv10);
853            Cv20 = madd(Av0, Bv2, Cv20);
854            Cv30 = madd(Av0, Bv3, Cv30);
855            Cv40 = madd(Av0, Bv4, Cv40);
856
857            V Av1 = load<V>(A + lda * (ii + 1) + l);
858            Cv01 = madd(Av1, Bv0, Cv01);
859            Cv11 = madd(Av1, Bv1, Cv11);
860            Cv21 = madd(Av1, Bv2, Cv21);
861            Cv31 = madd(Av1, Bv3, Cv31);
862            Cv41 = madd(Av1, Bv4, Cv41);
863
864            V Av2 = load<V>(A + lda * (ii + 2) + l);
865            Cv02 = madd(Av2, Bv0, Cv02);
866            Cv12 = madd(Av2, Bv1, Cv12);
867            Cv22 = madd(Av2, Bv2, Cv22);
868            Cv32 = madd(Av2, Bv3, Cv32);
869            Cv42 = madd(Av2, Bv4, Cv42);
870
871            V Av3 = load<V>(A + lda * (ii + 3) + l);
872            Cv03 = madd(Av3, Bv0, Cv03);
873            Cv13 = madd(Av3, Bv1, Cv13);
874            Cv23 = madd(Av3, Bv2, Cv23);
875            Cv33 = madd(Av3, Bv3, Cv33);
876            Cv43 = madd(Av3, Bv4, Cv43);
877        }
878
879        C[ldc * (jj + 0) + (ii + 0)] = hsum(Cv00);
880        C[ldc * (jj + 0) + (ii + 1)] = hsum(Cv01);
881        C[ldc * (jj + 0) + (ii + 2)] = hsum(Cv02);
882        C[ldc * (jj + 0) + (ii + 3)] = hsum(Cv03);
883        C[ldc * (jj + 1) + (ii + 0)] = hsum(Cv10);
884        C[ldc * (jj + 1) + (ii + 1)] = hsum(Cv11);
885        C[ldc * (jj + 1) + (ii + 2)] = hsum(Cv12);
886        C[ldc * (jj + 1) + (ii + 3)] = hsum(Cv13);
887        C[ldc * (jj + 2) + (ii + 0)] = hsum(Cv20);
888        C[ldc * (jj + 2) + (ii + 1)] = hsum(Cv21);
889        C[ldc * (jj + 2) + (ii + 2)] = hsum(Cv22);
890        C[ldc * (jj + 2) + (ii + 3)] = hsum(Cv23);
891        C[ldc * (jj + 3) + (ii + 0)] = hsum(Cv30);
892        C[ldc * (jj + 3) + (ii + 1)] = hsum(Cv31);
893        C[ldc * (jj + 3) + (ii + 2)] = hsum(Cv32);
894        C[ldc * (jj + 3) + (ii + 3)] = hsum(Cv33);
895        C[ldc * (jj + 4) + (ii + 0)] = hsum(Cv40);
896        C[ldc * (jj + 4) + (ii + 1)] = hsum(Cv41);
897        C[ldc * (jj + 4) + (ii + 2)] = hsum(Cv42);
898        C[ldc * (jj + 4) + (ii + 3)] = hsum(Cv43);
899    }
900
901    inline void gemm_bloc_4x4(int64_t ii, int64_t jj) {
902        size_t vl = vlmax<V>();
903        D Cv00 = set_zero<D>();
904        D Cv01 = set_zero<D>();
905        D Cv02 = set_zero<D>();
906        D Cv03 = set_zero<D>();
907        D Cv10 = set_zero<D>();
908        D Cv11 = set_zero<D>();
909        D Cv12 = set_zero<D>();
910        D Cv13 = set_zero<D>();
911        D Cv20 = set_zero<D>();
912        D Cv21 = set_zero<D>();
913        D Cv22 = set_zero<D>();
914        D Cv23 = set_zero<D>();
915        D Cv30 = set_zero<D>();
916        D Cv31 = set_zero<D>();
917        D Cv32 = set_zero<D>();
918        D Cv33 = set_zero<D>();
919
920        for (int64_t l = 0; l < k; l += vl) {
921            V Av0 = load<V>(A + lda * (ii + 0) + l);
922            V Av1 = load<V>(A + lda * (ii + 1) + l);
923            V Av2 = load<V>(A + lda * (ii + 2) + l);
924            V Av3 = load<V>(A + lda * (ii + 3) + l);
925
926            V Bv0 = load<V>(B + ldb * (jj + 0) + l);
927            Cv00 = madd(Av0, Bv0, Cv00);
928            Cv01 = madd(Av1, Bv0, Cv01);
929            Cv02 = madd(Av2, Bv0, Cv02);
930            Cv03 = madd(Av3, Bv0, Cv03);
931
932            V Bv1 = load<V>(B + ldb * (jj + 1) + l);
933            Cv10 = madd(Av0, Bv1, Cv10);
934            Cv11 = madd(Av1, Bv1, Cv11);
935            Cv12 = madd(Av2, Bv1, Cv12);
936            Cv13 = madd(Av3, Bv1, Cv13);
937
938            V Bv2 = load<V>(B + ldb * (jj + 2) + l);
939            Cv20 = madd(Av0, Bv2, Cv20);
940            Cv21 = madd(Av1, Bv2, Cv21);
941            Cv22 = madd(Av2, Bv2, Cv22);
942            Cv23 = madd(Av3, Bv2, Cv23);
943
944            V Bv3 = load<V>(B + ldb * (jj + 3) + l);
945            Cv30 = madd(Av0, Bv3, Cv30);
946            Cv31 = madd(Av1, Bv3, Cv31);
947            Cv32 = madd(Av2, Bv3, Cv32);
948            Cv33 = madd(Av3, Bv3, Cv33);
949        }
950
951        C[ldc * (jj + 0) + (ii + 0)] = hsum(Cv00);
952        C[ldc * (jj + 0) + (ii + 1)] = hsum(Cv01);
953        C[ldc * (jj + 0) + (ii + 2)] = hsum(Cv02);
954        C[ldc * (jj + 0) + (ii + 3)] = hsum(Cv03);
955        C[ldc * (jj + 1) + (ii + 0)] = hsum(Cv10);
956        C[ldc * (jj + 1) + (ii + 1)] = hsum(Cv11);
957        C[ldc * (jj + 1) + (ii + 2)] = hsum(Cv12);
958        C[ldc * (jj + 1) + (ii + 3)] = hsum(Cv13);
959        C[ldc * (jj + 2) + (ii + 0)] = hsum(Cv20);
960        C[ldc * (jj + 2) + (ii + 1)] = hsum(Cv21);
961        C[ldc * (jj + 2) + (ii + 2)] = hsum(Cv22);
962        C[ldc * (jj + 2) + (ii + 3)] = hsum(Cv23);
963        C[ldc * (jj + 3) + (ii + 0)] = hsum(Cv30);
964        C[ldc * (jj + 3) + (ii + 1)] = hsum(Cv31);
965        C[ldc * (jj + 3) + (ii + 2)] = hsum(Cv32);
966        C[ldc * (jj + 3) + (ii + 3)] = hsum(Cv33);
967    }
968
969    inline void gemm_bloc_4x3(int64_t ii, int64_t jj) {
970        size_t vl = vlmax<V>();
971        D Cv00 = set_zero<D>();
972        D Cv01 = set_zero<D>();
973        D Cv02 = set_zero<D>();
974        D Cv03 = set_zero<D>();
975        D Cv10 = set_zero<D>();
976        D Cv11 = set_zero<D>();
977        D Cv12 = set_zero<D>();
978        D Cv13 = set_zero<D>();
979        D Cv20 = set_zero<D>();
980        D Cv21 = set_zero<D>();
981        D Cv22 = set_zero<D>();
982        D Cv23 = set_zero<D>();
983
984        for (int64_t l = 0; l < k; l += vl) {
985            V Av0 = load<V>(A + lda * (ii + 0) + l);
986            V Av1 = load<V>(A + lda * (ii + 1) + l);
987            V Av2 = load<V>(A + lda * (ii + 2) + l);
988            V Av3 = load<V>(A + lda * (ii + 3) + l);
989
990            V Bv0 = load<V>(B + ldb * (jj + 0) + l);
991            Cv00 = madd(Av0, Bv0, Cv00);
992            Cv01 = madd(Av1, Bv0, Cv01);
993            Cv02 = madd(Av2, Bv0, Cv02);
994            Cv03 = madd(Av3, Bv0, Cv03);
995
996            V Bv1 = load<V>(B + ldb * (jj + 1) + l);
997            Cv10 = madd(Av0, Bv1, Cv10);
998            Cv11 = madd(Av1, Bv1, Cv11);
999            Cv12 = madd(Av2, Bv1, Cv12);
1000            Cv13 = madd(Av3, Bv1, Cv13);
1001
1002            V Bv2 = load<V>(B + ldb * (jj + 2) + l);
1003            Cv20 = madd(Av0, Bv2, Cv20);
1004            Cv21 = madd(Av1, Bv2, Cv21);
1005            Cv22 = madd(Av2, Bv2, Cv22);
1006            Cv23 = madd(Av3, Bv2, Cv23);
1007        }
1008
1009        C[ldc * (jj + 0) + (ii + 0)] = hsum(Cv00);
1010        C[ldc * (jj + 0) + (ii + 1)] = hsum(Cv01);
1011        C[ldc * (jj + 0) + (ii + 2)] = hsum(Cv02);
1012        C[ldc * (jj + 0) + (ii + 3)] = hsum(Cv03);
1013        C[ldc * (jj + 1) + (ii + 0)] = hsum(Cv10);
1014        C[ldc * (jj + 1) + (ii + 1)] = hsum(Cv11);
1015        C[ldc * (jj + 1) + (ii + 2)] = hsum(Cv12);
1016        C[ldc * (jj + 1) + (ii + 3)] = hsum(Cv13);
1017        C[ldc * (jj + 2) + (ii + 0)] = hsum(Cv20);
1018        C[ldc * (jj + 2) + (ii + 1)] = hsum(Cv21);
1019        C[ldc * (jj + 2) + (ii + 2)] = hsum(Cv22);
1020        C[ldc * (jj + 2) + (ii + 3)] = hsum(Cv23);
1021    }
1022
1023    inline void gemm_bloc_4x2(int64_t ii, int64_t jj) {
1024        size_t vl = vlmax<V>();
1025        D Cv00 = set_zero<D>();
1026        D Cv01 = set_zero<D>();
1027        D Cv02 = set_zero<D>();
1028        D Cv03 = set_zero<D>();
1029        D Cv10 = set_zero<D>();
1030        D Cv11 = set_zero<D>();
1031        D Cv12 = set_zero<D>();
1032        D Cv13 = set_zero<D>();
1033
1034        for (int64_t l = 0; l < k; l += vl) {
1035            V Av0 = load<V>(A + lda * (ii + 0) + l);
1036            V Av1 = load<V>(A + lda * (ii + 1) + l);
1037            V Av2 = load<V>(A + lda * (ii + 2) + l);
1038            V Av3 = load<V>(A + lda * (ii + 3) + l);
1039
1040            V Bv0 = load<V>(B + ldb * (jj + 0) + l);
1041            Cv00 = madd(Av0, Bv0, Cv00);
1042            Cv01 = madd(Av1, Bv0, Cv01);
1043            Cv02 = madd(Av2, Bv0, Cv02);
1044            Cv03 = madd(Av3, Bv0, Cv03);
1045
1046            V Bv1 = load<V>(B + ldb * (jj + 1) + l);
1047            Cv10 = madd(Av0, Bv1, Cv10);
1048            Cv11 = madd(Av1, Bv1, Cv11);
1049            Cv12 = madd(Av2, Bv1, Cv12);
1050            Cv13 = madd(Av3, Bv1, Cv13);
1051        }
1052
1053        C[ldc * (jj + 0) + (ii + 0)] = hsum(Cv00);
1054        C[ldc * (jj + 0) + (ii + 1)] = hsum(Cv01);
1055        C[ldc * (jj + 0) + (ii + 2)] = hsum(Cv02);
1056        C[ldc * (jj + 0) + (ii + 3)] = hsum(Cv03);
1057        C[ldc * (jj + 1) + (ii + 0)] = hsum(Cv10);
1058        C[ldc * (jj + 1) + (ii + 1)] = hsum(Cv11);
1059        C[ldc * (jj + 1) + (ii + 2)] = hsum(Cv12);
1060        C[ldc * (jj + 1) + (ii + 3)] = hsum(Cv13);
1061    }
1062
1063    inline void gemm_bloc_4x1(int64_t ii, int64_t jj) {
1064        size_t vl = vlmax<V>();
1065        D Cv00 = set_zero<D>();
1066        D Cv01 = set_zero<D>();
1067        D Cv02 = set_zero<D>();
1068        D Cv03 = set_zero<D>();
1069
1070        for (int64_t l = 0; l < k; l += vl) {
1071            V Av0 = load<V>(A + lda * (ii + 0) + l);
1072            V Av1 = load<V>(A + lda * (ii + 1) + l);
1073            V Av2 = load<V>(A + lda * (ii + 2) + l);
1074            V Av3 = load<V>(A + lda * (ii + 3) + l);
1075
1076            V Bv0 = load<V>(B + ldb * (jj + 0) + l);
1077            Cv00 = madd(Av0, Bv0, Cv00);
1078            Cv01 = madd(Av1, Bv0, Cv01);
1079            Cv02 = madd(Av2, Bv0, Cv02);
1080            Cv03 = madd(Av3, Bv0, Cv03);
1081        }
1082
1083        C[ldc * (jj + 0) + (ii + 0)] = hsum(Cv00);
1084        C[ldc * (jj + 0) + (ii + 1)] = hsum(Cv01);
1085        C[ldc * (jj + 0) + (ii + 2)] = hsum(Cv02);
1086        C[ldc * (jj + 0) + (ii + 3)] = hsum(Cv03);
1087    }
1088
1089    inline void gemm_bloc_2x2(int64_t ii, int64_t jj) {
1090        size_t vl = vlmax<V>();
1091        D Cv00 = set_zero<D>();
1092        D Cv01 = set_zero<D>();
1093        D Cv10 = set_zero<D>();
1094        D Cv11 = set_zero<D>();
1095
1096        for (int64_t l = 0; l < k; l += vl) {
1097            V Av0 = load<V>(A + lda * (ii + 0) + l);
1098            V Av1 = load<V>(A + lda * (ii + 1) + l);
1099
1100            V Bv0 = load<V>(B + ldb * (jj + 0) + l);
1101            Cv00 = madd(Av0, Bv0, Cv00);
1102            Cv01 = madd(Av1, Bv0, Cv01);
1103
1104            V Bv1 = load<V>(B + ldb * (jj + 1) + l);
1105            Cv10 = madd(Av0, Bv1, Cv10);
1106            Cv11 = madd(Av1, Bv1, Cv11);
1107        }
1108
1109        C[ldc * (jj + 0) + (ii + 0)] = hsum(Cv00);
1110        C[ldc * (jj + 0) + (ii + 1)] = hsum(Cv01);
1111        C[ldc * (jj + 1) + (ii + 0)] = hsum(Cv10);
1112        C[ldc * (jj + 1) + (ii + 1)] = hsum(Cv11);
1113    }
1114
1115    inline void gemm_bloc_2x1(int64_t ii, int64_t jj) {
1116        size_t vl = vlmax<V>();
1117        D Cv00 = set_zero<D>();
1118        D Cv01 = set_zero<D>();
1119
1120        for (int64_t l = 0; l < k; l += vl) {
1121            V Av0 = load<V>(A + lda * (ii + 0) + l);
1122            V Av1 = load<V>(A + lda * (ii + 1) + l);
1123
1124            V Bv0 = load<V>(B + ldb * (jj + 0) + l);
1125            Cv00 = madd(Av0, Bv0, Cv00);
1126            Cv01 = madd(Av1, Bv0, Cv01);
1127        }
1128
1129        C[ldc * (jj + 0) + (ii + 0)] = hsum(Cv00);
1130        C[ldc * (jj + 0) + (ii + 1)] = hsum(Cv01);
1131    }
1132
1133    template <int RM, int RN>
1134    inline void gemm_bloc(int64_t ii, int64_t jj) {
1135        if constexpr (RM == 4) {
1136            if constexpr (RN == 6) { return gemm_bloc_4x6(ii, jj); }
1137            if constexpr (RN == 5) { return gemm_bloc_4x5(ii, jj); }
1138            if constexpr (RN == 4) { return gemm_bloc_4x4(ii, jj); }
1139            if constexpr (RN == 3) { return gemm_bloc_4x3(ii, jj); }
1140            if constexpr (RN == 2) { return gemm_bloc_4x2(ii, jj); }
1141            if constexpr (RN == 1) { return gemm_bloc_4x1(ii, jj); }
1142        } else if constexpr (RM == 2) {
1143            if constexpr (RN == 2) { return gemm_bloc_2x2(ii, jj); }
1144            if constexpr (RN == 1) { return gemm_bloc_2x1(ii, jj); }
1145        }
1146    }
1147
1148    template <int RM, int RN, int BM>
1149    NOINLINE void gemm(int64_t m, int64_t n, int64_t BN) {
1150        GGML_ASSERT(m % (RM * BM) == 0);
1151        const int64_t ytiles = m / (RM * BM);
1152        const int64_t xtiles = (n + RN -1) / RN;
1153        const int64_t jj_RN = (xtiles - (xtiles * RN - n));
1154
1155        // "round" bloc_size to "nearest" BN
1156        const int64_t NB_BN = xtiles < BN ? 1 : (xtiles + BN / 2) / BN;
1157        const int64_t SIZE_BN = xtiles % NB_BN == 0 ? xtiles / NB_BN : xtiles / NB_BN + 1;
1158        const int64_t jj_BN = (NB_BN - (NB_BN * SIZE_BN - xtiles));
1159        const int64_t nb_job = ytiles * NB_BN;
1160
1161        if (params->ith == 0) {
1162            GGML_ASSERT( jj_BN * SIZE_BN + (NB_BN - jj_BN) * (SIZE_BN - 1) == xtiles);
1163            // Every thread starts at ith, so the first unprocessed chunk is nth.  This save a bit of coordination right at the start.
1164            ggml_threadpool_chunk_set(params->threadpool, params->nth);
1165        }
1166
1167        ggml_barrier(params->threadpool);
1168
1169        int64_t job = params->ith;
1170        while (job < nb_job) {
1171            const int64_t ii = (job % ytiles) * RM * BM;
1172            const int64_t jb =  job / ytiles;
1173            const int64_t jr0 = BLOC_POS(jb  , jj_BN, SIZE_BN);
1174            const int64_t jrN = BLOC_POS(jb+1, jj_BN, SIZE_BN);
1175
1176            const int64_t jj0 = BLOC_POS(jr0, jj_RN, RN);
1177            const int64_t jj2 = BLOC_POS(jrN, jj_RN, RN);
1178            const int64_t jj1 = jj2 < jj_RN * RN ? jj2 : jj_RN * RN;
1179
1180            for (int64_t bi = 0; bi < BM * RM; bi += RM) {
1181                int64_t jj = jj0;
1182                for (; jj < jj1; jj += RN) {
1183                    gemm_bloc<RM, RN>(ii + bi, jj);
1184                }
1185                if constexpr (RN > 1) {
1186                    for (; jj < jj2; jj += RN - 1) {
1187                        gemm_bloc<RM, RN-1>(ii + bi, jj);
1188                    }
1189                }
1190                GGML_ASSERT(jj == jj2);
1191            }
1192
1193            job = ggml_threadpool_chunk_add(params->threadpool, 1);
1194        }
1195
1196        ggml_barrier(params->threadpool);
1197        return;
1198    }
1199
1200    const ggml_compute_params * params;
1201    const TA *const A;
1202    const TB *const B;
1203    TC *const C;
1204    const int64_t k;
1205    const int64_t lda;
1206    const int64_t ldb;
1207    const int64_t ldc;
1208};
1209#endif
1210
1211//////////////////////////////////////////////////////////////////////////////////////////
1212// QUANT ZERO MATRIX MULTIPLICATION
1213
1214#if defined(__ARM_FEATURE_DOTPROD)
1215template <typename TA>
1216class tinyBLAS_Q0_ARM {
1217  public:
1218    tinyBLAS_Q0_ARM(int64_t k,
1219                    const TA *A, int64_t lda,
1220                    const block_q8_0 *B, int64_t ldb,
1221                    float *C, int64_t ldc,
1222                    int ith, int nth)
1223        : A(A), B(B), C(C), k(k), lda(lda), ldb(ldb), ldc(ldc), ith(ith), nth(nth) {
1224    }
1225
1226    void matmul(int64_t m, int64_t n) {
1227        mnpack(0, m, 0, n);
1228    }
1229
1230  private:
1231    NOINLINE void mnpack(int64_t m0, int64_t m, int64_t n0, int64_t n) {
1232        int64_t mc, nc, mp, np;
1233        switch ((MIN(m - m0, 3) << 4) | MIN(n - n0, 3ll)) {
1234        case 0x33:
1235            mc = 3;
1236            nc = 3;
1237            gemm<3, 3>(m0, m, n0, n);
1238            break;
1239        case 0x32:
1240            mc = 3;
1241            nc = 2;
1242            gemm<3, 2>(m0, m, n0, n);
1243            break;
1244        case 0x23:
1245            mc = 2;
1246            nc = 3;
1247            gemm<2, 3>(m0, m, n0, n);
1248            break;
1249        case 0x22:
1250            mc = 2;
1251            nc = 2;
1252            gemm<2, 2>(m0, m, n0, n);
1253            break;
1254        case 0x31:
1255            mc = 3;
1256            nc = 1;
1257            gemm<3, 1>(m0, m, n0, n);
1258            break;
1259        case 0x13:
1260            mc = 1;
1261            nc = 3;
1262            gemm<1, 3>(m0, m, n0, n);
1263            break;
1264        case 0x21:
1265            mc = 2;
1266            nc = 1;
1267            gemm<2, 1>(m0, m, n0, n);
1268            break;
1269        case 0x12:
1270            mc = 1;
1271            nc = 2;
1272            gemm<1, 2>(m0, m, n0, n);
1273            break;
1274        case 0x11:
1275            mc = 1;
1276            nc = 1;
1277            gemm<1, 1>(m0, m, n0, n);
1278            break;
1279        default:
1280            return;
1281        }
1282        mp = m0 + (m - m0) / mc * mc;
1283        np = n0 + (n - n0) / nc * nc;
1284        mnpack(mp, m, n0, np);
1285        mnpack(m0, m, np, n);
1286    }
1287
1288    template <int RM, int RN>
1289    NOINLINE void gemm(int64_t m0, int64_t m, int64_t n0, int64_t n) {
1290        int64_t ytiles = (m - m0) / RM;
1291        int64_t xtiles = (n - n0) / RN;
1292        int64_t tiles = xtiles * ytiles;
1293        int64_t duty = (tiles + nth - 1) / nth;
1294        int64_t start = duty * ith;
1295        int64_t end = start + duty;
1296        if (end > tiles)
1297            end = tiles;
1298        for (int64_t job = start; job < end; ++job) {
1299            int64_t ii = m0 + job / xtiles * RM;
1300            int64_t jj = n0 + job % xtiles * RN;
1301            float32x4_t Cv[RN][RM] = {};
1302            for (int64_t l = 0; l < k; ++l)
1303                for (int64_t j = 0; j < RN; ++j)
1304                    for (int64_t i = 0; i < RM; ++i)
1305                        Cv[j][i] = vmlaq_n_f32(Cv[j][i],
1306                                               vcvtq_f32_s32(vdotq_s32(
1307                                                   vdotq_s32(vdupq_n_s32(0),
1308                                                             load_lo(A + lda * (ii + i) + l),
1309                                                             load_lo(B + ldb * (jj + j) + l)),
1310                                                   load_hi(A + lda * (ii + i) + l),
1311                                                   load_hi(B + ldb * (jj + j) + l))),
1312                                               unhalf(A[lda * (ii + i) + l].d) *
1313                                               unhalf(B[ldb * (jj + j) + l].d));
1314            for (int64_t j = 0; j < RN; ++j)
1315                for (int64_t i = 0; i < RM; ++i)
1316                    C[ldc * (jj + j) + (ii + i)] = hsum(Cv[j][i]);
1317        }
1318    }
1319
1320    inline int8x16_t load_lo(const block_q8_0 *b) {
1321        return vld1q_s8(b->qs);
1322    }
1323
1324    inline int8x16_t load_hi(const block_q8_0 *b) {
1325        return vld1q_s8(b->qs + 16);
1326    }
1327
1328    inline int8x16_t load_lo(const block_q4_0 *b) {
1329        return vsubq_s8(vreinterpretq_s8_u8(vandq_u8(vld1q_u8(b->qs),
1330                                                     vdupq_n_u8(0x0f))),
1331                        vdupq_n_s8(0x8));
1332    }
1333
1334    inline int8x16_t load_hi(const block_q4_0 *b) {
1335        return vsubq_s8(vreinterpretq_s8_u8(vshrq_n_u8(vld1q_u8(b->qs), 4)),
1336                        vdupq_n_s8(0x8));
1337    }
1338
1339    const TA *const A;
1340    const block_q8_0 *const B;
1341    float *const C;
1342    const int64_t k;
1343    const int64_t lda;
1344    const int64_t ldb;
1345    const int64_t ldc;
1346    const int ith;
1347    const int nth;
1348};
1349#endif // __ARM_FEATURE_DOTPROD
1350
1351#if defined(__AVX2__) || defined(__AVX512F__) || defined(__AVX__)
1352template <typename TA, typename TB, typename TC>
1353class tinyBLAS_Q0_AVX {
1354  public:
1355    tinyBLAS_Q0_AVX(int64_t k,
1356                    const TA *A, int64_t lda,
1357                    const TB *B, int64_t ldb,
1358                    TC *C, int64_t ldc,
1359                    int ith, int nth)
1360        : A(A), B(B), C(C), k(k), lda(lda), ldb(ldb), ldc(ldc), ith(ith), nth(nth) {
1361        const int8_t kvalues_iq4nl[16] = {
1362            -127, -104, -83, -65,
1363            -49,  -35,  -22, -10,
1364              1,   13,   25,  38,
1365             53,   69,   89, 113
1366        };
1367
1368        iq4nlt = _mm_loadu_si128((const __m128i *)kvalues_iq4nl);
1369    }
1370
1371    void matmul(int64_t m, int64_t n) {
1372        mnpack(0, m, 0, n);
1373    }
1374
1375  private:
1376    void mnpack(int64_t m0, int64_t m, int64_t n0, int64_t n) {
1377        int64_t mc, nc, mp, np;
1378        switch ((MIN(m - m0, 4) << 4) | MIN(n - n0, 4)) {
1379#if VECTOR_REGISTERS == 32
1380        case 0x44:
1381            mc = 4;
1382            nc = 4;
1383#if defined(__AVX2__) && defined(__F16C__)
1384            gemm4xN<4>(m0, m, n0, n);
1385#else
1386            gemm<4, 4>(m0, m, n0, n);
1387#endif
1388            break;
1389        case 0x43:
1390            mc = 4;
1391            nc = 3;
1392#if defined(__AVX2__) && defined(__F16C__)
1393            gemm4xN<3>(m0, m, n0, n);
1394#else
1395            gemm<4, 3>(m0, m, n0, n);
1396#endif
1397            break;
1398        case 0x34:
1399            mc = 3;
1400            nc = 4;
1401#if defined(__AVX2__) && defined(__F16C__)
1402            gemmMx4<3>(m0, m, n0, n);
1403#else
1404            gemm<3, 4>(m0, m, n0, n);
1405#endif
1406            break;
1407        case 0x33:
1408            mc = 3;
1409            nc = 3;
1410            gemm<3, 3>(m0, m, n0, n);
1411            break;
1412        case 0x42:
1413            mc = 4;
1414            nc = 2;
1415#if defined(__AVX2__) && defined(__F16C__)
1416            gemm4xN<2>(m0, m, n0, n);
1417#else
1418            gemm<4, 2>(m0, m, n0, n);
1419#endif
1420            break;
1421        case 0x24:
1422            mc = 2;
1423            nc = 4;
1424#if defined(__AVX2__) && defined(__F16C__)
1425            gemmMx4<2>(m0, m, n0, n);
1426#else
1427            gemm<2, 4>(m0, m, n0, n);
1428#endif
1429            break;
1430#else
1431        case 0x44:
1432        case 0x43:
1433        case 0x42:
1434            mc = 4;
1435            nc = 2;
1436#if defined(__AVX2__) && defined(__F16C__)
1437            gemm4xN<2>(m0, m, n0, n);
1438#else
1439            gemm<4, 2>(m0, m, n0, n);
1440#endif
1441            break;
1442        case 0x34:
1443        case 0x24:
1444            mc = 2;
1445            nc = 4;
1446#if defined(__AVX2__) && defined(__F16C__)
1447            gemmMx4<2>(m0, m, n0, n);
1448#else
1449            gemm<2, 4>(m0, m, n0, n);
1450#endif
1451            break;
1452        case 0x33:
1453#endif
1454        case 0x32:
1455            mc = 3;
1456            nc = 2;
1457            gemm<3, 2>(m0, m, n0, n);
1458            break;
1459        case 0x23:
1460            mc = 2;
1461            nc = 3;
1462            gemm<2, 3>(m0, m, n0, n);
1463            break;
1464        case 0x41:
1465            mc = 4;
1466            nc = 1;
1467#if defined(__AVX2__) && defined(__F16C__)
1468            gemm4xN<1>(m0, m, n0, n);
1469#else
1470            gemm<4, 1>(m0, m, n0, n);
1471#endif
1472            break;
1473        case 0x22:
1474            mc = 2;
1475            nc = 2;
1476            gemm<2, 2>(m0, m, n0, n);
1477            break;
1478        case 0x14:
1479            mc = 1;
1480            nc = 4;
1481#if defined(__AVX2__) && defined(__F16C__)
1482            gemmMx4<1>(m0, m, n0, n);
1483#else
1484            gemm<1, 4>(m0, m, n0, n);
1485#endif
1486            break;
1487        case 0x31:
1488            mc = 3;
1489            nc = 1;
1490            gemm<3, 1>(m0, m, n0, n);
1491            break;
1492        case 0x13:
1493            mc = 1;
1494            nc = 3;
1495            gemm<1, 3>(m0, m, n0, n);
1496            break;
1497        case 0x21:
1498            mc = 2;
1499            nc = 1;
1500            gemm<2, 1>(m0, m, n0, n);
1501            break;
1502        case 0x12:
1503            mc = 1;
1504            nc = 2;
1505            gemm<1, 2>(m0, m, n0, n);
1506            break;
1507        case 0x11:
1508            mc = 1;
1509            nc = 1;
1510            gemm<1, 1>(m0, m, n0, n);
1511            break;
1512        default:
1513            return;
1514        }
1515        mp = m0 + (m - m0) / mc * mc;
1516        np = n0 + (n - n0) / nc * nc;
1517        mnpack(mp, m, n0, np);
1518        mnpack(m0, m, np, n);
1519    }
1520
1521#if defined(__AVX2__) && defined(__F16C__)
1522// Templated functions for gemm of dimensions 4xN
1523    template <int RN>
1524    NOINLINE void gemm4xN(int64_t m0, int64_t m, int64_t n0, int64_t n) {
1525        int64_t ytiles = (m - m0) / 4;
1526        int64_t xtiles = (n - n0) / RN;
1527        int64_t tiles = xtiles * ytiles;
1528        int64_t duty = (tiles + nth - 1) / nth;
1529        int64_t start = duty * ith;
1530        int64_t end = start + duty;
1531        if (end > tiles)
1532            end = tiles;
1533        for (int64_t job = start; job < end; ++job) {
1534            int64_t ii = m0 + job / xtiles * 4;
1535            int64_t jj = n0 + job % xtiles * RN;
1536            __m256 Cv[RN][4] = {};
1537            for (int64_t l = 0; l < k; ++l) {
1538                uint64_t a_delta = ((uint64_t)A[lda * (ii + 3) + l].d << 48) | ((uint64_t)A[lda * (ii + 2) + l].d << 32) | ((uint64_t)A[lda * (ii + 1) + l].d << 16) | (A[lda * (ii + 0) + l].d);
1539                // Convert delta values for four blocks to float values
1540                __m128 da = _mm_cvtph_ps(_mm_set_epi64x(0, a_delta));
1541                __m256i avec0 = load(A + lda * (ii + 0) + l);
1542                __m256i avec1 = load(A + lda * (ii + 1) + l);
1543                __m256i avec2 = load(A + lda * (ii + 2) + l);
1544                __m256i avec3 = load(A + lda * (ii + 3) + l);
1545                for (int64_t j = 0; j < RN; ++j) {
1546                        __m128 db = _mm_set1_ps(unhalf(B[ldb * (jj + j) + l].d));
1547                        // Computation of product of delta values for four blocks and replicate it across 256 bit lane
1548                        __m256 dvec =  _mm256_castps128_ps256(_mm_mul_ps(da, db));
1549                        dvec = _mm256_permute2f128_ps(dvec ,dvec, 0);
1550                        // Computation of dot product and multiplication with appropriate delta value products
1551                        Cv[j][0] = madd(_mm256_shuffle_ps(dvec, dvec, 0),
1552                                    updot(_mm256_sign_epi8(avec0, avec0),
1553                                          _mm256_sign_epi8(load(B + ldb * (jj + j) + l), avec0)),
1554                                    Cv[j][0]);
1555                        Cv[j][1] = madd(_mm256_shuffle_ps(dvec, dvec, 85),
1556                                    updot(_mm256_sign_epi8(avec1, avec1),
1557                                            _mm256_sign_epi8(load(B + ldb * (jj + j) + l), avec1)),
1558                                    Cv[j][1]);
1559                        Cv[j][2] = madd(_mm256_shuffle_ps(dvec, dvec, 170),
1560                                    updot(_mm256_sign_epi8(avec2, avec2),
1561                                            _mm256_sign_epi8(load(B + ldb * (jj + j) + l), avec2)),
1562                                    Cv[j][2]);
1563                        Cv[j][3] = madd(_mm256_shuffle_ps(dvec, dvec, 255),
1564                                    updot(_mm256_sign_epi8(avec3, avec3),
1565                                            _mm256_sign_epi8(load(B + ldb * (jj + j) + l), avec3)),
1566                                    Cv[j][3]);
1567                }
1568            }
1569
1570            for (int64_t j = 0; j < RN; ++j)
1571                for (int64_t i = 0; i < 4; ++i)
1572                    C[ldc * (jj + j) + (ii + i)] = hsum(Cv[j][i]);
1573        }
1574    }
1575
1576    // Templated functions for gemm of dimensions Mx4
1577    template <int RM>
1578    NOINLINE void gemmMx4(int64_t m0, int64_t m, int64_t n0, int64_t n) {
1579        int64_t ytiles = (m - m0) / RM;
1580        int64_t xtiles = (n - n0) / 4;
1581        int64_t tiles = xtiles * ytiles;
1582        int64_t duty = (tiles + nth - 1) / nth;
1583        int64_t start = duty * ith;
1584        int64_t end = start + duty;
1585        if (end > tiles)
1586            end = tiles;
1587        for (int64_t job = start; job < end; ++job) {
1588            int64_t ii = m0 + job / xtiles * RM;
1589            int64_t jj = n0 + job % xtiles * 4;
1590            __m256 Cv[4][RM] = {};
1591            for (int64_t l = 0; l < k; ++l) {
1592                uint64_t b_delta = ((uint64_t)B[ldb * (jj + 3) + l].d << 48) | ((uint64_t)B[ldb * (jj + 2) + l].d << 32) | ((uint64_t)B[ldb * (jj + 1) + l].d << 16) | (B[ldb * (jj + 0) + l].d);
1593                // Convert delta values for four blocks to float values
1594                __m128 db = _mm_cvtph_ps(_mm_set_epi64x(0, b_delta));
1595                __m256i bvec0 = load(B + ldb * (jj + 0) + l);
1596                __m256i bvec1 = load(B + ldb * (jj + 1) + l);
1597                __m256i bvec2 = load(B + ldb * (jj + 2) + l);
1598                __m256i bvec3 = load(B + ldb * (jj + 3) + l);
1599                for (int64_t i = 0; i < RM; ++i) {
1600                    __m128 da = _mm_set1_ps(unhalf((A[lda * (ii + i) + l].d)));
1601                    // Computation of product of delta values for four blocks and replicate it across 256 bit lane
1602                    __m256 dvec =  _mm256_castps128_ps256(_mm_mul_ps(da, db));
1603                    dvec = _mm256_permute2f128_ps(dvec ,dvec, 0);
1604                    // Computation of dot product and multiplication with appropriate delta value products
1605                    Cv[0][i] = madd(_mm256_shuffle_ps(dvec, dvec, 0),
1606                                    updot(_mm256_sign_epi8(load(A + lda * (ii + i) + l),
1607                                                            load(A + lda * (ii + i) + l)),
1608                                            _mm256_sign_epi8(bvec0, load(A + lda * (ii + i) + l))),
1609                                    Cv[0][i]);
1610                    Cv[1][i] = madd(_mm256_shuffle_ps(dvec, dvec, 85),
1611                                    updot(_mm256_sign_epi8(load(A + lda * (ii + i) + l),
1612                                                            load(A + lda * (ii + i) + l)),
1613                                            _mm256_sign_epi8(bvec1, load(A + lda * (ii + i) + l))),
1614                                    Cv[1][i]);
1615                    Cv[2][i] = madd(_mm256_shuffle_ps(dvec, dvec, 170),
1616                                    updot(_mm256_sign_epi8(load(A + lda * (ii + i) + l),
1617                                                            load(A + lda * (ii + i) + l)),
1618                                            _mm256_sign_epi8(bvec2, load(A + lda * (ii + i) + l))),
1619                                    Cv[2][i]);
1620                    Cv[3][i] = madd(_mm256_shuffle_ps(dvec, dvec, 255),
1621                                    updot(_mm256_sign_epi8(load(A + lda * (ii + i) + l),
1622                                                            load(A + lda * (ii + i) + l)),
1623                                            _mm256_sign_epi8(bvec3, load(A + lda * (ii + i) + l))),
1624                                    Cv[3][i]);
1625                }
1626            }
1627            for (int64_t j = 0; j < 4; ++j)
1628                for (int64_t i = 0; i < RM; ++i)
1629                    C[ldc * (jj + j) + (ii + i)] = hsum(Cv[j][i]);
1630        }
1631    }
1632#endif
1633
1634    template <int RM, int RN>
1635    NOINLINE void gemm(int64_t m0, int64_t m, int64_t n0, int64_t n) {
1636        int64_t ytiles = (m - m0) / RM;
1637        int64_t xtiles = (n - n0) / RN;
1638        int64_t tiles = xtiles * ytiles;
1639        int64_t duty = (tiles + nth - 1) / nth;
1640        int64_t start = duty * ith;
1641        int64_t end = start + duty;
1642        if (end > tiles)
1643            end = tiles;
1644        for (int64_t job = start; job < end; ++job) {
1645            int64_t ii = m0 + job / xtiles * RM;
1646            int64_t jj = n0 + job % xtiles * RN;
1647            __m256 Cv[RN][RM] = {};
1648            for (int64_t l = 0; l < k; ++l)
1649                for (int64_t j = 0; j < RN; ++j)
1650                    for (int64_t i = 0; i < RM; ++i) {
1651#if defined(__AVX2__)
1652                        __m256 udTmp = updot(_mm256_sign_epi8(load(A + lda * (ii + i) + l),
1653                                                              load(A + lda * (ii + i) + l)),
1654                                             _mm256_sign_epi8(load(B + ldb * (jj + j) + l),
1655                                                              load(A + lda * (ii + i) + l)));
1656#else
1657                        __m128i ali0 = load0(A + lda * (ii + i) + l);
1658                        __m128i ali1 = load1(A + lda * (ii + i) + l);
1659                        __m128i blj0 = load0(B + ldb * (jj + j) + l);
1660                        __m128i blj1 = load1(B + ldb * (jj + j) + l);
1661
1662                        __m128i sepAA0 = _mm_sign_epi8(ali0, ali0);
1663                        __m128i sepAA1 = _mm_sign_epi8(ali1, ali1);
1664                        __m128i sepBA0 = _mm_sign_epi8(blj0, ali0);
1665                        __m128i sepBA1 = _mm_sign_epi8(blj1, ali1);
1666
1667                        // updot
1668                        const __m128i oneFill = _mm_set1_epi16(1);
1669                        __m128i mad0 = _mm_maddubs_epi16(sepAA0, sepBA0);
1670                        __m128i mad1 = _mm_maddubs_epi16(sepAA1, sepBA1);
1671                        __m256 udTmp = _mm256_cvtepi32_ps(MM256_SET_M128I(_mm_madd_epi16(oneFill, mad1), _mm_madd_epi16(oneFill, mad0)));
1672#endif
1673                        Cv[j][i] = madd(_mm256_set1_ps(unhalf(A[lda * (ii + i) + l].d) *
1674                                                       unhalf(B[ldb * (jj + j) + l].d)),
1675                                                       udTmp,
1676                                                       Cv[j][i]);
1677                    }
1678            for (int64_t j = 0; j < RN; ++j)
1679                for (int64_t i = 0; i < RM; ++i)
1680                    C[ldc * (jj + j) + (ii + i)] = hsum(Cv[j][i]);
1681        }
1682    }
1683
1684    inline __m256i load(const block_q8_0 *b) {
1685        return _mm256_loadu_si256((const __m256i *)b->qs);
1686    }
1687
1688    inline __m128i load0(const block_q8_0 *b) {
1689        return _mm_loadu_si128((const __m128i *)b->qs);
1690    }
1691
1692    inline __m128i load1(const block_q8_0 *b) {
1693        return _mm_loadu_si128(((const __m128i *)b->qs) + 1);
1694    }
1695
1696    inline __m256i load(const block_q4_0 *b) {
1697        return _mm256_sub_epi8(denibble(b->qs), _mm256_set1_epi8(8));
1698    }
1699
1700    inline __m128i load0(const block_q4_0 *b) {
1701        const __m128i x = _mm_loadu_si128((const __m128i *)(b->qs));
1702        return _mm_sub_epi8(_mm_and_si128(_mm_set1_epi8(15), x), _mm_set1_epi8(8));
1703    }
1704
1705    inline __m128i load1(const block_q4_0 *b) {
1706        const __m128i x = _mm_loadu_si128((const __m128i *)(b->qs));
1707        return _mm_sub_epi8(_mm_and_si128(_mm_set1_epi8(15), _mm_srli_epi16(x, 4)), _mm_set1_epi8(8));
1708    }
1709
1710    inline __m256i load(const block_q5_0 *b) {
1711        return _mm256_or_si256(denibble(b->qs), bittobyte(b->qh));
1712    }
1713
1714    inline __m128i load0(const block_q5_0* b) {
1715        const __m128i x = _mm_loadu_si128((const __m128i *)(b->qs));
1716        uint32_t x32;
1717        memcpy(&x32, b->qh, sizeof(uint32_t));
1718        __m128i qxl = _mm_and_si128(_mm_set1_epi8(15), x);
1719        __m128i bytesl = _mm_cmpeq_epi8(_mm_set1_epi64x(-1),
1720                                        _mm_or_si128(_mm_set1_epi64x(0x7fbfdfeff7fbfdfe),
1721                                                     _mm_shuffle_epi8(_mm_set1_epi32(x32),
1722                                                                      _mm_set_epi64x(0x0101010101010101, 0x0000000000000000))));
1723        bytesl = _mm_andnot_si128(bytesl, _mm_set1_epi8((char)0xF0));
1724        return _mm_or_si128(qxl, bytesl);
1725    }
1726
1727    inline __m128i load1(const block_q5_0* b) {
1728        const __m128i x = _mm_loadu_si128((const __m128i *)(b->qs));
1729        uint32_t x32;
1730        memcpy(&x32, b->qh, sizeof(uint32_t));
1731        __m128i qxh = _mm_and_si128(_mm_set1_epi8(15), _mm_srli_epi16(x, 4));
1732        __m128i bytesh = _mm_cmpeq_epi8(_mm_set1_epi64x(-1),
1733                                        _mm_or_si128(_mm_set1_epi64x(0x7fbfdfeff7fbfdfe),
1734                                                     _mm_shuffle_epi8(_mm_set1_epi32(x32),
1735                                                                      _mm_set_epi64x(0x0303030303030303, 0x0202020202020202))));
1736        bytesh = _mm_andnot_si128(bytesh, _mm_set1_epi8((char)0xF0));
1737        return _mm_or_si128(qxh, bytesh);
1738    }
1739
1740    inline __m256i load(const block_iq4_nl *b) {
1741        return MM256_SET_M128I(load1(b), load0(b));
1742    }
1743
1744    inline __m128i load0(const block_iq4_nl *b) {
1745        const __m128i x = _mm_loadu_si128((const __m128i *)(b->qs));
1746        return _mm_shuffle_epi8(iq4nlt, _mm_and_si128(_mm_set1_epi8(15), x));
1747    }
1748
1749    inline __m128i load1(const block_iq4_nl *b) {
1750        const __m128i x = _mm_loadu_si128((const __m128i *)(b->qs));
1751        return _mm_shuffle_epi8(iq4nlt, _mm_and_si128(_mm_set1_epi8(15), _mm_srli_epi16(x, 4)));
1752    }
1753
1754    inline __m256 updot(__m256i u, __m256i s) {
1755        __m256i res;
1756#if defined(__AVX512VNNI__) && defined(__AVX512VL__)
1757        res = _mm256_dpbusd_epi32(_mm256_setzero_si256(), u, s);
1758#elif defined(__AVXVNNI__)
1759        res = _mm256_dpbusd_avx_epi32(_mm256_setzero_si256(), u, s);
1760#else
1761        res = _mm256_madd_epi16(_mm256_set1_epi16(1), _mm256_maddubs_epi16(u, s));
1762#endif
1763        return _mm256_cvtepi32_ps(res);
1764    }
1765
1766    static inline __m256i denibble(const uint8_t *p) {
1767        __m128i x = _mm_loadu_si128((const __m128i *)p);
1768        return _mm256_and_si256(_mm256_set1_epi8(15),
1769                                _mm256_insertf128_si256(_mm256_castsi128_si256(x),
1770                                                        _mm_srli_epi16(x, 4), 1));
1771    }
1772
1773    static inline __m256i bittobyte(const uint8_t *p) {
1774        uint32_t x32;
1775        memcpy(&x32, p, sizeof(uint32_t));
1776        __m256i bytes = _mm256_cmpeq_epi8(_mm256_set1_epi64x(-1),
1777                                          _mm256_or_si256(_mm256_set1_epi64x(0x7fbfdfeff7fbfdfe),
1778                                                          _mm256_shuffle_epi8(_mm256_set1_epi32(x32),
1779                                                                              _mm256_set_epi64x(0x0303030303030303, 0x0202020202020202,
1780                                                                                                0x0101010101010101, 0x0000000000000000))));
1781        return _mm256_andnot_si256(bytes, _mm256_set1_epi8((char)0xF0));
1782    }
1783
1784    const TA *const A;
1785    const TB *const B;
1786    TC *const C;
1787    const int64_t k;
1788    const int64_t lda;
1789    const int64_t ldb;
1790    const int64_t ldc;
1791    const int ith;
1792    const int nth;
1793    __m128i iq4nlt;
1794};
1795#endif // __AVX__
1796
1797//PPC Implementation
1798#if defined(__MMA__)
1799
1800template<typename T>
1801struct mma_instr;
1802
1803template<>
1804struct mma_instr<ggml_bf16_t> {
1805    static inline void outer_product(acc_t *acc, vec_t a, vec_t b) {
1806        __builtin_mma_xvbf16ger2pp(acc, a, b);
1807    }
1808};
1809
1810template<>
1811struct mma_instr<ggml_fp16_t> {
1812    static inline void outer_product(acc_t *acc, vec_t a, vec_t b) {
1813        __builtin_mma_xvf16ger2pp(acc, a, b);
1814    }
1815};
1816
1817template <typename TA, typename TB, typename TC>
1818class tinyBLAS_HP16_PPC {
1819  public:
1820    tinyBLAS_HP16_PPC(int64_t k,
1821                const TA *A, int64_t lda,
1822                const TB *B, int64_t ldb,
1823                TC *C, int64_t ldc,
1824                int ith, int nth)
1825        : A(A), B(B), C(C), k(k), lda(lda), ldb(ldb), ldc(ldc), ith(ith), nth(nth) {
1826    }
1827
1828    void matmul(int64_t m, int64_t n) {
1829        int64_t mc = 256;
1830        int64_t nc = 256;
1831        int64_t kc = 256;
1832    #if defined(_AIX) || defined(__BIG_ENDIAN__)
1833        mc = 128;
1834        nc = 128;
1835        kc = 128;
1836    #endif
1837        if (k < kc) {
1838            kc = k;
1839        }
1840        bool can_use_tiled = (m % mc == 0) && (n % nc == 0) && (k % kc == 0);
1841        if (can_use_tiled) {
1842            matmul_tiled(m, n, mc, nc, kc);
1843        } else {
1844            mnpack(0, m, 0, n);
1845        }
1846    }
1847
1848  private:
1849    __attribute__((always_inline))
1850    inline void save_acc(acc_t * ACC, int64_t ii, int64_t jj) {
1851        vec_t vec_C[4];
1852        __builtin_mma_disassemble_acc(vec_C, ACC);
1853        for (int I = 0; I < 4; I++) {
1854            for (int J = 0; J < 4; J++) {
1855                *((float *)(C+ii+((jj+J)*ldc)+I)) = *((float *)&vec_C[I]+J);
1856            }
1857        }
1858    }
1859
1860    __attribute__((always_inline))
1861    inline void add_save_acc(acc_t * ACC, int64_t ii, int64_t jj) {
1862        vec_t vec_C[4];
1863        __builtin_mma_disassemble_acc(vec_C, ACC);
1864        for (int I = 0; I < 4; I++) {
1865            for (int J = 0; J < 4; J++) {
1866                float * c_ptr = (float *)(C+ii+((jj+J)*ldc)+I);
1867                *c_ptr += *((float *)&vec_C[I]+J);
1868            }
1869        }
1870    }
1871
1872    void vector_permute_store(vec_t *c, int numVec, unsigned char *vecOffset) {
1873        vec_t t[8], s[8];
1874        vec_t swiz1 = {0, 1, 2, 3, 16, 17, 18, 19, 4, 5, 6, 7, 20, 21, 22, 23};
1875        vec_t swiz2 = {8, 9, 10, 11, 24, 25, 26, 27, 12, 13, 14, 15, 28, 29, 30, 31};
1876        vec_t swiz3 = {0, 1, 2, 3, 4, 5, 6, 7, 16, 17, 18, 19, 20, 21, 22, 23};
1877        vec_t swiz4 = {8, 9, 10, 11, 12, 13, 14, 15, 24, 25, 26, 27, 28, 29, 30, 31};
1878
1879        if (numVec == 2) {
1880            t[0] = vec_perm(c[0], c[1], swiz1);
1881            t[1] = vec_perm(c[2], c[3], swiz1);
1882            s[0] = vec_perm(t[0], t[1], swiz3);
1883            s[1] = vec_perm(t[0], t[1], swiz4);
1884            vec_xst(s[0], 0, (vec_t*)vecOffset);
1885            vec_xst(s[1], 0, (vec_t*)(vecOffset + 16));
1886        } else if (numVec == 4) {
1887            t[0] = vec_perm(c[0], c[1], swiz1);
1888            t[1] = vec_perm(c[0], c[1], swiz2);
1889            t[2] = vec_perm(c[2], c[3], swiz1);
1890            t[3] = vec_perm(c[2], c[3], swiz2);
1891            s[0] = vec_perm(t[0], t[2], swiz3);
1892            s[1] = vec_perm(t[0], t[2], swiz4);
1893            s[2] = vec_perm(t[1], t[3], swiz3);
1894            s[3] = vec_perm(t[1], t[3], swiz4);
1895            for (int i = 0; i < 4; ++i)
1896                vec_xst(s[i], 0, (vec_t*)(vecOffset + i * 16));
1897        } else if (numVec == 8) {
1898            for (int i = 0; i < 4; i += 2) {
1899                t[i+0] = vec_perm(c[i+0], c[i+1], swiz1);
1900                t[i+1] = vec_perm(c[i+0], c[i+1], swiz2);
1901            }
1902            for (int i = 4; i < 8; i += 2) {
1903                t[i+0] = vec_perm(c[i+0], c[i+1], swiz1);
1904                t[i+1] = vec_perm(c[i+0], c[i+1], swiz2);
1905            }
1906            s[0] = vec_perm(t[0], t[2], swiz3);
1907            s[1] = vec_perm(t[0], t[2], swiz4);
1908            s[2] = vec_perm(t[1], t[3], swiz3);
1909            s[3] = vec_perm(t[1], t[3], swiz4);
1910            s[4] = vec_perm(t[4], t[6], swiz3);
1911            s[5] = vec_perm(t[4], t[6], swiz4);
1912            s[6] = vec_perm(t[5], t[7], swiz3);
1913            s[7] = vec_perm(t[5], t[7], swiz4);
1914            for (int i = 0; i < 8; ++i)
1915                vec_xst(s[i], 0, (vec_t*)(vecOffset + i * 16));
1916        }
1917    }
1918
1919    void packNormal(const TA* a, int64_t lda, int rows, int cols, unsigned char* vec) {
1920        int64_t i, j;
1921        TA *aoffset = NULL;
1922        unsigned char *vecOffset = NULL;
1923        TA * aoffsets[8];
1924        vector unsigned char c_arr[8];
1925        aoffset = const_cast<TA*>(a);
1926        vecOffset = vec;
1927        j = (rows >> 3);
1928        if (j > 0) {
1929            do {
1930                aoffsets[0] = aoffset;
1931                if (cols == 4) {
1932                    aoffsets[0] = aoffset;
1933                    for (int it = 1; it < 4; ++it)
1934                        aoffsets[it] = aoffsets[it-1] + lda;
1935                    aoffset += 4 * lda;
1936                    for (int i = 0; i < 4; ++i)
1937                        c_arr[i] = vec_xl(0, (vector unsigned char*)aoffsets[i]);
1938                    vector_permute_store(c_arr, 4, vecOffset);
1939                    for (int i = 0; i<4; i++)
1940                        aoffsets[i] = aoffsets[i]+lda;
1941                    vecOffset +=64;
1942                }
1943                i = (cols >> 3);
1944                if (i > 0) {
1945                    for (int it = 1; it < 8; ++it) {
1946                        aoffsets[it] = aoffsets[it-1] + lda;
1947                    }
1948                    aoffset += 8 * lda;
1949
1950                    do {
1951                        for (int it = 0; it < 8; ++it)
1952                            c_arr[it] = vec_xl(0, (vector unsigned char*)aoffsets[it]);
1953                        vector_permute_store(c_arr, 8, vecOffset);
1954                        for (int it = 0; it < 8; ++it)
1955                            aoffsets[it] = aoffsets[it] + 8;
1956                        vecOffset += 128;
1957                        i--;
1958                    } while(i > 0);
1959                }
1960                j--;
1961            } while(j > 0);
1962        }
1963        if (rows & 4) {
1964            aoffsets[0] = aoffset;
1965            for (int it = 1; it < 4; ++it)
1966                aoffsets[it] = aoffsets[it-1] + lda;
1967            aoffset += 4 * lda;
1968            if (cols == 4) {
1969                for (int it = 0; it < 4; ++it)
1970                    c_arr[it] = vec_xl(0, (vector unsigned char*)aoffsets[it]);
1971                vector_permute_store(c_arr, 2, vecOffset);
1972                for (int it = 0; it< 4; it++)
1973                    aoffsets[it] = aoffsets[it] + lda;
1974                vecOffset += 32;
1975            }
1976            i = (cols >> 3);
1977            if (i > 0) {
1978                do {
1979                    for (int it = 0; it < 4; ++it)
1980                        c_arr[it] = vec_xl(0, (vector unsigned char*)aoffsets[it]);
1981                    vector_permute_store(c_arr, 4, vecOffset);
1982                    for (int it = 0; it< 4; it++)
1983                        aoffsets[it] = aoffsets[it] + 8*lda;
1984                    vecOffset += 64;
1985                    i--;
1986                } while(i > 0);
1987            }
1988        }
1989        if (rows & 3) {
1990            aoffsets[0] = aoffset;
1991            for (int it = 1; it < 4; ++it)
1992                aoffsets[it] = aoffsets[it-1] + lda;
1993            if (cols == 4) {
1994                switch(rows) {
1995                    case 3: c_arr[2] = vec_xl(0, (vector unsigned char*)aoffsets[2]);
1996                    case 2: c_arr[1] = vec_xl(0, (vector unsigned char*)aoffsets[1]);
1997                    case 1: c_arr[0] = vec_xl(0, (vector unsigned char*)aoffsets[0]);
1998                        break;
1999                }
2000                vector_permute_store(c_arr, 2, vecOffset);
2001                for (int it = 0; it< 4; it++)
2002                     aoffsets[it] = aoffsets[it] + lda;
2003                vecOffset += 32;
2004            }
2005            i = (cols >> 3);
2006            if (i > 0) {
2007                do {
2008                    switch(rows) {
2009                        case 3: c_arr[2] = vec_xl(0, (vector unsigned char*)aoffsets[2]);
2010                        case 2: c_arr[1] = vec_xl(0, (vector unsigned char*)aoffsets[1]);
2011                        case 1: c_arr[0] = vec_xl(0, (vector unsigned char*)aoffsets[0]);
2012                            break;
2013                    }
2014                    vector_permute_store(c_arr, 4, vecOffset);
2015                    for (int it = 0; it <4; it++)
2016                         aoffsets[it] = aoffsets[it] + 8* lda;
2017                    vecOffset += 64;
2018                    i--;
2019                } while(i > 0);
2020            }
2021        }
2022    }
2023
2024    void mnpack(int64_t m0, int64_t m, int64_t n0, int64_t n) {
2025        int64_t mc, nc, mp, np;
2026        int m_rem = MIN(m - m0, 8);
2027        int n_rem = MIN(n - n0, 8);
2028
2029        if (m_rem >= 8 && n_rem >= 8) {
2030            mc = 8;
2031            nc = 8;
2032            gemm<8,8>(m0, m, n0, n);
2033        } else if (m_rem >= 4 && n_rem >= 8) {
2034            mc = 4;
2035            nc = 8;
2036            gemm<4,8>(m0, m, n0, n);
2037        } else if (m_rem >=8 && n_rem >=4){
2038                mc = 8;
2039                nc = 4;
2040                gemm<8,4>(m0, m, n0, n);
2041        } else if ((m_rem < 4) && (n_rem >= 8)) {
2042            nc = 8;
2043            switch(m_rem) {
2044                case 1:
2045                    mc = 1;
2046                    gemm_Mx8<1>(m0, m, n0, n);
2047                    break;
2048                case 2:
2049                    mc = 2;
2050                    gemm_Mx8<2>(m0, m, n0, n);
2051                    break;
2052                case 3:
2053                    mc = 3;
2054                    gemm_Mx8<3>(m0, m, n0, n);
2055                    break;
2056                default:
2057                    return;
2058            }
2059        } else if (m_rem >= 4 && n_rem >= 4) {
2060            mc = 4;
2061            nc = 4;
2062            gemm_small<4, 4>(m0, m, n0, n);
2063        } else if ((m_rem > 4) && (n_rem < 4)) {
2064            mc = 4;
2065            switch(n_rem) {
2066                case 1:
2067                    nc = 1;
2068                    gemm_small<4, 1>(m0, m, n0, n);
2069                    break;
2070                case 2:
2071                    nc = 2;
2072                    gemm_small<4, 2>(m0, m, n0, n);
2073                    break;
2074                case 3:
2075                    nc = 3;
2076                    gemm_small<4, 3>(m0, m, n0, n);
2077                    break;
2078
2079                default:
2080                    return;
2081            }
2082        } else {
2083            switch((m_rem << 4) | n_rem) {
2084                case 0x43:
2085                    mc = 4;
2086                    nc = 3;
2087                    gemm_small<4, 3>(m0, m, n0, n);
2088                    break;
2089                case 0x42:
2090                    mc = 4;
2091                    nc = 2;
2092                    gemm_small<4, 2>(m0, m, n0, n);
2093                    break;
2094                case 0x41:
2095                    mc = 4;
2096                    nc = 1;
2097                    gemm_small<4, 1>(m0, m, n0, n);
2098                    break;
2099                case 0x34:
2100                    mc = 3;
2101                    nc = 4;
2102                    gemm_small<3, 4>(m0, m, n0, n);
2103                    break;
2104                case 0x33:
2105                    mc = 3;
2106                    nc = 3;
2107                    gemm_small<3, 3>(m0, m, n0, n);
2108                    break;
2109                case 0x32:
2110                    mc = 3;
2111                    nc = 2;
2112                    gemm_small<3, 2>(m0, m, n0, n);
2113                    break;
2114                case 0x31:
2115                    mc = 3;
2116                    nc = 1;
2117                    gemm_small<3, 1>(m0, m, n0, n);
2118                    break;
2119                case 0x24:
2120                    mc = 2;
2121                    nc = 4;
2122                    gemm_small<2,4>(m0, m, n0, n);
2123                    break;
2124                case 0x23:
2125                    mc = 2;
2126                    nc = 3;
2127                    gemm_small<2, 3>(m0, m, n0, n);
2128                    break;
2129                case 0x22:
2130                    mc = 2;
2131                    nc = 2;
2132                    gemm_small<2, 2>(m0, m, n0, n);
2133                    break;
2134                case 0x21:
2135                    mc = 2;
2136                    nc = 1;
2137                    gemm_small<2, 1>(m0, m, n0, n);
2138                    break;
2139                case 0x14:
2140                    mc = 1;
2141                    nc = 4;
2142                    gemm_small<1, 4>(m0, m, n0, n);
2143                    break;
2144                case 0x13:
2145                    mc = 1;
2146                    nc = 3;
2147                    gemm_small<1, 3>(m0, m, n0, n);
2148                    break;
2149                case 0x12:
2150                    mc = 1;
2151                    nc = 2;
2152                    gemm_small<1, 2>(m0, m, n0, n);
2153                    break;
2154                case 0x11:
2155                    mc = 1;
2156                    nc = 1;
2157                    gemm_small<1, 1>(m0, m, n0, n);
2158                    break;
2159                default:
2160                    return;
2161            }
2162        }
2163        mp = m0 + (m - m0) / mc * mc;
2164        np = n0 + (n - n0) / nc * nc;
2165        mnpack(mp, m, n0, np);
2166        mnpack(m0, m, np, n);
2167    }
2168
2169    void KERNEL_4x8(int64_t ii, int64_t jj) {
2170        vec_t vec_A[4], vec_B[8] , vec_C[4];
2171        acc_t acc_0, acc_1;
2172        __builtin_mma_xxsetaccz(&acc_0);
2173        __builtin_mma_xxsetaccz(&acc_1);
2174        for (int l = 0; l < k; l+=8) {
2175            packNormal((A+(ii*lda)+l), lda, 4, 8, (uint8_t*)vec_A);
2176            packNormal((B+(jj*ldb)+l), ldb, 8, 8, (uint8_t*)vec_B);
2177            for (int x = 0; x < 4; x++) {
2178                mma_instr<TA>::outer_product(&acc_0, vec_A[x], vec_B[x]);
2179                mma_instr<TA>::outer_product(&acc_1, vec_A[x], vec_B[x+4]);
2180            }
2181        }
2182        save_acc(&acc_0, ii, jj);
2183        save_acc(&acc_1, ii, jj+4);
2184    }
2185
2186    void KERNEL_8x4(int64_t ii, int64_t jj) {
2187        vec_t vec_A[8], vec_B[4] , vec_C[4];
2188        acc_t acc_0, acc_1;
2189        __builtin_mma_xxsetaccz(&acc_0);
2190        __builtin_mma_xxsetaccz(&acc_1);
2191        for (int l = 0; l < k; l+=8) {
2192            packNormal((A+(ii*lda)+l), lda, 8, 8, (uint8_t*)vec_A);
2193            packNormal((B+(jj*ldb)+l), ldb, 8, 4, (uint8_t*)vec_B);
2194            for (int x = 0; x < 4; x++) {
2195                mma_instr<TA>::outer_product(&acc_0, vec_A[x], vec_B[x]);
2196                mma_instr<TA>::outer_product(&acc_1, vec_A[x+4], vec_B[x]);
2197            }
2198        }
2199        save_acc(&acc_0, ii, jj);
2200        save_acc(&acc_1, ii+4, jj);
2201    }
2202
2203
2204    void KERNEL_8x8(int64_t ii, int64_t jj) {
2205        vec_t vec_A[8], vec_B[8], vec_C[4];
2206        acc_t acc_0, acc_1, acc_2, acc_3;
2207        __builtin_mma_xxsetaccz(&acc_0);
2208        __builtin_mma_xxsetaccz(&acc_1);
2209        __builtin_mma_xxsetaccz(&acc_2);
2210        __builtin_mma_xxsetaccz(&acc_3);
2211        for (int l = 0; l < k; l+=8) {
2212            packNormal(A+(ii*lda)+l, lda, 8, 8, (uint8_t*)vec_A);
2213            packNormal(B+(jj*ldb)+l, ldb, 8, 8, (uint8_t*)vec_B);
2214            for (int x = 0; x < 4; x++) {
2215                mma_instr<TA>::outer_product(&acc_0, vec_A[x], vec_B[x]);
2216                mma_instr<TA>::outer_product(&acc_1, vec_A[x], vec_B[x+4]);
2217                mma_instr<TA>::outer_product(&acc_2, vec_A[x+4], vec_B[x]);
2218                mma_instr<TA>::outer_product(&acc_3, vec_A[x+4], vec_B[x+4]);
2219            }
2220        }
2221        save_acc(&acc_0, ii, jj);
2222        save_acc(&acc_1, ii, jj+4);
2223        save_acc(&acc_2, ii+4, jj);
2224        save_acc(&acc_3, ii+4, jj+4);
2225    }
2226
2227    inline void MMA_16x8(vec_t * vec_A0, vec_t * vec_A1, vec_t * vec_B, acc_t * acc) {
2228        for (int x = 0; x < 4; x ++) {
2229             mma_instr<TA>::outer_product(&acc[0], vec_A0[x], vec_B[x]);
2230             mma_instr<TA>::outer_product(&acc[1], vec_A0[x], vec_B[x+4]);
2231             mma_instr<TA>::outer_product(&acc[2], vec_A0[x+4], vec_B[x]);
2232             mma_instr<TA>::outer_product(&acc[3], vec_A0[x+4], vec_B[x+4]);
2233             mma_instr<TA>::outer_product(&acc[4], vec_A1[x], vec_B[x]);
2234             mma_instr<TA>::outer_product(&acc[5], vec_A1[x], vec_B[x+4]);
2235             mma_instr<TA>::outer_product(&acc[6], vec_A1[x+4], vec_B[x]);
2236             mma_instr<TA>::outer_product(&acc[7], vec_A1[x+4], vec_B[x+4]);
2237        }
2238    }
2239    void KERNEL(int64_t ii, int64_t jj, int64_t mc, int64_t nc, int64_t kc, vec_t * vec_A, vec_t * vec_B, int64_t kk) {
2240        for (int64_t i = 0; i < mc; i += 16) {
2241            int A_base_addr = (mc / 8) * (i / 8) * 8;
2242            for (int64_t j = 0; j < nc; j += 8) {
2243                 int B_base_addr = (nc / 8) * (j / 8) * 8;
2244                 acc_t acc[8];
2245                 vec_t A0_block[8]; vec_t A1_block[8];
2246                 for (int x = 0; x < 8; x++)
2247                     __builtin_mma_xxsetaccz(&acc[x]);
2248                 for (int64_t l = 0; l < kc; l += 8) {
2249                     int A0_block_idx = A_base_addr + (l / 8) * 8;
2250                     int A1_block_idx = A0_block_idx + (mc / 8) * 8;
2251                     int B_block_idx = B_base_addr + (l / 8) * 8;
2252                     vec_t* A0_block = &vec_A[A0_block_idx];
2253                     vec_t* A1_block = &vec_A[A1_block_idx];
2254                     vec_t* B_block = &vec_B[B_block_idx];
2255                     MMA_16x8(A0_block, A1_block, B_block, acc);
2256                 }
2257                 if (kk == 0) {
2258                     save_acc(&acc[0], ii + i, jj + j);
2259                     save_acc(&acc[1], ii + i, jj + j + 4);
2260                     save_acc(&acc[2], ii + i + 4, jj + j);
2261                     save_acc(&acc[3], ii + i + 4, jj + j + 4);
2262                     save_acc(&acc[4], ii + i + 8, jj + j);
2263                     save_acc(&acc[5], ii + i + 8, jj + j + 4);
2264                     save_acc(&acc[6], ii + i + 12, jj + j);
2265                     save_acc(&acc[7], ii + i + 12, jj + j + 4);
2266                 } else {
2267                     add_save_acc(&acc[0], ii + i, jj + j);
2268                     add_save_acc(&acc[1], ii + i, jj + j + 4);
2269                     add_save_acc(&acc[2], ii + i + 4, jj + j);
2270                     add_save_acc(&acc[3], ii + i + 4, jj + j + 4);
2271                     add_save_acc(&acc[4], ii + i + 8, jj + j);
2272                     add_save_acc(&acc[5], ii + i + 8, jj + j + 4);
2273                     add_save_acc(&acc[6], ii + i + 12, jj + j);
2274                     add_save_acc(&acc[7], ii + i + 12, jj + j + 4);
2275                 }
2276            }
2277        }
2278    }
2279    template<int RM, int RN>
2280    void gemm_small(int64_t m0, int64_t m, int64_t n0, int64_t n) {
2281        int64_t ytiles = (m - m0) / RM;
2282        int64_t xtiles = (n - n0) / RN;
2283        int64_t tiles = xtiles * ytiles;
2284        int64_t duty = (tiles + nth - 1) / nth;
2285        int64_t start = duty * ith;
2286        int64_t end = start + duty;
2287        if (end > tiles)
2288            end = tiles;
2289        for (int64_t job = start; job < end; ++job) {
2290            int64_t ii = m0 + job / xtiles * RM;
2291            int64_t jj = n0 + job % xtiles * RN;
2292            vec_t vec_C[4];
2293            acc_t acc_0;
2294            __builtin_mma_xxsetaccz(&acc_0);
2295            vec_t vec_A[2], vec_B[2];
2296            for (int l=0; l<k; l+=4) {
2297                packNormal(A+(ii*lda)+l, lda, RM, 4, (uint8_t*)vec_A);
2298                packNormal(B+(jj*ldb)+l, ldb, RN, 4, (uint8_t*)vec_B);
2299                for (int x = 0; x<2; x++) {
2300                    mma_instr<TA>::outer_product(&acc_0, vec_A[x], vec_B[x]);
2301                }
2302            }
2303            __builtin_mma_disassemble_acc(vec_C, &acc_0);
2304            for (int I = 0; I < RM; I++) {
2305                for (int J = 0; J < RN; J++) {
2306                    *((TC*)(C+ii+((jj+J)*ldc)+I)) = *((TC*)&vec_C[I]+J);
2307                }
2308            }
2309        }
2310    }
2311
2312    template<int RM>
2313    void gemm_Mx8(int64_t m0, int64_t m, int64_t n0, int64_t n) {
2314        int RN = 8;
2315        int64_t ytiles = (m - m0) / RM;
2316        int64_t xtiles = (n - n0) / RN;
2317        int64_t tiles = xtiles * ytiles;
2318        int64_t duty = (tiles + nth - 1) / nth;
2319        int64_t start = duty * ith;
2320        int64_t end = start + duty;
2321        if (end > tiles)
2322            end = tiles;
2323        for (int64_t job = start; job < end; ++job) {
2324            int64_t ii = m0 + job / xtiles * RM;
2325            int64_t jj = n0 + job % xtiles * RN;
2326            vec_t vec_C[4];
2327            acc_t acc_0, acc_1;
2328            __builtin_mma_xxsetaccz(&acc_0);
2329            __builtin_mma_xxsetaccz(&acc_1);
2330            vec_t vec_A[4], vec_B[8];
2331            for (int l=0; l<k; l+=8) {
2332                packNormal(A+(ii*lda)+l, lda, RM, 8, (uint8_t*)vec_A);
2333                packNormal(B+(jj*ldb)+l, ldb, RN, 8, (uint8_t*)vec_B);
2334                for (int x = 0; x<4; x++) {
2335                    mma_instr<TA>::outer_product(&acc_0, vec_A[x], vec_B[x]);
2336                    mma_instr<TA>::outer_product(&acc_1, vec_A[x], vec_B[x+4]);
2337                }
2338            }
2339            __builtin_mma_disassemble_acc(vec_C, &acc_0);
2340            for (int I = 0; I < RM; I++) {
2341                for (int J = 0; J < 4; J++) {
2342                    *((TC*)(C+ii+((jj+J)*ldc)+I)) = *((TC*)&vec_C[I]+J);
2343                }
2344            }
2345            __builtin_mma_disassemble_acc(vec_C, &acc_1);
2346            for (int I = 0; I < RM; I++) {
2347                for (int J = 0; J < 4; J++) {
2348                    *((TC*)(C+ii+((jj+4+J)*ldc)+I)) = *((TC*)&vec_C[I]+J);
2349                }
2350            }
2351        }
2352    }
2353
2354    template<int RM, int RN>
2355    inline void kernel(int64_t ii, int64_t jj) {
2356       if constexpr(RM == 4 && RN == 8) {
2357          KERNEL_4x8(ii,jj);
2358       } else if constexpr(RM == 8 && RN == 8) {
2359          KERNEL_8x8(ii,jj);
2360       } else if constexpr(RM == 8 && RN == 4) {
2361          KERNEL_8x4(ii,jj);
2362       } else {
2363          assert(false && "RN/RM values not supported");
2364       }
2365    }
2366
2367    void matmul_tiled(int64_t m, int64_t n, int64_t mc, int64_t nc, int64_t kc) {
2368        int64_t ytiles = m / mc;
2369        int64_t xtiles = n / nc;
2370        int64_t tiles = xtiles * ytiles;
2371        int64_t duty = (tiles + nth - 1) / nth;
2372        int64_t start = duty * ith;
2373        int64_t end = start + duty;
2374        if (end > tiles) {
2375            end = tiles;
2376        }
2377        for (int64_t job = start; job < end; ++job) {
2378            int64_t ii = (job / xtiles) * mc;
2379            int64_t jj = (job % xtiles) * nc;
2380            for (int64_t kk = 0; kk < k; kk += kc) {
2381                 vec_t A_pack[kc * mc / 8];
2382                 vec_t B_pack[kc * nc / 8];
2383                 packNormal(A + (ii * lda) + kk, lda, kc, mc, (uint8_t *)A_pack);
2384                 packNormal(B + (jj * ldb) + kk, ldb, kc, nc, (uint8_t *)B_pack);
2385                 KERNEL(ii, jj, mc, nc, kc, A_pack, B_pack, kk);
2386            }
2387        }
2388    }
2389
2390    template <int RM, int RN>
2391    NOINLINE void gemm(int64_t m0, int64_t m, int64_t n0, int64_t n) {
2392        int64_t ytiles = (m - m0) / RM;
2393        int64_t xtiles = (n - n0) / RN;
2394        int64_t tiles = xtiles * ytiles;
2395        int64_t duty = (tiles + nth - 1) / nth;
2396        int64_t start = duty * ith;
2397        int64_t end = start + duty;
2398        if (end > tiles)
2399            end = tiles;
2400        for (int64_t job = start; job < end; ++job) {
2401            int64_t ii = m0 + job / xtiles * RM;
2402            int64_t jj = n0 + job % xtiles * RN;
2403            kernel<RM, RN>(ii, jj);
2404        }
2405    }
2406
2407    const TA *const A;
2408    const TB *const B;
2409    TC *C;
2410    const int64_t k;
2411    const int64_t lda;
2412    const int64_t ldb;
2413    const int64_t ldc;
2414    const int ith;
2415    const int nth;
2416};
2417
2418template <typename TA>
2419class tinyBLAS_Q0_PPC {
2420  public:
2421    tinyBLAS_Q0_PPC(int64_t k,
2422             const TA * A, int64_t lda,
2423             const block_q8_0 * B, int64_t ldb,
2424             float * C, int64_t ldc,
2425             int ith, int nth)
2426        : A(A), B(B), C(C), k(k), lda(lda), ldb(ldb), ldc(ldc), ith(ith), nth(nth) {
2427    }
2428
2429    void matmul(int64_t m, int64_t n) {
2430        int64_t mc = 64;
2431        int64_t nc = 64;
2432        int64_t kc = 64;
2433        int64_t n_chunk = 64;
2434    #if defined(_AIX) || defined(__BIG_ENDIAN__)
2435        mc = 32;
2436        nc = 32;
2437        kc = 32;
2438        n_chunk = 32;
2439    #endif
2440        int64_t n_aligned = 0;
2441        if (n % n_chunk == 0) {
2442            n_aligned = n;
2443        } else if (n == 4) {
2444            n_aligned = 4;
2445        } else if (n < n_chunk) {
2446            n_aligned = (n / 8) * 8;
2447        } else {
2448            n_aligned = (n / n_chunk) * n_chunk;
2449        }
2450        if (n_aligned > 0) {
2451            if (n_aligned % n_chunk == 0) nc = n_chunk;
2452            else if (n_aligned == n)      nc = n;
2453            else if (n_aligned % 32 == 0) nc = 32;
2454            else if (n_aligned % 24 == 0) nc = 24;
2455            else if (n_aligned % 16 == 0) nc = 16;
2456            else                          nc = 8;
2457        }
2458        bool can_use_tiled = n_aligned > 0 && (m % mc == 0);
2459        if (can_use_tiled) {
2460            matmul_tiled(m, n_aligned, mc, nc, kc);
2461            if (n > n_aligned) {
2462                mnpack(0, m, n_aligned, n);
2463            }
2464        } else {
2465            mnpack(0, m, 0, n);
2466        }
2467    }
2468
2469  private:
2470    inline void save_res(int ii, int jj, int idx, vector float * fin_res, int RM = 4, int RN = 4) {
2471        for (int I = 0; I < RM; I++) {
2472            for (int J = 0; J < RN; J++) {
2473                *((float *)(C + ii + ((jj + J) * ldc) + I)) = *((float *)&fin_res[idx + I] + J);
2474            }
2475        }
2476    }
2477
2478    inline void save_acc(acc_t * ACC, int64_t ii, int64_t jj) {
2479        vec_t vec_C[4];
2480        __builtin_mma_disassemble_acc(vec_C, ACC);
2481        for (int I = 0; I < 4; I++) {
2482            for (int J = 0; J < 4; J++) {
2483                *((float *)(C + ii + ((jj + J) * ldc) + I)) = *((float *)&vec_C[I] + J);
2484            }
2485        }
2486    }
2487
2488    inline void add_save_acc(acc_t * ACC, int64_t ii, int64_t jj) {
2489        vec_t vec_C[4];
2490        __builtin_mma_disassemble_acc(vec_C, ACC);
2491        for (int I = 0; I < 4; I++) {
2492            for (int J = 0; J < 4; J++) {
2493                float * c_ptr = (float *)(C + ii+ ((jj + J) * ldc) + I);
2494                *c_ptr += *((float *)&vec_C[I] + J);
2495            }
2496        }
2497    }
2498
2499    template<typename ArrayType>
2500    inline void compute(acc_t * ACC, int c_idx, int s_idx, ArrayType & comparray, vector float * vs, vector float * fin_res) {
2501        vector signed int vec_C[4];
2502        vector float CA[4] = {0};
2503        vector float res[4] = {0};
2504        __builtin_mma_disassemble_acc(vec_C, ACC);
2505        for (int i = 0; i < 4; i++) {
2506            CA[i] = vec_splats((float)(((double)comparray[c_idx + i]) * -128.0));
2507            res[i] = vec_add(vec_ctf(vec_C[i], 0), CA[i]);
2508            fin_res[s_idx + i] = vec_madd(res[i], vs[s_idx + i], fin_res[s_idx + i]);
2509        }
2510    }
2511
2512    inline void process_q4_elements(vector signed char (&c)[2], int * ca) {
2513        const vector signed char lowMask = vec_splats((signed char)0xF);
2514        const vector unsigned char v4 = vec_splats((unsigned char)0x4);
2515        const vector signed char v8 = vec_splats((signed char)0x8);
2516        vector signed int vsum = {0};
2517        vector signed int vsum2 = {0};
2518        c[0] = vec_and(c[1], lowMask);
2519        c[1] = vec_sr(c[1], v4);
2520        c[0] = vec_sub(c[0], v8);
2521        c[1] = vec_sub(c[1], v8);
2522        vsum = vec_sum4s(c[0], vsum);
2523        vsum2 = vec_sum4s(c[1], vsum2);
2524        vsum = vec_add(vsum, vsum2);
2525        *(ca) = vsum[0] + vsum[1] + vsum[2] + vsum[3];
2526    }
2527
2528    template <typename V1, typename V2>
2529    inline void vector_permute_store(V2 & s1, V2 & s2, V2 & s3, V2 & s4, V1 * vecOffset, bool flip) {
2530        vector unsigned char swiz1 = {0, 1, 2, 3, 4, 5, 6, 7, 16, 17, 18, 19, 20, 21, 22, 23};
2531        vector unsigned char swiz2 = {8, 9, 10, 11, 12, 13, 14, 15, 24, 25, 26, 27, 28, 29, 30, 31};
2532        vector unsigned char swiz3 = {0, 1, 2, 3, 8, 9, 10, 11, 16, 17, 18, 19, 24, 25, 26, 27};
2533        vector unsigned char swiz4 = {4, 5, 6, 7, 12, 13, 14, 15, 20, 21, 22, 23, 28, 29, 30, 31};
2534        V2 t1, t2, t3, t4, t5, t6, t7, t8;
2535        vector unsigned char xor_vector;
2536        uint8_t flip_vec = 0x80;
2537        xor_vector = vec_splats(flip_vec);
2538        t1 = vec_perm(s1, s2, swiz1);
2539        t2 = vec_perm(s1, s2, swiz2);
2540        t3 = vec_perm(s3, s4, swiz1);
2541        t4 = vec_perm(s3, s4, swiz2);
2542        t5 = vec_perm(t1, t3, swiz3);
2543        t6 = vec_perm(t1, t3, swiz4);
2544        t7 = vec_perm(t2, t4, swiz3);
2545        t8 = vec_perm(t2, t4, swiz4);
2546        if (flip == true) {
2547            t5 = vec_xor(t5, xor_vector);
2548            t6 = vec_xor(t6, xor_vector);
2549            t7 = vec_xor(t7, xor_vector);
2550            t8 = vec_xor(t8, xor_vector);
2551        }
2552        vec_xst(t5, 0, vecOffset);
2553        vec_xst(t6, 0, vecOffset + 16);
2554        vec_xst(t7, 0, vecOffset + 32);
2555        vec_xst(t8, 0, vecOffset + 48);
2556    }
2557
2558    inline void unpack_q4_to_q8(vector signed char packed, vector signed char & lo, vector signed char & hi) {
2559        const vector signed char lowMask = vec_splats((signed char)0x0F);
2560        const vector signed char v8      = vec_splats((signed char)0x08);
2561        const vector unsigned char v4    = vec_splats((unsigned char)4);
2562        lo = vec_and(packed, lowMask);
2563        hi = vec_sr(packed, v4);
2564        lo = vec_sub(lo, v8);
2565        hi = vec_sub(hi, v8);
2566    }
2567
2568    inline void vector_permute_store_fp16(vec_t * c, unsigned char * vecOffset) {
2569        vec_t t[8], s[8];
2570        vec_t swiz1 = {0, 1, 2, 3, 16, 17, 18, 19, 4, 5, 6, 7, 20, 21, 22, 23};
2571        vec_t swiz2 = {8, 9, 10, 11, 24, 25, 26, 27, 12, 13, 14, 15, 28, 29, 30, 31};
2572        vec_t swiz3 = {0, 1, 2, 3, 4, 5, 6, 7, 16, 17, 18, 19, 20, 21, 22, 23};
2573        vec_t swiz4 = {8, 9, 10, 11, 12, 13, 14, 15, 24, 25, 26, 27, 28, 29, 30, 31};
2574        for (int i = 0; i < 4; i += 2) {
2575            t[i + 0] = vec_perm(c[i + 0], c[i + 1], swiz1);
2576            t[i + 1] = vec_perm(c[i + 0], c[i + 1], swiz2);
2577        }
2578        for (int i = 4; i < 8; i += 2) {
2579            t[i + 0] = vec_perm(c[i + 0], c[i + 1], swiz1);
2580            t[i + 1] = vec_perm(c[i + 0], c[i + 1], swiz2);
2581        }
2582        s[0] = vec_perm(t[0], t[2], swiz3);
2583        s[1] = vec_perm(t[0], t[2], swiz4);
2584        s[2] = vec_perm(t[1], t[3], swiz3);
2585        s[3] = vec_perm(t[1], t[3], swiz4);
2586        s[4] = vec_perm(t[4], t[6], swiz3);
2587        s[5] = vec_perm(t[4], t[6], swiz4);
2588        s[6] = vec_perm(t[5], t[7], swiz3);
2589        s[7] = vec_perm(t[5], t[7], swiz4);
2590        for (int i = 0; i < 8; ++i) {
2591            vec_xst(s[i], 0, (vec_t *)(vecOffset + i * 16));
2592        }
2593    }
2594
2595    static inline void convert_and_scale_q8(vector signed char raw, vector float v_scale, vector unsigned short & out_hi, vector unsigned short & out_lo) {
2596        vector signed short i16_hi = vec_unpackh(raw);
2597        vector signed short i16_lo = vec_unpackl(raw);
2598
2599        vector float f_hi_h = vec_ctf(vec_unpackh(i16_hi), 0);
2600        vector float f_hi_l = vec_ctf(vec_unpackl(i16_hi), 0);
2601        vector float f_lo_h = vec_ctf(vec_unpackh(i16_lo), 0);
2602        vector float f_lo_l = vec_ctf(vec_unpackl(i16_lo), 0);
2603        out_hi = vec_pack_to_short_fp32(vec_mul(f_hi_h, v_scale), vec_mul(f_hi_l, v_scale));
2604        out_lo = vec_pack_to_short_fp32(vec_mul(f_lo_h, v_scale), vec_mul(f_lo_l, v_scale));
2605    }
2606
2607    void packNormal_q4_fp16(const block_q4_0 * a, int64_t lda, int rows, int blocks, unsigned char * vec) {
2608        unsigned char * vecOffset = vec;
2609        for (int i = 0; i < rows; i += 8) {
2610            const block_q4_0 * rows_base[8];
2611            for (int r = 0; r < 8; r++) {
2612                rows_base[r] = a + (i + r) * lda;
2613            }
2614            for (int blk = 0; blk < blocks; blk++) {
2615                vector unsigned short hp_res[8][4];
2616                for (int r = 0; r < 8; r++) {
2617                    const block_q4_0 * current_blk = rows_base[r] + blk;
2618                    vector float v_scale = vec_extract_fp32_from_shorth(vec_splats(current_blk->d));
2619                    vector signed char v_qs = vec_xl(0, (const vector signed char *)current_blk->qs);
2620                    vector signed char c1, c2;
2621                    unpack_q4_to_q8(v_qs, c1, c2);
2622                    convert_and_scale_q8(c1, v_scale, hp_res[r][0], hp_res[r][1]);
2623                    convert_and_scale_q8(c2, v_scale, hp_res[r][2], hp_res[r][3]);
2624                }
2625                for (int c = 0; c < 4; c++) {
2626                    vector unsigned char c_arr[8];
2627                    for (int r = 0; r < 8; r++) {
2628                        c_arr[r] = (vector unsigned char)hp_res[r][c];
2629                    }
2630                    vector_permute_store_fp16((vec_t *)c_arr, vecOffset);
2631                    vecOffset += 128;
2632                }
2633            }
2634        }
2635    }
2636
2637    template <int chunk_size>
2638    static inline void pack_q8_block(const block_q8_0 * a, int64_t lda, int rows, int blocks, unsigned char * vec) {
2639        unsigned char * vecOffset = vec;
2640        const vec_t swiz1 = {0, 1, 2, 3, 16, 17, 18, 19, 4, 5, 6, 7, 20, 21, 22, 23};
2641        const vec_t swiz2 = {8, 9, 10, 11, 24, 25, 26, 27, 12, 13, 14, 15, 28, 29, 30, 31};
2642        const vec_t swiz3 = {0, 1, 2, 3, 4, 5, 6, 7, 16, 17, 18, 19, 20, 21, 22, 23};
2643        const vec_t swiz4 = {8, 9, 10, 11, 12, 13, 14, 15, 24, 25, 26, 27, 28, 29, 30, 31};
2644
2645        for (int i = 0; i < rows; i += chunk_size) {
2646            const block_q8_0 * rows_base[chunk_size];
2647            for (int r = 0; r < chunk_size; r++) {
2648                rows_base[r] = a + (i + r) * lda;
2649            }
2650            for (int blk = 0; blk < blocks; blk++) {
2651                vector unsigned short hp_res[chunk_size][4];
2652                for (int r = 0; r < chunk_size; r++) {
2653                    const block_q8_0 * b = rows_base[r] + blk;
2654                    vector float v_scale = vec_extract_fp32_from_shorth(vec_splats(b->d));
2655                    vector signed char c[2];
2656                    __vector_pair pair = __builtin_vsx_lxvp(0, (__vector_pair *)b->qs);
2657                    __builtin_vsx_disassemble_pair(c, & pair);
2658                    convert_and_scale_q8(c[0], v_scale, hp_res[r][0], hp_res[r][1]);
2659                    convert_and_scale_q8(c[1], v_scale, hp_res[r][2], hp_res[r][3]);
2660                }
2661                for (int col = 0; col < 4; col++) {
2662                    if constexpr (chunk_size == 8) {
2663                        vec_t t[8];
2664                        t[0] = vec_perm((vec_t)hp_res[0][col], (vec_t)hp_res[1][col], swiz1);
2665                        t[1] = vec_perm((vec_t)hp_res[0][col], (vec_t)hp_res[1][col], swiz2);
2666                        t[2] = vec_perm((vec_t)hp_res[2][col], (vec_t)hp_res[3][col], swiz1);
2667                        t[3] = vec_perm((vec_t)hp_res[2][col], (vec_t)hp_res[3][col], swiz2);
2668                        t[4] = vec_perm((vec_t)hp_res[4][col], (vec_t)hp_res[5][col], swiz1);
2669                        t[5] = vec_perm((vec_t)hp_res[4][col], (vec_t)hp_res[5][col], swiz2);
2670                        t[6] = vec_perm((vec_t)hp_res[6][col], (vec_t)hp_res[7][col], swiz1);
2671                        t[7] = vec_perm((vec_t)hp_res[6][col], (vec_t)hp_res[7][col], swiz2);
2672
2673                        vec_xst(vec_perm(t[0], t[2], swiz3), 0, (vec_t *)(vecOffset + 0));
2674                        vec_xst(vec_perm(t[0], t[2], swiz4), 0, (vec_t *)(vecOffset + 16));
2675                        vec_xst(vec_perm(t[1], t[3], swiz3), 0, (vec_t *)(vecOffset + 32));
2676                        vec_xst(vec_perm(t[1], t[3], swiz4), 0, (vec_t *)(vecOffset + 48));
2677                        vec_xst(vec_perm(t[4], t[6], swiz3), 0, (vec_t *)(vecOffset + 64));
2678                        vec_xst(vec_perm(t[4], t[6], swiz4), 0, (vec_t *)(vecOffset + 80));
2679                        vec_xst(vec_perm(t[5], t[7], swiz3), 0, (vec_t *)(vecOffset + 96));
2680                        vec_xst(vec_perm(t[5], t[7], swiz4), 0, (vec_t *)(vecOffset + 112));
2681                        vecOffset += 128;
2682                    } else {
2683                        vec_t t0 = vec_perm((vec_t)hp_res[0][col], (vec_t)hp_res[1][col], swiz1);
2684                        vec_t t1 = vec_perm((vec_t)hp_res[0][col], (vec_t)hp_res[1][col], swiz2);
2685                        vec_t t2 = vec_perm((vec_t)hp_res[2][col], (vec_t)hp_res[3][col], swiz1);
2686                        vec_t t3 = vec_perm((vec_t)hp_res[2][col], (vec_t)hp_res[3][col], swiz2);
2687
2688                        vec_xst(vec_perm(t0, t2, swiz3), 0, (vec_t *)(vecOffset + 0));
2689                        vec_xst(vec_perm(t0, t2, swiz4), 0, (vec_t *)(vecOffset + 16));
2690                        vec_xst(vec_perm(t1, t3, swiz3), 0, (vec_t *)(vecOffset + 32));
2691                        vec_xst(vec_perm(t1, t3, swiz4), 0, (vec_t *)(vecOffset + 48));
2692                        vecOffset += 64;
2693                    }
2694                }
2695            }
2696        }
2697    }
2698
2699    void packNormal_q8_fp16(const block_q8_0 * a, int64_t lda, int rows, int blocks, unsigned char * vec) {
2700        if (rows == 4) {
2701            pack_q8_block<4>(a, lda, rows, blocks, vec);
2702        } else {
2703            pack_q8_block<8>(a, lda, rows, blocks, vec);
2704        }
2705    }
2706
2707    template<int size>
2708    void packNormalInt4(const TA * a, int64_t lda, int rows, int cols, int8_t * vec, std::array<int, size> & comparray) {
2709        int64_t i, j;
2710        TA * aoffset = NULL;
2711        int8_t * vecOffset = NULL;
2712        TA * aoffset1 = NULL, * aoffset2 = NULL, * aoffset3 = NULL, * aoffset4 = NULL;
2713        TA * aoffset5 = NULL, * aoffset6 = NULL, * aoffset7 = NULL, * aoffset8 = NULL;
2714        vector signed char c1[2] = {0}, c2[2] = {0}, c3[2] = {0}, c4[2] = {0};
2715        vector signed char c5[2] = {0}, c6[2] = {0}, c7[2] = {0}, c8[2] = {0};
2716        aoffset = const_cast<TA *>(a);
2717        vecOffset = vec;
2718        j = (rows >> 3);
2719        if (j > 0) {
2720            do {
2721                aoffset1 = aoffset;
2722                aoffset2 = aoffset1 + lda;
2723                aoffset3 = aoffset2 + lda;
2724                aoffset4 = aoffset3 + lda;
2725                aoffset5 = aoffset4 + lda;
2726                aoffset6 = aoffset5 + lda;
2727                aoffset7 = aoffset6 + lda;
2728                aoffset8 = aoffset7 + lda;
2729                aoffset += 8 * lda;
2730                i = (cols >> 2);
2731                if (i > 0) {
2732                    do {
2733                        c1[1] = vec_xl(0, (const vector signed char *)aoffset1->qs);
2734                        c2[1] = vec_xl(0, (const vector signed char *)aoffset2->qs);
2735                        c3[1] = vec_xl(0, (const vector signed char *)aoffset3->qs);
2736                        c4[1] = vec_xl(0, (const vector signed char *)aoffset4->qs);
2737                        c5[1] = vec_xl(0, (const vector signed char *)aoffset5->qs);
2738                        c6[1] = vec_xl(0, (const vector signed char *)aoffset6->qs);
2739                        c7[1] = vec_xl(0, (const vector signed char *)aoffset7->qs);
2740                        c8[1] = vec_xl(0, (const vector signed char *)aoffset8->qs);
2741
2742                        process_q4_elements(c1, & comparray[0]);
2743                        process_q4_elements(c2, & comparray[1]);
2744                        process_q4_elements(c3, & comparray[2]);
2745                        process_q4_elements(c4, & comparray[3]);
2746                        process_q4_elements(c5, & comparray[4]);
2747                        process_q4_elements(c6, & comparray[5]);
2748                        process_q4_elements(c7, & comparray[6]);
2749                        process_q4_elements(c8, & comparray[7]);
2750                        vector_permute_store<int8_t, vector signed char>(c1[0], c2[0], c3[0], c4[0], vecOffset, false);
2751                        vector_permute_store<int8_t, vector signed char>(c1[1], c2[1], c3[1], c4[1], vecOffset + 64, false);
2752                        vector_permute_store<int8_t, vector signed char>(c5[0], c6[0], c7[0], c8[0], vecOffset + 128, false);
2753                        vector_permute_store<int8_t, vector signed char>(c5[1], c6[1], c7[1], c8[1], vecOffset + 192, false);
2754                        aoffset1 += lda;
2755                        aoffset2 += lda;
2756                        aoffset3 += lda;
2757                        aoffset4 += lda;
2758                        aoffset5 += lda;
2759                        aoffset6 += lda;
2760                        aoffset7 += lda;
2761                        aoffset8 += lda;
2762                        vecOffset += 256;
2763                        i--;
2764                    } while (i > 0);
2765                }
2766                j--;
2767            } while (j > 0);
2768        }
2769
2770        if (rows & 4) {
2771            aoffset1 = aoffset;
2772            aoffset2 = aoffset1 + lda;
2773            aoffset3 = aoffset2 + lda;
2774            aoffset4 = aoffset3 + lda;
2775            aoffset += 4 * lda;
2776            i = (cols >> 2);
2777            if (i > 0) {
2778                do {
2779                    c1[1] = vec_xl(0, (const vector signed char *)aoffset1->qs);
2780                    c2[1] = vec_xl(0, (const vector signed char *)aoffset2->qs);
2781                    c3[1] = vec_xl(0, (const vector signed char *)aoffset3->qs);
2782                    c4[1] = vec_xl(0, (const vector signed char *)aoffset4->qs);
2783
2784                    process_q4_elements(c1, & comparray[0]);
2785                    process_q4_elements(c2, & comparray[1]);
2786                    process_q4_elements(c3, & comparray[2]);
2787                    process_q4_elements(c4, & comparray[3]);
2788                    vector_permute_store<int8_t, vector signed char>(c1[0], c2[0], c3[0], c4[0], vecOffset, false);
2789                    vector_permute_store<int8_t, vector signed char>(c1[1], c2[1], c3[1], c4[1], vecOffset + 64, false);
2790                    aoffset1 += lda;
2791                    aoffset2 += lda;
2792                    aoffset3 += lda;
2793                    aoffset4 += lda;
2794                    vecOffset += 128;
2795                    i--;
2796                } while (i > 0);
2797            }
2798        }
2799
2800        if (rows & 3) {
2801            aoffset1 = aoffset;
2802            aoffset2 = aoffset1 + lda;
2803            aoffset3 = aoffset2 + lda;
2804            i = (cols >> 2);
2805            if (i > 0) {
2806                do {
2807                    switch(rows) {
2808                        case 3: c3[1] = vec_xl(0, (const vector signed char *)aoffset3->qs);
2809                        case 2: c2[1] = vec_xl(0, (const vector signed char *)aoffset2->qs);
2810                        case 1: c1[1] = vec_xl(0, (const vector signed char *)aoffset1->qs);
2811                            break;
2812                    }
2813                    process_q4_elements(c1, & comparray[0]);
2814                    process_q4_elements(c2, & comparray[1]);
2815                    process_q4_elements(c3, & comparray[2]);
2816                    process_q4_elements(c4, & comparray[3]);
2817                    vector_permute_store<int8_t, vector signed char>(c1[0], c2[0], c3[0], c4[0], vecOffset, false);
2818                    vector_permute_store<int8_t, vector signed char>(c1[1], c2[1], c3[1], c4[1], vecOffset + 64, false);
2819                    aoffset1 += lda;
2820                    aoffset2 += lda;
2821                    aoffset3 += lda;
2822                    vecOffset += 128;
2823                    i--;
2824                } while(i > 0);
2825            }
2826        }
2827    }
2828
2829    template<typename VA, typename VB>
2830    void packNormal(const block_q8_0 * a, int64_t lda, int rows, int cols, VA * vec, bool flip) {
2831        int64_t i, j;
2832        block_q8_0 * aoffset = NULL;
2833        VA * vecOffset = NULL;
2834        block_q8_0 * aoffsets[8];
2835        __vector_pair arr[8];
2836        VB c[8][2] = {0};
2837        VB c1[8] = {0}; VB c2[8] = {0};
2838        aoffset = const_cast<block_q8_0 *>(a);
2839        vecOffset = vec;
2840        j = (rows >> 3);
2841        if (j > 0) {
2842            do {
2843                aoffsets[0] = aoffset;
2844                for (int it = 1; it < 8; it++)
2845                    aoffsets[it] = aoffsets[it - 1] + lda;
2846                aoffset += 8 * lda;
2847
2848                i = (cols >> 3);
2849                if (i > 0) {
2850                do {
2851                    for (int it = 0; it < 8; it++) {
2852                        arr[it] = __builtin_vsx_lxvp(0, (__vector_pair *)aoffsets[it]->qs);
2853                        __builtin_vsx_disassemble_pair(c[it], & arr[it]);
2854                        c1[it] = c[it][0];
2855                        c2[it] = c[it][1];
2856                    }
2857                    vector_permute_store<VA, VB>(c1[0], c1[1], c1[2], c1[3], vecOffset, flip);
2858                    vector_permute_store<VA, VB>(c2[0], c2[1], c2[2], c2[3], vecOffset + 64, flip);
2859                    vector_permute_store<VA, VB>(c1[4], c1[5], c1[6], c1[7], vecOffset + 128, flip);
2860                    vector_permute_store<VA, VB>(c2[4], c2[5], c2[6], c2[7], vecOffset + 192, flip);
2861                    for (int it = 0; it < 8; it++)
2862                        aoffsets[it] += lda;
2863                    vecOffset += 256;
2864                    i--;
2865               } while(i > 0);
2866            }
2867            j--;
2868        } while(j > 0);
2869    }
2870    if (rows & 4) {
2871            aoffsets[0]  = aoffset;
2872            for (int it = 1; it < 4; it++ )
2873                aoffsets[it] = aoffsets[it-1] + lda;
2874            aoffset += 4 * lda;
2875        i = (cols >> 3);
2876            if (i > 0) {
2877               do {
2878                    for (int it = 0; it < 4; it++) {
2879                        arr[it] = __builtin_vsx_lxvp(0, (__vector_pair *)aoffsets[it]->qs);
2880                        __builtin_vsx_disassemble_pair(c[it], & arr[it]);
2881                        c1[it] = c[it][0];
2882                        c2[it] = c[it][1];
2883                    }
2884                    vector_permute_store<VA, VB>(c1[0], c1[1], c1[2], c1[3], vecOffset, flip);
2885                    vector_permute_store<VA, VB>(c2[0], c2[1], c2[2], c2[3], vecOffset + 64, flip);
2886                    for (int it = 0; it < 4; it++) {
2887                        aoffsets[it] += lda;
2888                    }
2889                    vecOffset += 128;
2890                    i--;
2891               } while(i > 0);
2892            }
2893        }
2894
2895        if (rows & 3) {
2896            aoffsets[0]  = aoffset;
2897            for (int it = 1; it < 3; it++ )
2898                aoffsets[it] = aoffsets[it - 1] + lda;
2899            i = (cols >> 3);
2900            if (i > 0) {
2901                do {
2902                    switch(rows) {
2903                        case 3: arr[2] = __builtin_vsx_lxvp(0, (__vector_pair *)aoffsets[2]->qs);
2904                                __builtin_vsx_disassemble_pair(c[2], & arr[2]);
2905                                c1[2] = c[2][0]; c2[2] = c[2][1];
2906                        case 2: arr[1] = __builtin_vsx_lxvp(0, (__vector_pair *)aoffsets[1]->qs);
2907                                __builtin_vsx_disassemble_pair(c[1], & arr[1]);
2908                                c1[1] = c[1][0]; c2[1] = c[1][1];
2909                        case 1: arr[0] = __builtin_vsx_lxvp(0, (__vector_pair *)aoffsets[0]->qs);
2910                                __builtin_vsx_disassemble_pair(c[0], & arr[0]);
2911                                c1[0] = c[0][0]; c2[0] = c[0][1];
2912                                break;
2913                    }
2914                    vector_permute_store<VA, VB>(c1[0], c1[1], c1[2], c1[3], vecOffset, flip);
2915                    vector_permute_store<VA, VB>(c2[0], c2[1], c2[2], c2[3], vecOffset + 64, flip);
2916                    for (int it = 0; it < 3; it++)
2917                         aoffsets[it] += lda;
2918                    vecOffset += 128;
2919                    i--;
2920               } while(i > 0);
2921            }
2922        }
2923    }
2924
2925    void mnpack(int64_t m0, int64_t m, int64_t n0, int64_t n) {
2926        int m_rem = MIN(m - m0, 16);
2927        int n_rem = MIN(n - n0, 16);
2928
2929        int mc = 0, nc = 0;
2930
2931        if (m_rem >= 8 && n_rem >= 8) {
2932           mc = 8;
2933           nc = 8;
2934           gemm<8, 8>(m0, m, n0, n);
2935        } else if (m_rem >= 4 && n_rem >= 8) {
2936            mc = 4;
2937            nc = 8;
2938            gemm<4, 8>(m0, m, n0, n);
2939        } else if (m_rem >= 8 && n_rem >= 4) {
2940            mc = 8;
2941            nc = 4;
2942            gemm<8, 4>(m0, m, n0, n);
2943        } else if (m_rem >= 4 && n_rem >= 4) {
2944            mc = 4;
2945            nc = 4;
2946            gemm_small(m0, m, n0, n, mc, nc);
2947        } else {
2948            mc = (m_rem >= 4) ? 4 : m_rem;
2949            nc = (n_rem >= 4) ? 4 : n_rem;
2950            if (mc == 0 || nc == 0)
2951               return;
2952            gemm_small(m0, m, n0, n, mc, nc);
2953        }
2954
2955        int64_t mp = m0 + ((m - m0) / mc) * mc;
2956        int64_t np = n0 + ((n - n0) / nc) * nc;
2957        mnpack(mp, m, n0, np);
2958        mnpack(m0, m, np, n);
2959    }
2960
2961
2962    void KERNEL_4x8(int64_t ii, int64_t jj) {
2963        vec_t vec_A[8], vec_B[16] = {0};
2964        acc_t acc_0, acc_1;
2965        std::array<int, 4> comparray {};
2966        vector float fin_res[8] = {0};
2967        vector float vs[8] = {0};
2968        bool isAblock_q4 = std::is_same_v<TA, block_q4_0>;
2969        for (int l = 0; l < k; l++) {
2970            __builtin_mma_xxsetaccz(& acc_0);
2971            __builtin_mma_xxsetaccz(& acc_1);
2972            if (std::is_same_v<TA, block_q4_0>) {
2973               packNormalInt4<4>((A + (ii * lda) + l), lda, 4, 4, (int8_t *)vec_A, comparray);
2974            } else {
2975               packNormal<int8_t, vector signed char>((const block_q8_0 *)(A + (ii * lda) + l), lda, 4, 8, (int8_t *)vec_A, false);
2976            }
2977            packNormal<uint8_t, vector unsigned char>((B + (jj * ldb) + l), ldb, 8, 8, (uint8_t *)vec_B, true);
2978            for(int x = 0; x < 8; x++) {
2979                __builtin_mma_xvi8ger4pp(& acc_0, vec_A[x], vec_B[x]);
2980                __builtin_mma_xvi8ger4pp(& acc_1, vec_A[x], vec_B[x+8]);
2981            }
2982            for (int I = 0; I<4; I++) {
2983                for (int J = 0; J<4; J++) {
2984                    *((float *)& vs[I] + J) = (unhalf((A + ((ii + I) * lda) + l)->d) * unhalf((B + ((jj + J) * ldb) + l)->d));
2985                    *((float *)& vs[I + 4] + J) = (unhalf((A +((ii + I) * lda) + l)->d) * unhalf((B + ((jj + J + 4) * ldb) + l)->d));
2986                }
2987            }
2988            if (!isAblock_q4) {
2989                auto aoffset = A + (ii * lda) + l;
2990                for (int i = 0; i < 4; i++) {
2991                    comparray[i] = 0;
2992                    int ca = 0;
2993                    auto *at = aoffset->qs;
2994                    for (int j = 0; j < 32; j++)
2995                        ca += (int)*at++;
2996                    comparray[i] = ca;
2997                    aoffset += lda;
2998                }
2999            }
3000            compute(& acc_0, 0, 0, comparray, vs, fin_res);
3001            compute(& acc_1, 0, 4, comparray, vs, fin_res);
3002        }
3003        save_res(ii, jj, 0, fin_res);
3004        save_res(ii, jj + 4, 4, fin_res);
3005    }
3006
3007    void KERNEL_8x4(int64_t ii, int64_t jj) {
3008        vec_t vec_A[16], vec_B[8] = {0};
3009        acc_t acc_0, acc_1;
3010        std::array<int, 8> comparray {};
3011        vector float fin_res[8] = {0};
3012        vector float vs[8] = {0};
3013        bool isAblock_q4 = std::is_same_v<TA, block_q4_0>;
3014        for (int l = 0; l < k; l++) {
3015            __builtin_mma_xxsetaccz(& acc_0);
3016            __builtin_mma_xxsetaccz(& acc_1);
3017            if (std::is_same_v<TA, block_q4_0>) {
3018               packNormalInt4<8>((A + (ii * lda) + l), lda, 8, 4, (int8_t *)vec_A, comparray);
3019            } else {
3020               packNormal<int8_t, vector signed char>((const block_q8_0 *)(A + (ii * lda) + l), lda, 8, 8, (int8_t *)vec_A, false);
3021            }
3022            packNormal<uint8_t, vector unsigned char>((B + (jj * ldb) + l), ldb, 4, 8, (uint8_t *)vec_B, true);
3023            for(int x = 0; x < 8; x++) {
3024                __builtin_mma_xvi8ger4pp(& acc_0, vec_A[x], vec_B[x]);
3025                __builtin_mma_xvi8ger4pp(& acc_1, vec_A[x + 8], vec_B[x]);
3026            }
3027            for (int I = 0; I < 8; I++) {
3028                for (int J = 0; J < 4; J++) {
3029                    *((float *)&vs[I] + J) = (unhalf((A + ((ii + I) * lda) + l)->d) * unhalf((B + ((jj + J) * ldb) + l)->d));
3030                }
3031            }
3032            if (!isAblock_q4) {
3033                auto aoffset = A + (ii * lda) + l;
3034                for (int i = 0; i < 8; i++) {
3035                    comparray[i] = 0;
3036                    int ca = 0;
3037                    auto *at = aoffset->qs;
3038                    for (int j = 0; j < 32; j++)
3039                        ca += (int)*at++;
3040                    comparray[i] = ca;
3041                    aoffset += lda;
3042                }
3043            }
3044            compute(& acc_0, 0, 0, comparray, vs, fin_res);
3045            compute(& acc_1, 4, 4, comparray, vs, fin_res);
3046        }
3047        save_res(ii, jj, 0, fin_res);
3048        save_res(ii + 4, jj, 4, fin_res);
3049    }
3050
3051    void KERNEL_8x8(int64_t ii, int64_t jj) {
3052        vec_t vec_A[16], vec_B[16] = {0};
3053        acc_t acc_0, acc_1, acc_2, acc_3;
3054        acc_t acc_4, acc_5, acc_6, acc_7;
3055        std::array<int, 8> comparray {};
3056        vector float fin_res[16] = {0};
3057        vector float vs[16] = {0};
3058        bool isAblock_q4 = std::is_same_v<TA, block_q4_0>;
3059        for (int l = 0; l < k; l++) {
3060            __builtin_mma_xxsetaccz(& acc_0);
3061            __builtin_mma_xxsetaccz(& acc_1);
3062            __builtin_mma_xxsetaccz(& acc_2);
3063            __builtin_mma_xxsetaccz(& acc_3);
3064            if (std::is_same_v<TA, block_q4_0>) {
3065               packNormalInt4<8>((A + (ii * lda) + l), lda, 8, 4, (int8_t *)vec_A, comparray);
3066            } else {
3067               packNormal<int8_t, vector signed char>((const block_q8_0 *)(A + (ii * lda) + l), lda, 8, 8, (int8_t *)vec_A, false);
3068            }
3069            packNormal<uint8_t, vector unsigned char>((B + (jj * ldb) + l), ldb, 8, 8, (uint8_t *)vec_B, true);
3070            for(int x = 0; x < 8; x++) {
3071                __builtin_mma_xvi8ger4pp(& acc_0, vec_A[x], vec_B[x]);
3072                __builtin_mma_xvi8ger4pp(& acc_1, vec_A[x + 8], vec_B[x]);
3073                __builtin_mma_xvi8ger4pp(& acc_2, vec_A[x], vec_B[x + 8]);
3074                __builtin_mma_xvi8ger4pp(& acc_3, vec_A[x + 8], vec_B[x + 8]);
3075            }
3076            for (int I = 0; I < 8 ; I++) {
3077                for (int J = 0; J < 4; J++) {
3078                    *((float *)& vs[I] + J) = (unhalf((A + ((ii + I) * lda) + l)->d) * unhalf((B + ((jj + J) * ldb) + l)->d));
3079                    *((float *)& vs[I + 8] + J) = (unhalf((A + ((ii + I) * lda) + l)->d) * unhalf((B + ((jj + J + 4) * ldb) + l)->d));
3080                }
3081            }
3082            if (!isAblock_q4) {
3083                auto aoffset = A + (ii * lda) + l;
3084                for (int i = 0; i < 8; i++) {
3085                    comparray[i] = 0;
3086                    int ca = 0;
3087                    auto *at = aoffset->qs;
3088                    for (int j = 0; j < 32; j++)
3089                        ca += (int)*at++;
3090                    comparray[i] = ca;
3091                    aoffset += lda;
3092                }
3093            }
3094            compute(& acc_0, 0, 0, comparray, vs, fin_res);
3095            compute(& acc_1, 4, 4, comparray, vs, fin_res);
3096            compute(& acc_2, 0, 8, comparray, vs, fin_res);
3097            compute(& acc_3, 4, 12, comparray, vs, fin_res);
3098        }
3099        save_res(ii, jj, 0, fin_res);
3100        save_res(ii + 4, jj, 4, fin_res);
3101        save_res(ii, jj + 4, 8, fin_res);
3102        save_res(ii + 4, jj + 4, 12, fin_res);
3103    }
3104
3105    void KERNEL_Q0(int64_t ii, int64_t jj, int64_t mc, int64_t nc, int64_t kc, int64_t l, vec_t * vec_A, vec_t * vec_B) {
3106        acc_t acc[8];
3107        for (int i = 0; i < mc ; i += 16) {
3108            for (int j = 0; j < nc; j += 8) {
3109                int A0_base = (i / 16) * (2 * 32 * kc);
3110                int B0_base = (j / 8) * (32 * kc);
3111                for (int x = 0; x < 8; x++) {
3112                     __builtin_mma_xxsetaccz(&acc[x]);
3113                }
3114                for (int64_t kk = 0; kk < kc; kk++) {
3115                    int A0_block_idx = A0_base + kk * 32;
3116                    int B0_block_idx = B0_base + kk * 32;
3117                    int A1_block_idx = A0_block_idx + 32 * kc;
3118                    int B1_block_idx = B0_block_idx + 32 * kc;
3119                    vec_t * A0_block = & vec_A[A0_block_idx];
3120                    vec_t * B0_block = & vec_B[B0_block_idx];
3121                    vec_t * A1_block = & vec_A[A1_block_idx];
3122                    for (int it = 0; it < 4; it++) {
3123                        for (int x = 0; x < 4; x++) {
3124                            __builtin_mma_xvf16ger2pp(& acc[0], A0_block[8 * it + x], B0_block[8 * it + x]);
3125                            __builtin_mma_xvf16ger2pp(& acc[1], A0_block[8 * it + x], B0_block[8 * it + x + 4]);
3126                            __builtin_mma_xvf16ger2pp(& acc[2], A0_block[8 * it + x + 4], B0_block[8 * it + x]);
3127                            __builtin_mma_xvf16ger2pp(& acc[3], A0_block[8 * it + x + 4], B0_block[8 * it + x + 4]);
3128                            __builtin_mma_xvf16ger2pp(& acc[4], A1_block[8 * it + x], B0_block[8 * it + x]);
3129                            __builtin_mma_xvf16ger2pp(& acc[5], A1_block[8 * it + x], B0_block[8 * it+ x + 4]);
3130                            __builtin_mma_xvf16ger2pp(& acc[6], A1_block[8 * it + x + 4], B0_block[8 * it + x]);
3131                            __builtin_mma_xvf16ger2pp(& acc[7], A1_block[8 * it + x + 4], B0_block[8 * it + x + 4]);
3132                        }
3133                    }
3134                }
3135                if (l == 0) {
3136                    save_acc(& acc[0], ii + i, jj + j);
3137                    save_acc(& acc[1], ii + i, jj + j + 4);
3138                    save_acc(& acc[2], ii + i + 4, jj + j);
3139                    save_acc(& acc[3], ii + i + 4, jj + j + 4);
3140                    save_acc(& acc[4], ii + i + 8, jj + j);
3141                    save_acc(& acc[5], ii + i + 8, jj + j + 4);
3142                    save_acc(& acc[6], ii + i + 12, jj + j);
3143                    save_acc(& acc[7], ii + i + 12, jj + j + 4);
3144                } else {
3145                    add_save_acc(& acc[0], ii + i, jj + j);
3146                    add_save_acc(& acc[1], ii + i, jj + j + 4);
3147                    add_save_acc(& acc[2], ii + i + 4, jj + j);
3148                    add_save_acc(& acc[3], ii + i + 4, jj + j + 4);
3149                    add_save_acc(& acc[4], ii + i + 8, jj + j);
3150                    add_save_acc(& acc[5], ii + i + 8, jj + j + 4);
3151                    add_save_acc(& acc[6], ii + i + 12, jj + j);
3152                    add_save_acc(& acc[7], ii + i + 12, jj + j + 4);
3153                }
3154            }
3155        }
3156    }
3157
3158    void matmul_tiled(int64_t m, int64_t n, int64_t mc, int64_t nc, int64_t kc) {
3159        vec_t A_pack[mc * kc * 4];
3160        vec_t B_pack[nc * kc * 4];
3161        constexpr bool is_Ablock_q4 = std::is_same_v<TA, block_q4_0>;
3162        int64_t ytiles = m / mc;
3163        int64_t xtiles = n / nc;
3164        int64_t tiles  = xtiles * ytiles;
3165        int64_t duty = (tiles + nth - 1) / nth;
3166        int64_t start = duty * ith;
3167        int64_t end = start + duty;
3168        if (end > tiles) {
3169            end = tiles;
3170        }
3171        for (int64_t job = start; job < end; ++job) {
3172            int64_t ii = (job / xtiles) * mc;
3173            int64_t jj = (job % xtiles) * nc;
3174            for (int64_t kk = 0; kk < k; kk += kc) {
3175                int64_t k_cur = MIN(kc, k - kk);
3176                if constexpr(is_Ablock_q4) {
3177                    packNormal_q4_fp16(A + ii * lda + kk, lda, mc, k_cur, (uint8_t *)A_pack);
3178                } else {
3179                    packNormal_q8_fp16(A + ii * lda + kk, lda, mc, k_cur, (uint8_t *)A_pack);
3180                }
3181                packNormal_q8_fp16(B + jj * ldb + kk, ldb, nc, k_cur, (uint8_t *)B_pack);
3182                KERNEL_Q0(ii, jj, mc, nc, k_cur, kk, A_pack, B_pack);
3183            }
3184        }
3185    }
3186
3187    void gemm_small(int64_t m0, int64_t m, int64_t n0, int64_t n, int RM, int RN) {
3188        int64_t ytiles = (m - m0) / RM;
3189        int64_t xtiles = (n - n0) / RN;
3190        int64_t tiles = xtiles * ytiles;
3191        int64_t duty = (tiles + nth - 1) / nth;
3192        int64_t start = duty * ith;
3193        int64_t end = start + duty;
3194        vec_t vec_A[8] = {0}, vec_B[8] = {0};
3195        vector signed int vec_C[4];
3196        acc_t acc_0;
3197        bool isAblock_q4 = std::is_same_v<TA, block_q4_0>;
3198
3199        if (end > tiles)
3200            end = tiles;
3201        for (int64_t job = start; job < end; ++job) {
3202            int64_t ii = m0 + job / xtiles * RM;
3203            int64_t jj = n0 + job % xtiles * RN;
3204            std::array<int, 4> comparray{};
3205            vector float res[4] = {0};
3206            vector float fin_res[4] = {0};
3207            vector float vs[4] = {0};
3208            vector float CA[4] = {0};
3209            __builtin_prefetch((A + (ii * lda) + 0)->qs, 0, 1); // prefetch first value
3210            __builtin_prefetch((B + (jj * ldb) + 0)->qs, 0, 1); // prefetch first value
3211            for (int l = 0; l < k; l++) {
3212                __builtin_prefetch((A + (ii * lda) + (l + 1))->qs, 0, 1); // prefetch one loop ahead
3213                __builtin_prefetch((B + (jj * ldb) + (l + 1))->qs, 0, 1); // prefetch one loop ahead
3214                __builtin_mma_xxsetaccz(& acc_0);
3215                if (isAblock_q4) {
3216                    packNormalInt4<4>((A + (ii * lda) + l), lda, RM, 4, (int8_t *)vec_A, comparray);
3217                } else {
3218                    packNormal<int8_t, vector signed char>((const block_q8_0 *)(A + (ii * lda) + l), lda, RM, 8, (int8_t *)vec_A, false);
3219                }
3220                packNormal<uint8_t, vector unsigned char>((B + (jj * ldb) + l), ldb, RN, 8, (uint8_t *)vec_B, true);
3221                for (int x = 0; x < 8; x += 4) {
3222                    __builtin_mma_xvi8ger4pp(& acc_0, vec_A[x], vec_B[x]);
3223                    __builtin_mma_xvi8ger4pp(& acc_0, vec_A[x + 1], vec_B[x + 1]);
3224                    __builtin_mma_xvi8ger4pp(& acc_0, vec_A[x + 2], vec_B[x + 2]);
3225                    __builtin_mma_xvi8ger4pp(& acc_0, vec_A[x + 3], vec_B[x + 3]);
3226                }
3227                for (int I = 0; I < RM; I++) {
3228                    for (int J = 0; J < RN; J++) {
3229                        *((float*)&vs[I] + J) = (unhalf((A + ((ii + I) * lda) + l)->d) * unhalf((B + ((jj + J) * ldb) + l)->d));
3230                    }
3231                }
3232                __builtin_mma_disassemble_acc(vec_C, & acc_0);
3233                if (!isAblock_q4) {
3234                    auto aoffset = A + (ii * lda) + l;
3235                    for (int i = 0; i < RM; i++) {
3236                        comparray[i] = 0;
3237                        int ca = 0;
3238                        auto *at = aoffset->qs;
3239                        for (int j = 0; j < 32; j++)
3240                            ca += (int)*at++;
3241                        comparray[i] = ca;
3242                        aoffset += lda;
3243                    }
3244                }
3245                for (int i = 0; i < RM; i++) {
3246                    CA[i] = vec_splats((float)(((double)comparray[i]) * -128.0));
3247                    res[i] = vec_add(vec_ctf(vec_C[i], 0), CA[i]);
3248                    fin_res[i] = vec_madd(res[i], vs[i], fin_res[i]);
3249                }
3250            }
3251            save_res(ii, jj, 0, fin_res, RM, RN);
3252        }
3253    }
3254
3255    template<int RM, int RN>
3256    inline void kernel(int64_t ii, int64_t jj) {
3257        if constexpr(RM == 4 && RN == 8) {
3258            KERNEL_4x8(ii,jj);
3259        } else if constexpr(RM == 8 && RN == 4) {
3260            KERNEL_8x4(ii,jj);
3261        } else if constexpr(RM == 8 && RN == 8) {
3262            KERNEL_8x8(ii,jj);
3263        } else {
3264            assert(false && "RN/RM values not supported");
3265        }
3266    }
3267
3268    template <int RM, int RN>
3269    NOINLINE void gemm(int64_t m0, int64_t m, int64_t n0, int64_t n) {
3270        int64_t ytiles = (m - m0) / RM;
3271        int64_t xtiles = (n - n0) / RN;
3272        int64_t tiles = xtiles * ytiles;
3273        int64_t duty = (tiles + nth - 1) / nth;
3274        int64_t start = duty * ith;
3275        int64_t end = start + duty;
3276        if (end > tiles)
3277            end = tiles;
3278        for (int64_t job = start; job < end; ++job) {
3279            int64_t ii = m0 + job / xtiles * RM;
3280            int64_t jj = n0 + job % xtiles * RN;
3281            kernel<RM, RN>(ii, jj);
3282        }
3283    }
3284    const TA * const A;
3285    const block_q8_0 * const B;
3286    float * C;
3287    const int64_t k;
3288    int64_t kc;
3289    const int64_t lda;
3290    const int64_t ldb;
3291    const int64_t ldc;
3292    const int ith;
3293    const int nth;
3294};
3295
3296class tinyBLAS_PPC {
3297  public:
3298    tinyBLAS_PPC(int64_t k,
3299                const float * A, int64_t lda,
3300                const float * B, int64_t ldb,
3301                float * C, int64_t ldc,
3302                int ith, int nth)
3303        : A(A), B(B), C(C), k(k), lda(lda), ldb(ldb), ldc(ldc), ith(ith), nth(nth) {
3304    }
3305
3306    void matmul(int64_t m, int64_t n) {
3307        int64_t mc = 256;
3308        int64_t nc = 256;
3309        int64_t kc = 256;
3310    #if defined(_AIX) || defined(__BIG_ENDIAN__)
3311        mc = 128;
3312        nc = 128;
3313        kc = 128;
3314    #endif
3315        if (m % mc == 0 && n % nc == 0 && k % kc == 0) {
3316            matmul_tiled(m, n, mc, nc, kc);
3317        } else {
3318            mnpack(0, m, 0, n);
3319        }
3320    }
3321
3322  private:
3323
3324    __attribute__((always_inline))
3325    inline void save_acc(acc_t * ACC, int64_t ii, int64_t jj) {
3326        vec_t vec_C[4];
3327        __builtin_mma_disassemble_acc(vec_C, ACC);
3328        for (int I = 0; I < 4; I++) {
3329            for (int J = 0; J < 4; J++) {
3330                *((float *)(C+ii+((jj+J)*ldc)+I)) = *((float *)&vec_C[I]+J);
3331            }
3332        }
3333    }
3334
3335    __attribute__((always_inline))
3336    inline void add_save_acc(acc_t * ACC, int64_t ii, int64_t jj) {
3337        vec_t vec_C[4];
3338        __builtin_mma_disassemble_acc(vec_C, ACC);
3339        for (int I = 0; I < 4; I++) {
3340            for (int J = 0; J < 4; J++) {
3341                float * c_ptr = (float *)(C+ii+((jj+J)*ldc)+I);
3342                *c_ptr += *((float *)&vec_C[I]+J);
3343            }
3344        }
3345    }
3346
3347    inline void vector_permute_store_4(vector float * src, float * vecOffset) {
3348        vector float t1, t2, t3, t4, t5, t6, t7, t8;
3349        t1 = vec_mergeh(src[0], src[1]);
3350        t2 = vec_mergeh(src[2], src[3]);
3351        t3 = vec_mergel(src[0], src[1]);
3352        t4 = vec_mergel(src[2], src[3]);
3353
3354        t5 = vec_xxpermdi(t1, t2, 0);
3355        t6 = vec_xxpermdi(t1, t2, 3);
3356        t7 = vec_xxpermdi(t3, t4, 0);
3357        t8 = vec_xxpermdi(t3, t4, 3);
3358
3359        vec_xst(t5, 0, vecOffset);
3360        vec_xst(t6, 0, vecOffset + 4);
3361        vec_xst(t7, 0, vecOffset + 8);
3362        vec_xst(t8, 0, vecOffset + 12);
3363    }
3364
3365    inline void vector_permute_store_8(vector float * src, float * vecOffset) {
3366        vector float t1, t2, t3, t4, t5, t6, t7, t8;
3367        t1 = vec_mergeh(src[0], src[1]);
3368        t2 = vec_mergeh(src[2], src[3]);
3369        t3 = vec_mergeh(src[4], src[5]);
3370        t4 = vec_mergeh(src[6], src[7]);
3371
3372        t5 = vec_xxpermdi(t1, t2, 0);
3373        t6 = vec_xxpermdi(t3, t4, 0);
3374        t7 = vec_xxpermdi(t1, t2, 3);
3375        t8 = vec_xxpermdi(t3, t4, 3);
3376
3377        vec_xst(t5, 0, vecOffset);
3378        vec_xst(t6, 0, vecOffset + 4);
3379        vec_xst(t7, 0, vecOffset + 8);
3380        vec_xst(t8, 0, vecOffset + 12);
3381
3382        t1 = vec_mergel(src[0], src[1]);
3383        t2 = vec_mergel(src[2], src[3]);
3384        t3 = vec_mergel(src[4], src[5]);
3385        t4 = vec_mergel(src[6], src[7]);
3386
3387        t5 = vec_xxpermdi(t1, t2, 0);
3388        t6 = vec_xxpermdi(t3, t4, 0);
3389        t7 = vec_xxpermdi(t1, t2, 3);
3390        t8 = vec_xxpermdi(t3, t4, 3);
3391
3392        vec_xst(t5, 0, vecOffset + 16);
3393        vec_xst(t6, 0, vecOffset + 20);
3394        vec_xst(t7, 0, vecOffset + 24);
3395        vec_xst(t8, 0, vecOffset + 28);
3396    }
3397
3398    void packTranspose(const float * a, int64_t lda, int rows, int cols, float * vec) {
3399        int64_t i, j;
3400        float * aoffsets[8];
3401        float * aoffset = NULL, * boffset = NULL;
3402        __vector_pair arr[8];
3403        vector float c[8][2] = {0};
3404        vector float c1[8] = {0};
3405        vector float c2[8] = {0};
3406        aoffset = const_cast<float *>(a);
3407        boffset = vec;
3408        j = (rows >> 3);
3409        if (j > 0) {
3410            do {
3411                aoffsets[0] = aoffset;
3412                for (int it = 1; it < 8; it++)
3413                    aoffsets[it] = aoffsets[it-1] + lda;
3414                aoffset += 8 * lda;
3415                i = (cols >> 3);
3416                if (i > 0) {
3417                    do {
3418                        for (int it = 0; it < 8; it++) {
3419                            arr[it] = __builtin_vsx_lxvp(0, (__vector_pair*)aoffsets[it]);
3420                            __builtin_vsx_disassemble_pair(c[it], &arr[it]);
3421                            c1[it] = c[it][0];
3422                            c2[it] = c[it][1];
3423                        }
3424
3425                        vector_permute_store_8(c1, boffset);
3426                        vector_permute_store_8(c2, boffset + 32);
3427                        boffset += 64;
3428                        i--;
3429                        if (i > 0) {
3430                           for (int it = 0; it < 8; it++) {
3431                               aoffsets[it] = aoffsets[it] + 8;
3432                           }
3433                        }
3434                    } while(i > 0);
3435                }
3436                if (cols & 4) {
3437                    for (int it = 0; it < 8 ; it++)
3438                        c1[it] = vec_xl(0, aoffsets[it]);
3439                    vector_permute_store_8(c1, boffset);
3440                }
3441            j--;
3442            } while(j > 0);
3443        }
3444
3445        if (rows & 4) {
3446            aoffsets[0] = aoffset;
3447            for (int it = 1; it < 4; it++)
3448                aoffsets[it] = aoffsets[it-1] + lda;
3449            aoffset += 4 * lda;
3450            i = (cols >> 3);
3451            if (i > 0) {
3452                do {
3453                    for (int it = 0; it < 4; it++) {
3454                        arr[it] = __builtin_vsx_lxvp(0, (__vector_pair*)aoffsets[it]);
3455                        __builtin_vsx_disassemble_pair(c[it], &arr[it]);
3456                        c1[it] = c[it][0];
3457                        c2[it] = c[it][1];
3458                    }
3459                    vector_permute_store_4(c1, boffset);
3460                    vector_permute_store_4(c2, boffset + 16);
3461                    for (int it = 0; it < 4; it++)
3462                        aoffsets[it] += 8 * lda;
3463                    boffset += 32;
3464                    i--;
3465                } while(i > 0);
3466            }
3467
3468            if (cols & 4) {
3469               for (int it = 0; it < 4; it++)
3470                   c1[it] = vec_xl(0, aoffsets[it]);
3471                vector_permute_store_4(c1, boffset);
3472            }
3473        }
3474        if (rows & 3) {
3475            aoffsets[0] = aoffset;
3476            for (int it = 1; it < 3; it++)
3477                aoffsets[it] = aoffsets[it-1] + lda;
3478            if (cols & 4) {
3479                for (int it = 0; it < 3; it++)
3480                    c1[it] = vec_xl(0, aoffsets[it]);
3481                vector_permute_store_4(c1, boffset);
3482            }
3483        }
3484    }
3485
3486    void KERNEL_4x4(int64_t ii, int64_t jj) {
3487        vec_t vec_A[4], vec_B[4], vec_C[4];
3488        acc_t acc_0;
3489        __builtin_mma_xxsetaccz(&acc_0);
3490        for (int l = 0; l < k; l += 4) {
3491            packTranspose(A + (ii * lda) + l, lda, 4, 4, (float *)vec_A);
3492            packTranspose(B + (jj * ldb) + l, ldb, 4, 4, (float *)vec_B);
3493            __builtin_mma_xvf32gerpp(&acc_0, vec_A[0], vec_B[0]);
3494            __builtin_mma_xvf32gerpp(&acc_0, vec_A[1], vec_B[1]);
3495            __builtin_mma_xvf32gerpp(&acc_0, vec_A[2], vec_B[2]);
3496            __builtin_mma_xvf32gerpp(&acc_0, vec_A[3], vec_B[3]);
3497        }
3498        save_acc(&acc_0, ii, jj);
3499    }
3500
3501    void KERNEL_4x8(int64_t ii, int64_t jj) {
3502        vec_t vec_A[4], vec_B[8], vec_C[4];
3503        acc_t acc_0, acc_1;
3504        __builtin_mma_xxsetaccz(&acc_0);
3505        __builtin_mma_xxsetaccz(&acc_1);
3506        for (int64_t l = 0; l < k; l += 4) {
3507            packTranspose(A + (ii * lda) + l, lda, 4, 4, (float *)vec_A);
3508            packTranspose(B + (jj * ldb) + l, ldb, 8, 4, (float *)vec_B);
3509            __builtin_mma_xvf32gerpp(&acc_0, vec_A[0], (vec_t)vec_B[0]);
3510            __builtin_mma_xvf32gerpp(&acc_1, vec_A[0], (vec_t)vec_B[1]);
3511            __builtin_mma_xvf32gerpp(&acc_0, vec_A[1], (vec_t)vec_B[2]);
3512            __builtin_mma_xvf32gerpp(&acc_1, vec_A[1], (vec_t)vec_B[3]);
3513            __builtin_mma_xvf32gerpp(&acc_0, vec_A[2], (vec_t)vec_B[4]);
3514            __builtin_mma_xvf32gerpp(&acc_1, vec_A[2], (vec_t)vec_B[5]);
3515            __builtin_mma_xvf32gerpp(&acc_0, vec_A[3], (vec_t)vec_B[6]);
3516            __builtin_mma_xvf32gerpp(&acc_1, vec_A[3], (vec_t)vec_B[7]);
3517        }
3518        save_acc(&acc_0, ii, jj);
3519        save_acc(&acc_1, ii, jj + 4);
3520    }
3521
3522    void KERNEL_8x4(int64_t ii, int64_t jj) {
3523        vec_t vec_A[8], vec_B[4], vec_C[4];
3524        acc_t acc_0, acc_1;
3525        __builtin_mma_xxsetaccz(&acc_0);
3526        __builtin_mma_xxsetaccz(&acc_1);
3527        for (int64_t l = 0; l < k; l += 4) {
3528            packTranspose(A + (ii * lda) + l, lda, 8, 4, (float *)vec_A);
3529            packTranspose(B + (jj * ldb) + l, ldb, 4, 4, (float *)vec_B);
3530            __builtin_mma_xvf32gerpp(&acc_0, (vec_t)vec_A[0], vec_B[0]);
3531            __builtin_mma_xvf32gerpp(&acc_1, (vec_t)vec_A[1], vec_B[0]);
3532            __builtin_mma_xvf32gerpp(&acc_0, (vec_t)vec_A[2], vec_B[1]);
3533            __builtin_mma_xvf32gerpp(&acc_1, (vec_t)vec_A[3], vec_B[1]);
3534            __builtin_mma_xvf32gerpp(&acc_0, (vec_t)vec_A[4], vec_B[2]);
3535            __builtin_mma_xvf32gerpp(&acc_1, (vec_t)vec_A[5], vec_B[2]);
3536            __builtin_mma_xvf32gerpp(&acc_0, (vec_t)vec_A[6], vec_B[3]);
3537            __builtin_mma_xvf32gerpp(&acc_1, (vec_t)vec_A[7], vec_B[3]);
3538        }
3539        save_acc(&acc_0, ii, jj);
3540        save_acc(&acc_1, ii + 4, jj);
3541    }
3542
3543    void KERNEL_8x8(int64_t ii, int64_t jj) {
3544        vec_t vec_A[16], vec_B[16], vec_C[4];
3545        acc_t acc_0, acc_1, acc_2, acc_3;
3546        __builtin_mma_xxsetaccz(&acc_0);
3547        __builtin_mma_xxsetaccz(&acc_1);
3548        __builtin_mma_xxsetaccz(&acc_2);
3549        __builtin_mma_xxsetaccz(&acc_3);
3550        for (int l = 0; l < k; l+=8) {
3551            packTranspose(A + (ii * lda) + l, lda, 8, 8, (float *)vec_A);
3552            packTranspose(B + (jj * ldb) + l, ldb, 8, 8, (float *)vec_B);
3553            for(int x = 0; x < 16; x+=2) {
3554                __builtin_mma_xvf32gerpp(&acc_0, (vec_t)vec_A[x], vec_B[x]);
3555                __builtin_mma_xvf32gerpp(&acc_1, (vec_t)vec_A[x], vec_B[x + 1]);
3556                __builtin_mma_xvf32gerpp(&acc_2, (vec_t)vec_A[x + 1], vec_B[x]);
3557                __builtin_mma_xvf32gerpp(&acc_3, (vec_t)vec_A[x + 1], vec_B[x + 1]);
3558            }
3559        }
3560        save_acc(&acc_0, ii, jj);
3561        save_acc(&acc_1, ii, jj + 4);
3562        save_acc(&acc_2, ii + 4, jj);
3563        save_acc(&acc_3, ii + 4, jj + 4);
3564    }
3565
3566    inline void MMA_16x8(vec_t * vec_A0, vec_t * vec_A1, vec_t * vec_B, acc_t * acc) {
3567        for (int x = 0; x < 16; x += 2) {
3568            __builtin_mma_xvf32gerpp(&acc[0], vec_A0[x + 0], vec_B[x]);
3569            __builtin_mma_xvf32gerpp(&acc[1], vec_A0[x + 0], vec_B[x + 1]);
3570            __builtin_mma_xvf32gerpp(&acc[2], vec_A0[x + 1], vec_B[x]);
3571            __builtin_mma_xvf32gerpp(&acc[3], vec_A0[x + 1], vec_B[x + 1]);
3572            __builtin_mma_xvf32gerpp(&acc[4], vec_A1[x + 0], vec_B[x]);
3573            __builtin_mma_xvf32gerpp(&acc[5], vec_A1[x + 0], vec_B[x + 1]);
3574            __builtin_mma_xvf32gerpp(&acc[6], vec_A1[x + 1], vec_B[x]);
3575            __builtin_mma_xvf32gerpp(&acc[7], vec_A1[x + 1], vec_B[x + 1]);
3576        }
3577    }
3578
3579    void KERNEL(int64_t ii, int64_t jj, int64_t mc, int64_t nc, int64_t kc, vec_t * vec_A, vec_t * vec_B, int64_t kk) {
3580        for (int64_t i = 0; i < mc; i += 16) {
3581            int A_base_addr = (mc / 8) * (i / 8) * 16;
3582            for (int64_t j = 0; j < nc; j += 8) {
3583                 int B_base_addr = (nc / 8) * (j / 8) * 16;
3584                 acc_t acc[8];
3585                 vec_t A0_block[16]; vec_t A1_block[16];
3586                 for (int x = 0; x < 8; x++)
3587                     __builtin_mma_xxsetaccz(&acc[x]);
3588                 for (int64_t l = 0; l < kc; l += 8) {
3589                     int A0_block_idx = A_base_addr + (l / 8) * 16;
3590                     int A1_block_idx = A0_block_idx + (mc / 8) * 16;
3591                     int B_block_idx = B_base_addr + (l / 8) * 16;
3592                     vec_t* A0_block = &vec_A[A0_block_idx];
3593                     vec_t* A1_block = &vec_A[A1_block_idx];
3594                     vec_t* B_block = &vec_B[B_block_idx];
3595                     MMA_16x8(A0_block, A1_block, B_block, acc);
3596                 }
3597                 if (kk == 0) {
3598                     save_acc(&acc[0], ii + i, jj + j);
3599                     save_acc(&acc[1], ii + i, jj + j + 4);
3600                     save_acc(&acc[2], ii + i + 4, jj + j);
3601                     save_acc(&acc[3], ii + i + 4, jj + j + 4);
3602                     save_acc(&acc[4], ii + i + 8, jj + j);
3603                     save_acc(&acc[5], ii + i + 8, jj + j + 4);
3604                     save_acc(&acc[6], ii + i + 12, jj + j);
3605                     save_acc(&acc[7], ii + i + 12, jj + j + 4);
3606                 } else {
3607                     add_save_acc(&acc[0], ii + i, jj + j);
3608                     add_save_acc(&acc[1], ii + i, jj + j + 4);
3609                     add_save_acc(&acc[2], ii + i + 4, jj + j);
3610                     add_save_acc(&acc[3], ii + i + 4, jj + j + 4);
3611                     add_save_acc(&acc[4], ii + i + 8, jj + j);
3612                     add_save_acc(&acc[5], ii + i + 8, jj + j + 4);
3613                     add_save_acc(&acc[6], ii + i + 12, jj + j);
3614                     add_save_acc(&acc[7], ii + i + 12, jj + j + 4);
3615                 }
3616            }
3617        }
3618    }
3619
3620    void matmul_tiled(int64_t m , int64_t n, int64_t mc, int64_t nc, int64_t kc) {
3621        int64_t ytiles = m / mc;
3622        int64_t xtiles = n / nc;
3623        int64_t tiles = xtiles * ytiles;
3624        int64_t duty = (tiles + nth - 1) / nth;
3625        int64_t start = duty * ith;
3626        int64_t end = start + duty;
3627        if (end > tiles) {
3628            end = tiles;
3629        }
3630        for (int64_t job = start; job < end; ++job) {
3631            int64_t ii = (job / xtiles) * mc;
3632            int64_t jj = (job % xtiles) * nc;
3633            for (int64_t kk = 0; kk < k; kk += kc) {
3634                 vec_t A_pack[kc * mc / 4];
3635                 vec_t B_pack[kc * nc / 4];
3636                 packTranspose(A + (ii * lda) + kk, lda, kc, mc, (float *)A_pack);
3637                 packTranspose(B + (jj * ldb) + kk, ldb, kc, nc, (float *)B_pack);
3638                 KERNEL(ii, jj, mc, nc, kc, A_pack, B_pack, kk);
3639            }
3640        }
3641    }
3642
3643    void mnpack(int64_t m0, int64_t m, int64_t n0, int64_t n) {
3644        int m_rem = MIN(m - m0, 8);
3645        int n_rem = MIN(n - n0, 8);
3646        int mc = 0, nc = 0;
3647        if (m_rem >= 8 && n_rem >= 8) {
3648            mc = 8;
3649            nc = 8;
3650            gemm<8, 8>(m0, m, n0, n);
3651        } else if (m_rem >= 4 && n_rem >= 8) {
3652            mc = 4;
3653            nc = 8;
3654            gemm<4, 8>(m0, m, n0, n);
3655        } else if (m_rem >= 8 && n_rem >= 4) {
3656            mc = 8;
3657            nc = 4;
3658            gemm<8, 4>(m0, m, n0, n);
3659        } else if (m_rem >= 4 && n_rem >= 4) {
3660            mc = 4;
3661            nc = 4;
3662            gemm<4, 4>(m0, m, n0, n);
3663        } else {
3664            mc = (m_rem >= 4) ? 4 : m_rem;
3665            nc = (n_rem >= 4) ? 4 : n_rem;
3666            if (mc == 0 || nc == 0)
3667                return;
3668            gemm_small(m0, m, n0, n, mc, nc);
3669        }
3670        int64_t mp = m0 + ((m - m0) / mc) * mc;
3671        int64_t np = n0 + ((n - n0) / nc) * nc;
3672        mnpack(mp, m, n0, np);
3673        mnpack(m0, m, np, n);
3674    }
3675
3676    void gemm_small(int64_t m0, int64_t m, int64_t n0, int64_t n, int RM, int RN) {
3677        int64_t ytiles = (m - m0) / RM;
3678        int64_t xtiles = (n - n0) / RN;
3679        int64_t tiles = xtiles * ytiles;
3680        int64_t duty = (tiles + nth - 1) / nth;
3681        int64_t start = duty * ith;
3682        int64_t end = start + duty;
3683        if (end > tiles)
3684            end = tiles;
3685        for (int64_t job = start; job < end; ++job) {
3686            int64_t ii = m0 + job / xtiles * RM;
3687            int64_t jj = n0 + job % xtiles * RN;
3688            vec_t vec_C[4];
3689            acc_t acc_0;
3690            __builtin_mma_xxsetaccz(&acc_0);
3691            vec_t vec_A[4] = {0}, vec_B[4] = {0};
3692            for (int l = 0; l < k; l += 4) {
3693                /* 'GEMV Forwarding' concept is used in first two conditional loops.
3694                 * when one of the matrix has a single row/column, the elements are
3695                 * broadcasted, instead of using packing routine to prepack the
3696                 * matrix elements.
3697                 */
3698                if (RM == 1) {
3699                    float * a = const_cast<float *>(A + (ii) * lda + l);
3700                    packTranspose(B + (jj * ldb) + l, ldb, RN, 4, (float *)vec_B);
3701                    vec_A[0] = (vec_t)vec_xl(0,a);
3702                    vec_A[1] = (vec_t)vec_splats(*((float *)&vec_A+1));
3703                    vec_A[2] = (vec_t)vec_splats(*((float *)&vec_A+2));
3704                    vec_A[3] = (vec_t)vec_splats(*((float *)&vec_A+3));
3705                } else if (RN == 1) {
3706                    packTranspose(A + (ii * lda) + l, lda, RM, 4, (float *)vec_A);
3707                    float * b = const_cast<float *>(B + (jj) * ldb + l);
3708                    vec_B[0] = (vec_t)vec_xl(0,b);
3709                    vec_B[1] = (vec_t)vec_splats(*((float *)&vec_B+1));
3710                    vec_B[2] = (vec_t)vec_splats(*((float *)&vec_B+2));
3711                    vec_B[3] = (vec_t)vec_splats(*((float *)&vec_B+3));
3712                } else {
3713                    packTranspose(A + (ii * lda) + l, lda, RM, 4, (float *)vec_A);
3714                    packTranspose(B + (jj * ldb) + l, ldb, RN, 4, (float *)vec_B);
3715                }
3716                __builtin_mma_xvf32gerpp(&acc_0, vec_A[0], vec_B[0]);
3717                __builtin_mma_xvf32gerpp(&acc_0, vec_A[1], vec_B[1]);
3718                __builtin_mma_xvf32gerpp(&acc_0, vec_A[2], vec_B[2]);
3719                __builtin_mma_xvf32gerpp(&acc_0, vec_A[3], vec_B[3]);
3720            }
3721            __builtin_mma_disassemble_acc(vec_C, &acc_0);
3722            for (int I = 0; I < RM; I++) {
3723                for (int J = 0; J < RN; J++) {
3724                    *((float *)(C+ii+((jj+J)*ldc)+I)) = *((float *)&vec_C[I]+J);
3725                }
3726            }
3727       }
3728    }
3729
3730    template<int RM, int RN>
3731    inline void kernel(int64_t ii, int64_t jj) {
3732        if constexpr(RM == 4 && RN == 4) {
3733            KERNEL_4x4(ii, jj);
3734        } else if constexpr(RM == 4 && RN == 8) {
3735            KERNEL_4x8(ii, jj);
3736        } else if constexpr(RM == 8 && RN == 4) {
3737            KERNEL_8x4(ii, jj);
3738        } else if constexpr(RM == 8 && RN == 8) {
3739            KERNEL_8x8(ii, jj);
3740        } else {
3741            static_assert(false, "RN/RM values not supported");
3742        }
3743    }
3744
3745    template <int RM, int RN>
3746    NOINLINE void gemm(int64_t m0, int64_t m, int64_t n0, int64_t n) {
3747        int64_t ytiles = (m - m0) / RM;
3748        int64_t xtiles = (n - n0) / RN;
3749        int64_t tiles = xtiles * ytiles;
3750        int64_t duty = (tiles + nth - 1) / nth;
3751        int64_t start = duty * ith;
3752        int64_t end = start + duty;
3753        if (end > tiles)
3754            end = tiles;
3755        for (int64_t job = start; job < end; ++job) {
3756            int64_t ii = m0 + job / xtiles * RM;
3757            int64_t jj = n0 + job % xtiles * RN;
3758            kernel<RM, RN>(ii, jj);
3759        }
3760    }
3761
3762    const float * const A;
3763    const float * const B;
3764    float * C;
3765    const int64_t k;
3766    const int64_t lda;
3767    const int64_t ldb;
3768    const int64_t ldc;
3769    const int ith;
3770    const int nth;
3771};
3772#endif
3773} // namespace
3774
3775/**
3776 * Performs optimized matrix multiplication on CPU.
3777 *
3778 * This subroutine may compute C = Aᵀ * B with column major ordering.
3779 * Despite its name, this isn't a generalized implementation. Work is
3780 * only performed when a handwritten kernel is written and available.
3781 * Otherwise the caller should fall back to a general matmul routine.
3782 *
3783 * For example, for single-threaded single-precision GEMM you can say
3784 *
3785 *     llamafile_sgemm(m, n, k, A, lda, B, ldb, C, ldc,
3786 *                     0, 1,
3787 *                     GGML_TYPE_F32, GGML_TYPE_F32, GGML_TYPE_F32);
3788 *
3789 * @param m is rows in `A` and `C`
3790 * @param n is cols in `B` and `C`
3791 * @param k is cols in `A` and rows in `B`
3792 * @param A is first input matrix (always transposed)
3793 * @param lda is row stride of `A`
3794 * @param B is second input matrix (never transposed)
3795 * @param ldb is row stride of `B`
3796 * @param C is input/output array of output matrices
3797 * @param ldc is row stride of `C`
3798 * @param ith is thread id (must be less than `nth`)
3799 * @param nth is number of threads (must be greater than zero)
3800 * @param Atype is GGML data type of `A`
3801 * @param Btype is GGML data type of `B`
3802 * @param Ctype is GGML data type of `C`
3803 * @return true if this function was able to service the matmul request
3804 */
3805bool llamafile_sgemm(const struct ggml_compute_params * params, int64_t m, int64_t n, int64_t k,
3806                     const void *A, int64_t lda, const void *B, int64_t ldb, void *C,
3807                     int64_t ldc, int Atype, int Btype, int Ctype) {
3808
3809    assert(m >= 0);
3810    assert(n >= 0);
3811    assert(k >= 0);
3812    assert(lda >= k);
3813    assert(ldb >= k);
3814    assert(ldc >= m);
3815    assert(params->nth > 0);
3816    assert(params->ith < params->nth);
3817
3818    // only enable sgemm for prompt processing
3819#if !defined(__MMA__)
3820    if (n < 2)
3821        return false;
3822#endif
3823
3824    if (Ctype != GGML_TYPE_F32)
3825        return false;
3826
3827    switch (Atype) {
3828
3829    case GGML_TYPE_F32: {
3830        if (Btype != GGML_TYPE_F32)
3831            return false;
3832#if defined(__AVX512F__)
3833        tinyBLAS<16, __m512, __m512, float, float, float> tb{ params,
3834            k, (const float *)A, lda,
3835            (const float *)B, ldb,
3836            (float *)C, ldc};
3837        return tb.matmul(m, n);
3838#elif defined(__AVX__) || defined(__AVX2__)
3839        tinyBLAS<8, __m256, __m256, float, float, float> tb{ params,
3840            k, (const float *)A, lda,
3841            (const float *)B, ldb,
3842            (float *)C, ldc};
3843        return tb.matmul(m, n);
3844#elif defined(__ARM_NEON)
3845        if (n < 4)
3846            return false;
3847        tinyBLAS<4, float32x4_t, float32x4_t, float, float, float> tb{ params,
3848            k, (const float *)A, lda,
3849            (const float *)B, ldb,
3850            (float *)C, ldc};
3851        return tb.matmul(m, n);
3852#elif defined(__VXE__) || defined(__VXE2__)
3853        if (n < 4)
3854            return false;
3855        tinyBLAS<4, float32x4_t, float32x4_t, float, float, float> tb{ params,
3856            k, (const float *)A, lda,
3857            (const float *)B, ldb,
3858            (float *)C, ldc};
3859        return tb.matmul(m, n);
3860#elif defined(__MMA__)
3861        if (k % 8)
3862            return false;
3863        tinyBLAS_PPC tb{
3864            k, (const float *)A, lda,
3865            (const float *)B, ldb,
3866            (float *)C, ldc,
3867            params->ith, params->nth};
3868        tb.matmul(m, n);
3869        return true;
3870#elif defined(__riscv_v_intrinsic)
3871    #if LMUL == 1
3872        tinyBLAS_RVV<vfloat32m1_t, vfloat32m1_t, float, float, float> tb{ params,
3873            k, (const float *)A, lda,
3874            (const float *)B, ldb,
3875            (float *)C, ldc};
3876    #elif LMUL == 2
3877        tinyBLAS_RVV<vfloat32m2_t, vfloat32m2_t, float, float, float> tb{ params,
3878            k, (const float *)A, lda,
3879            (const float *)B, ldb,
3880            (float *)C, ldc};
3881    #else // LMUL = 4
3882        tinyBLAS_RVV<vfloat32m4_t, vfloat32m4_t, float, float, float> tb{ params,
3883            k, (const float *)A, lda,
3884            (const float *)B, ldb,
3885            (float *)C, ldc};
3886    #endif
3887        return tb.matmul(m, n);
3888#else
3889        return false;
3890#endif
3891    }
3892
3893    case GGML_TYPE_BF16: {
3894#if defined(__AVX512BF16__)
3895        if (Btype == GGML_TYPE_BF16) {
3896            tinyBLAS<32, __m512, __m512bh, ggml_bf16_t, ggml_bf16_t, float> tb{ params, k,
3897                (const ggml_bf16_t *)A, lda,
3898                (const ggml_bf16_t *)B, ldb,
3899                (float *)C, ldc};
3900            return tb.matmul(m, n);
3901        }
3902#elif defined(__AVX512F__)
3903        if (Btype == GGML_TYPE_BF16) {
3904            tinyBLAS<16, __m512, __m512, ggml_bf16_t, ggml_bf16_t, float> tb{ params, k,
3905                (const ggml_bf16_t *)A, lda,
3906                (const ggml_bf16_t *)B, ldb,
3907                (float *)C, ldc};
3908            return tb.matmul(m, n);
3909        }
3910#elif defined(__AVX2__)
3911        if (Btype == GGML_TYPE_BF16) {
3912            tinyBLAS<8, __m256, __m256, ggml_bf16_t, ggml_bf16_t, float> tb{ params, k,
3913                (const ggml_bf16_t *)A, lda,
3914                (const ggml_bf16_t *)B, ldb,
3915                (float *)C, ldc};
3916            return tb.matmul(m, n);
3917        }
3918#elif defined(__MMA__)
3919        if (k % 8) {
3920            return false;
3921        }
3922
3923        if (Btype == GGML_TYPE_BF16) {
3924            tinyBLAS_HP16_PPC<ggml_bf16_t, ggml_bf16_t, float> tb{ k,
3925                (const ggml_bf16_t *)A, lda,
3926                (const ggml_bf16_t *)B, ldb,
3927                (float *)C, ldc,
3928                params->ith, params->nth };
3929
3930            tb.matmul(m, n);
3931            return true;
3932        }
3933#elif defined(__riscv_zvfbfwma)
3934        if (Btype == GGML_TYPE_BF16) {
3935            #if LMUL == 1
3936                tinyBLAS_RVV<vfloat32m1_t, vbfloat16mf2_t, ggml_bf16_t, ggml_bf16_t, float> tb{ params,
3937                    k, (const ggml_bf16_t *)A, lda,
3938                    (const ggml_bf16_t *)B, ldb,
3939                    (float *)C, ldc};
3940            #elif LMUL == 2
3941                tinyBLAS_RVV<vfloat32m2_t, vbfloat16m1_t, ggml_bf16_t, ggml_bf16_t, float> tb{ params,
3942                    k, (const ggml_bf16_t *)A, lda,
3943                    (const ggml_bf16_t *)B, ldb,
3944                    (float *)C, ldc};
3945            #else // LMUL = 4
3946                tinyBLAS_RVV<vfloat32m4_t, vbfloat16m2_t, ggml_bf16_t, ggml_bf16_t, float> tb{ params,
3947                    k, (const ggml_bf16_t *)A, lda,
3948                    (const ggml_bf16_t *)B, ldb,
3949                    (float *)C, ldc};
3950            #endif
3951                return tb.matmul(m, n);
3952        }
3953#endif
3954        return false;
3955    }
3956
3957    case GGML_TYPE_F16: {
3958#if defined(__AVX512F__)
3959        if (Btype == GGML_TYPE_F16) {
3960            tinyBLAS<16, __m512, __m512, ggml_fp16_t, ggml_fp16_t, float> tb{ params, k,
3961                (const ggml_fp16_t *)A, lda,
3962                (const ggml_fp16_t *)B, ldb,
3963                (float *)C, ldc};
3964            return tb.matmul(m, n);
3965        }
3966#elif (defined(__AVX__) || defined(__AVX2__)) && defined(__F16C__)
3967        if (Btype == GGML_TYPE_F16) {
3968            tinyBLAS<8, __m256, __m256, ggml_fp16_t, ggml_fp16_t, float> tb{ params, k,
3969                (const ggml_fp16_t *)A, lda,
3970                (const ggml_fp16_t *)B, ldb,
3971                (float *)C, ldc};
3972            return tb.matmul(m, n);
3973        }
3974#elif defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC) && !defined(_MSC_VER)
3975        if (n < 8)
3976            return false;
3977        if (Btype == GGML_TYPE_F16) {
3978            tinyBLAS<8, float16x8_t, float16x8_t, ggml_fp16_t, ggml_fp16_t, float> tb{ params,
3979                k, (const ggml_fp16_t *)A, lda,
3980                (const ggml_fp16_t *)B, ldb,
3981                (float *)C, ldc};
3982            return tb.matmul(m, n);
3983        }
3984#elif defined(__ARM_NEON) && !defined(_MSC_VER)
3985        if (Btype == GGML_TYPE_F32) {
3986            tinyBLAS<4, float32x4_t, float32x4_t, ggml_fp16_t, float, float> tb{ params,
3987                k, (const ggml_fp16_t *)A, lda,
3988                (const float *)B, ldb,
3989                (float *)C, ldc};
3990            return tb.matmul(m, n);
3991        }
3992#elif defined(__VXE__) || defined(__VXE2__)
3993        if (n < 4)
3994            return false;
3995        if (Btype == GGML_TYPE_F16) {
3996            tinyBLAS<4, float32x4_t, float32x4_t, ggml_fp16_t, ggml_fp16_t, float> tb{ params,
3997                k, (const ggml_fp16_t *)A, lda,
3998                (const ggml_fp16_t *)B, ldb,
3999                (float *)C, ldc};
4000            return tb.matmul(m, n);
4001        }
4002#elif defined(__riscv_zvfh)
4003        if (Btype == GGML_TYPE_F16) {
4004        #if LMUL == 1
4005            tinyBLAS_RVV<vfloat32m1_t, vfloat16mf2_t, ggml_fp16_t, ggml_fp16_t, float> tb{ params,
4006                k, (const ggml_fp16_t *)A, lda,
4007                (const ggml_fp16_t *)B, ldb,
4008                (float *)C, ldc};
4009        #elif LMUL == 2
4010            tinyBLAS_RVV<vfloat32m2_t, vfloat16m1_t, ggml_fp16_t, ggml_fp16_t, float> tb{ params,
4011                k, (const ggml_fp16_t *)A, lda,
4012                (const ggml_fp16_t *)B, ldb,
4013                (float *)C, ldc};
4014        #else // LMUL = 4
4015            tinyBLAS_RVV<vfloat32m4_t, vfloat16m2_t, ggml_fp16_t, ggml_fp16_t, float> tb{ params,
4016                k, (const ggml_fp16_t *)A, lda,
4017                (const ggml_fp16_t *)B, ldb,
4018                (float *)C, ldc};
4019        #endif
4020            return tb.matmul(m, n);
4021        }
4022#elif defined(__MMA__)
4023        if (k % 8) {
4024            return false;
4025        }
4026
4027        if (Btype == GGML_TYPE_F16) {
4028            tinyBLAS_HP16_PPC<ggml_fp16_t, ggml_fp16_t, float> tb{ k,
4029                (const ggml_fp16_t *)A, lda,
4030                (const ggml_fp16_t *)B, ldb,
4031                (float *)C, ldc,
4032                params->ith, params->nth };
4033
4034            tb.matmul(m, n);
4035            return true;
4036        }
4037#endif
4038        return false;
4039    }
4040
4041    case GGML_TYPE_Q8_0: {
4042        if (Btype != GGML_TYPE_Q8_0)
4043           return false;
4044#if defined(__AVX2__) || defined(__AVX512F__) || defined(__AVX__)
4045        tinyBLAS_Q0_AVX<block_q8_0, block_q8_0, float> tb{
4046            k, (const block_q8_0 *)A, lda,
4047            (const block_q8_0 *)B, ldb,
4048            (float *)C, ldc,
4049            params->ith, params->nth};
4050        tb.matmul(m, n);
4051        return true;
4052#elif defined(__ARM_FEATURE_DOTPROD)
4053        tinyBLAS_Q0_ARM<block_q8_0> tb{
4054            k, (const block_q8_0 *)A, lda,
4055            (const block_q8_0 *)B, ldb,
4056            (float *)C, ldc,
4057            params->ith, params->nth};
4058        tb.matmul(m, n);
4059        return true;
4060#elif defined(__MMA__)
4061    //TO-DO: Remove this condition once gemv forwarding is enabled.
4062        if (n < 8 && n != 4)
4063           return false;
4064        if (m < 8 && m != 4)
4065           return false;
4066        tinyBLAS_Q0_PPC<block_q8_0> tb{
4067            k, (const block_q8_0 *)A, lda,
4068            (const block_q8_0 *)B, ldb,
4069            (float *)C, ldc,
4070            params->ith, params->nth};
4071        tb.matmul(m, n);
4072        return true;
4073#else
4074        return false;
4075#endif
4076    }
4077
4078    case GGML_TYPE_Q4_0: {
4079        if (Btype != GGML_TYPE_Q8_0)
4080            return false;
4081#if defined(__AVX2__) || defined(__AVX512F__) || defined(__AVX__)
4082        tinyBLAS_Q0_AVX<block_q4_0, block_q8_0, float> tb{
4083            k, (const block_q4_0 *)A, lda,
4084            (const block_q8_0 *)B, ldb,
4085            (float *)C, ldc,
4086            params->ith, params->nth};
4087        tb.matmul(m, n);
4088        return true;
4089#elif defined(__ARM_FEATURE_DOTPROD)
4090        tinyBLAS_Q0_ARM<block_q4_0> tb{
4091            k, (const block_q4_0 *)A, lda,
4092            (const block_q8_0 *)B, ldb,
4093            (float *)C, ldc,
4094            params->ith, params->nth};
4095        tb.matmul(m, n);
4096        return true;
4097#elif defined(__MMA__)
4098    //TO-DO: Remove this condition once gemv forwarding is enabled.
4099        if (n < 8 && n != 4)
4100           return false;
4101        if (m < 8 && m != 4)
4102           return false;
4103        tinyBLAS_Q0_PPC<block_q4_0> tb{
4104            k, (const block_q4_0 *)A, lda,
4105            (const block_q8_0 *)B, ldb,
4106            (float *)C, ldc,
4107            params->ith, params->nth};
4108        tb.matmul(m, n);
4109        return true;
4110#else
4111        return false;
4112#endif
4113    }
4114
4115    case GGML_TYPE_Q5_0: {
4116        if (Btype != GGML_TYPE_Q8_0)
4117            return false;
4118#if defined(__AVX2__) || defined(__AVX512F__) || defined(__AVX__)
4119        tinyBLAS_Q0_AVX<block_q5_0, block_q8_0, float> tb{
4120            k, (const block_q5_0 *)A, lda,
4121            (const block_q8_0 *)B, ldb,
4122            (float *)C, ldc,
4123            params->ith, params->nth};
4124        tb.matmul(m, n);
4125        return true;
4126#else
4127        return false;
4128#endif
4129    }
4130
4131    case GGML_TYPE_IQ4_NL: {
4132        if (Btype != GGML_TYPE_Q8_0)
4133            return false;
4134#if defined(__AVX2__) || defined(__AVX512F__) || defined(__AVX__)
4135        tinyBLAS_Q0_AVX<block_iq4_nl, block_q8_0, float> tb{
4136            k, (const block_iq4_nl *)A, lda,
4137            (const block_q8_0 *)B, ldb,
4138            (float *)C, ldc,
4139            params->ith, params->nth};
4140        tb.matmul(m, n);
4141        return true;
4142#else
4143        return false;
4144#endif
4145    }
4146
4147    default:
4148        return false;
4149    }
4150
4151    (void)params;
4152    (void)m;
4153    (void)n;
4154    (void)k;
4155    (void)A;
4156    (void)lda;
4157    (void)B;
4158    (void)ldb;
4159    (void)C;
4160    (void)ldc;
4161    (void)Atype;
4162    (void)Btype;
4163    (void)Ctype;
4164}