diff --git a/c/backend_cuda.cu b/c/backend_cuda.cu index 84c04de..81499ea 100644 --- a/c/backend_cuda.cu +++ b/c/backend_cuda.cu @@ -1122,7 +1122,7 @@ extern "C" int coli_cuda_attention_project_ragged(ColiCudaTensor *w,ColiCudaTens attention_absorb_ragged_kernel<<stream>>>(dc->ac,dc->aq,ddl,ddr, dn,w->weights,w->scales,w->fmt,S,H,Q,R,V,K,T,scale); quant_matmul<<O,S),256,0,dc->stream>>>(dc->y,dc->ac,proj->weights, - proj->scales,proj->fmt,S,proj->I,proj->O,row_bytes(proj->fmt,proj->I)); + proj->scales,proj->fmt,S,proj->I,proj->O,row_bytes(proj->fmt,proj->I),proj->gs,proj->ng); return cuda_ok(cudaGetLastError(),"ragged attention launch")&& cuda_ok(cudaMemcpyAsync(out,dc->y,ob,cudaMemcpyDeviceToHost,dc->stream),"ragged output download")&& cuda_ok(cudaStreamSynchronize(dc->stream),"ragged attention synchronize"); diff --git a/c/colibri.c b/c/colibri.c index 334c2df..7fc9fff 100644 --- a/c/colibri.c +++ b/c/colibri.c @@ -256,9 +256,17 @@ static inline int spec_pinned(void){ return g_spec_pin && g_spec_live; } static void matmul_qt_ex(float *y, const float *x, QT *w, int S, int allow_idot); static void matmul_qt(float *y, const float *x, QT *w, int S){ matmul_qt_ex(y,x,w,S,1); } +/* fmt=4 fused gate+up (defined later, after the quant kernels) */ +static void matmul_i4_grouped_pair(float *yg, float *yu, const float *x, + const uint8_t *qg, const float *sg, + const uint8_t *qu, const float *su, + int S, int I, int O, int gs); + static void expert_gate_up(float *g,float *u,const float *x,QT *wg,QT *wu,int S){ if(!g_no_fused_pair&&!spec_pinned()&&S==1&&wg->fmt==2&&wu->fmt==2&&wg->I==wu->I&&wg->O==wu->O) matmul_i4_pair(g,u,x,wg->q4,wg->s,wu->q4,wu->s,wg->I,wg->O); + else if(!g_no_fused_pair&&S==1&&wg->fmt==4&&wu->fmt==4&&wg->I==wu->I&&wg->O==wu->O&&wg->gs==wu->gs) + matmul_i4_grouped_pair(g,u,x,wg->q4,wg->s,wu->q4,wu->s,S,wg->I,wg->O,wg->gs); else { matmul_qt(g,x,wg,S); matmul_qt(u,x,wu,S); } } @@ -448,116 +456,9 @@ static float *falloc(int64_t n){ float *p=malloc((size_t)n*sizeof(float)); if(!p){fprintf(stderr,"OOM\n");exit(1);} return p; } -/* y[S,O] = x[S,I] @ W^T, W[O,I] f32 */ -static void matmul(float *y, const float *x, const float *W, int S, int I, int O){ - #pragma omp parallel for schedule(static) - for (int o=0;o>1))); /* 8 byte=16 nibble */ - __m128i lo=_mm_and_si128(by,m4), hi=_mm_and_si128(_mm_srli_epi16(by,4),m4); - __m128i nib=_mm_unpacklo_epi8(lo,hi); /* nibble in ordine */ - __m256 w0=_mm256_cvtepi32_ps(_mm256_sub_epi32(_mm256_cvtepu8_epi32(nib),b8)); - __m256 w1=_mm256_cvtepi32_ps(_mm256_sub_epi32(_mm256_cvtepu8_epi32(_mm_srli_si128(nib,8)),b8)); - acc=_mm256_fmadd_ps(_mm256_loadu_ps(xs+i), w0, acc); - acc=_mm256_fmadd_ps(_mm256_loadu_ps(xs+i+8), w1, acc); } - a=hsum256(acc); -#elif defined(__ARM_NEON) - const uint8x8_t m4=vdup_n_u8(0x0F); const int8x8_t b8=vdup_n_s8(8); - float32x4_t ac0=vdupq_n_f32(0), ac1=vdupq_n_f32(0); - for(;i+16<=I;i+=16){ uint8x8_t by=vld1_u8(w+(i>>1)); /* 8 byte=16 nibble */ - uint8x8x2_t z=vzip_u8(vand_u8(by,m4), vshr_n_u8(by,4)); /* nibble in ordine */ - int16x8_t w0=vmovl_s8(vsub_s8(vreinterpret_s8_u8(z.val[0]),b8)); - int16x8_t w1=vmovl_s8(vsub_s8(vreinterpret_s8_u8(z.val[1]),b8)); - ac0=vfmaq_f32(ac0, vld1q_f32(xs+i), vcvtq_f32_s32(vmovl_s16(vget_low_s16(w0)))); - ac1=vfmaq_f32(ac1, vld1q_f32(xs+i+4), vcvtq_f32_s32(vmovl_s16(vget_high_s16(w0)))); - ac0=vfmaq_f32(ac0, vld1q_f32(xs+i+8), vcvtq_f32_s32(vmovl_s16(vget_low_s16(w1)))); - ac1=vfmaq_f32(ac1, vld1q_f32(xs+i+12), vcvtq_f32_s32(vmovl_s16(vget_high_s16(w1)))); } - a=vaddvq_f32(vaddq_f32(ac0,ac1)); -#endif -#if defined(__AVX512F__) && defined(__AVX512BW__) - } -#endif - for(;i+1>1]; int lo=(int)(byte&0xF)-8, hi=(int)(byte>>4)-8; - a += xs[i]*(float)lo + xs[i+1]*(float)hi; } - if(i>1]; int lo=(int)(byte&0xF)-8; a += xs[i]*(float)lo; } - y[(int64_t)s*O+o]=a*sc; } } -} -/* y[S,O] = x[S,I] @ W^T with W int4 packed (2/byte) + per-GROUP scales (fmt=4). - * Same nibble math as matmul_i4, but the scale changes every `gs` elements along I. - * The accumulator resets at each group boundary: dot(x[grp], w[grp]) * scale[grp]. - * gs MUST be a multiple of 16 (the AVX2 vector width). */ -static void matmul_i4_grouped(float *y, const float *x, const uint8_t *q4, const float *scale, - int S, int I, int O, int gs){ - int rb=(I+1)/2; int ng=(I+gs-1)/gs; - #pragma omp parallel for schedule(static) - for(int o=0;oI) glen=I-base; - float sc=scl[g]; - int i=base; -#ifdef __AVX2__ - const __m128i m4=_mm_set1_epi8(0x0F); const __m256i b8=_mm256_set1_epi32(8); - __m256 acc=_mm256_setzero_ps(); - for(; i+16<=base+glen; i+=16){ __m128i by=_mm_loadl_epi64((const __m128i*)(w+(i>>1))); - __m128i lo=_mm_and_si128(by,m4),hi=_mm_and_si128(_mm_srli_epi16(by,4),m4); - __m128i nib=_mm_unpacklo_epi8(lo,hi); - __m256 w0=_mm256_cvtepi32_ps(_mm256_sub_epi32(_mm256_cvtepu8_epi32(nib),b8)); - __m256 w1=_mm256_cvtepi32_ps(_mm256_sub_epi32(_mm256_cvtepu8_epi32(_mm_srli_si128(nib,8)),b8)); - acc=_mm256_fmadd_ps(_mm256_loadu_ps(xs+i), w0, acc); - acc=_mm256_fmadd_ps(_mm256_loadu_ps(xs+i+8), w1, acc); } - a+=hsum256(acc)*sc; -#endif - /* scalar tail for the group remainder */ - for(; i>1]; - a+=(xs[i]*(float)((int)(byte&0xF)-8)+xs[i+1]*(float)((int)(byte>>4)-8))*sc; } - else { uint8_t byte=w[i>>1]; a+=xs[i]*(float)((int)(byte&0xF)-8)*sc; } - } - } - y[(int64_t)s*O+o]=a; - } - } -} + /* Fused gate+up for grouped int4 (fmt=4): computes both yg[S,O] and yu[S,O] from - * the same x[S,I], with per-group scales. One OpenMP dispatch covers both matrices, - * reading x once instead of twice — saves ~33% of expert-matmul time at decode. + * the same x[S,I], reading x once instead of twice — saves ~33% of expert-matmul time at decode. * The per-group scale logic matches matmul_i4_grouped exactly. */ static void matmul_i4_grouped_pair(float *yg, float *yu, const float *x, const uint8_t *qg, const float *sg, @@ -611,540 +512,6 @@ static void matmul_i4_grouped_pair(float *yg, float *yu, const float *x, } } } -/* Decode hot path for gate+up: same exact q4 dot products as matmul_i4, but one - * OpenMP dispatch covers both matrices. KTransformers uses persistent pools; - * this keeps colibri dependency-free while removing one team launch/expert. */ -static void matmul_i4_pair(float *yg, float *yu, const float *x, - const uint8_t *qg, const float *sg, - const uint8_t *qu, const float *su, int I, int O){ - int rb=(I+1)/2; - #pragma omp parallel for schedule(static) - for(int z=0;z<2*O;z++){ - int o=z>1))); - __m128i lo=_mm_and_si128(by,m4),hi=_mm_and_si128(_mm_srli_epi16(by,4),m4); - __m128i nib=_mm_unpacklo_epi8(lo,hi); - __m256 w0=_mm256_cvtepi32_ps(_mm256_sub_epi32(_mm256_cvtepu8_epi32(nib),b8)); - __m256 w1=_mm256_cvtepi32_ps(_mm256_sub_epi32(_mm256_cvtepu8_epi32(_mm_srli_si128(nib,8)),b8)); - acc=_mm256_fmadd_ps(_mm256_loadu_ps(x+i),w0,acc); - acc=_mm256_fmadd_ps(_mm256_loadu_ps(x+i+8),w1,acc); } - a=hsum256(acc); -#elif defined(__ARM_NEON) - const uint8x8_t m4=vdup_n_u8(0x0F); const int8x8_t b8=vdup_n_s8(8); - float32x4_t ac0=vdupq_n_f32(0),ac1=vdupq_n_f32(0); - for(;i+16<=I;i+=16){ uint8x8_t by=vld1_u8(w+(i>>1)); - uint8x8x2_t n=vzip_u8(vand_u8(by,m4),vshr_n_u8(by,4)); - int16x8_t w0=vmovl_s8(vsub_s8(vreinterpret_s8_u8(n.val[0]),b8)); - int16x8_t w1=vmovl_s8(vsub_s8(vreinterpret_s8_u8(n.val[1]),b8)); - ac0=vfmaq_f32(ac0,vld1q_f32(x+i),vcvtq_f32_s32(vmovl_s16(vget_low_s16(w0)))); - ac1=vfmaq_f32(ac1,vld1q_f32(x+i+4),vcvtq_f32_s32(vmovl_s16(vget_high_s16(w0)))); - ac0=vfmaq_f32(ac0,vld1q_f32(x+i+8),vcvtq_f32_s32(vmovl_s16(vget_low_s16(w1)))); - ac1=vfmaq_f32(ac1,vld1q_f32(x+i+12),vcvtq_f32_s32(vmovl_s16(vget_high_s16(w1)))); } - a=vaddvq_f32(vaddq_f32(ac0,ac1)); -#endif -#if defined(__AVX512F__) && defined(__AVX512BW__) - } -#endif - for(;i+1>1]; a+=x[i]*(float)((b&15)-8)+x[i+1]*(float)((b>>4)-8); } - if(i>1]&15)-8); - (z=2) non calcolano la STESSA funzione. Tre interruttori dipendono da S: il gate - * int4-IDOT (S>=g_i4s — asimmetrico proprio dove g_i4s>1), la fusione gate+up solo-S==1, - * e la soglia righe del GEMM Metal. Con SPEC_PIN=1 (default) ogni forward emesso mentre - * i draft del modello sono attivi resta sulla famiglia di kernel di S=1: draft e verifica - * coincidono per costruzione. Prefill e decode non speculativo sono intoccati. - * EN: MTP acceptance collapses when the draft (S=1) and verify (S>=2) forwards do not - * compute the SAME function. Three switches are S-dependent: the int4 IDOT gate - * (S>=g_i4s — asymmetric exactly on ISAs where g_i4s>1), the S==1-only gate+up fusion, - * and the Metal GEMM row threshold. SPEC_PIN=1 (default) pins every forward issued - * while model drafts are live to the platform's S=1 kernel family, so draft and verify - * agree by construction; prefill and non-speculative decode are untouched. - * SPEC_PIN=0 restores the S-dependent gates (A/B). */ -static int g_spec_pin=1; -static int g_spec_live=0; /* set by spec_decode while drafts are live */ -static inline int spec_pinned(void){ return g_spec_pin && g_spec_live; } -static void expert_gate_up(float *g,float *u,const float *x,QT *wg,QT *wu,int S){ - if(!g_no_fused_pair&&!spec_pinned()&&S==1&&wg->fmt==2&&wu->fmt==2&&wg->I==wu->I&&wg->O==wu->O) - matmul_i4_pair(g,u,x,wg->q4,wg->s,wu->q4,wu->s,wg->I,wg->O); - else if(!g_no_fused_pair&&S==1&&wg->fmt==4&&wu->fmt==4&&wg->I==wu->I&&wg->O==wu->O&&wg->gs==wu->gs) - matmul_i4_grouped_pair(g,u,x,wg->q4,wg->s,wu->q4,wu->s,S,wg->I,wg->O,wg->gs); - else { matmul_qt(g,x,wg,S); matmul_qt(u,x,wu,S); } -} -/* y[S,O] = x[S,I] @ W^T con W int2 impacchettato (4 valori/byte) + scala[O]. nibble 2-bit -> [-2,1]. */ -static void matmul_i2(float *y, const float *x, const uint8_t *q2, const float *scale, int S, int I, int O){ - int rb=(I+3)/4; - #pragma omp parallel for schedule(static) - for (int o=0;o>2))); /* 4 byte=16 valori */ - __m128i p0=_mm_and_si128(by,m2), p1=_mm_and_si128(_mm_srli_epi16(by,2),m2); - __m128i p2=_mm_and_si128(_mm_srli_epi16(by,4),m2), p3=_mm_and_si128(_mm_srli_epi16(by,6),m2); - __m128i lo=_mm_unpacklo_epi8(p0,p1), hi=_mm_unpacklo_epi8(p2,p3); - __m128i nib=_mm_unpacklo_epi16(lo,hi); /* 16 valori in ordine */ - __m256 w0=_mm256_cvtepi32_ps(_mm256_sub_epi32(_mm256_cvtepu8_epi32(nib),b2)); - __m256 w1=_mm256_cvtepi32_ps(_mm256_sub_epi32(_mm256_cvtepu8_epi32(_mm_srli_si128(nib,8)),b2)); - acc=_mm256_fmadd_ps(_mm256_loadu_ps(xs+i), w0, acc); - acc=_mm256_fmadd_ps(_mm256_loadu_ps(xs+i+8), w1, acc); } - a=hsum256(acc); -#elif defined(__ARM_NEON) - const uint8x8_t m2v=vdup_n_u8(3); const int8x8_t b2v=vdup_n_s8(2); - float32x4_t ac0=vdupq_n_f32(0), ac1=vdupq_n_f32(0); - for(;i+16<=I;i+=16){ uint32_t wd; memcpy(&wd, w+(i>>2), 4); /* 4 byte=16 valori */ - uint8x8_t by=vreinterpret_u8_u32(vdup_n_u32(wd)); - uint8x8x2_t z01=vzip_u8(vand_u8(by,m2v), vand_u8(vshr_n_u8(by,2),m2v)); - uint8x8x2_t z23=vzip_u8(vand_u8(vshr_n_u8(by,4),m2v), vshr_n_u8(by,6)); - uint16x4x2_t zz=vzip_u16(vreinterpret_u16_u8(z01.val[0]), vreinterpret_u16_u8(z23.val[0])); - int16x8_t w0=vmovl_s8(vsub_s8(vreinterpret_s8_u16(zz.val[0]),b2v)); /* 16 valori in ordine */ - int16x8_t w1=vmovl_s8(vsub_s8(vreinterpret_s8_u16(zz.val[1]),b2v)); - ac0=vfmaq_f32(ac0, vld1q_f32(xs+i), vcvtq_f32_s32(vmovl_s16(vget_low_s16(w0)))); - ac1=vfmaq_f32(ac1, vld1q_f32(xs+i+4), vcvtq_f32_s32(vmovl_s16(vget_high_s16(w0)))); - ac0=vfmaq_f32(ac0, vld1q_f32(xs+i+8), vcvtq_f32_s32(vmovl_s16(vget_low_s16(w1)))); - ac1=vfmaq_f32(ac1, vld1q_f32(xs+i+12), vcvtq_f32_s32(vmovl_s16(vget_high_s16(w1)))); } - a=vaddvq_f32(vaddq_f32(ac0,ac1)); -#endif - for(;i>2]; int sh=(i&3)*2; a += xs[i]*(float)((int)((byte>>sh)&3)-2); } - y[(int64_t)s*O+o]=a*sc; } } -} -/* ---- KERNEL INTERI (IDOT): attivazioni quantizzate a int8 per riga (absmax/127, - * stile Q8_0), prodotto scalare INTERO via maddubs/madd AVX2 — niente conversione - * f32 dei pesi nel ciclo caldo. ~2-3x sui matmul quantizzati; errore aggiunto ~0.3% - * RMS per matmul (attivazione int8), IDOT=0 torna al percorso f32 esatto. */ -#if defined(__AVX512VNNI__) && defined(__AVX512BW__) -#define IDOT_KERNEL "avx512-vnni" -#elif defined(__AVXVNNI__) && defined(__AVX2__) -#define IDOT_KERNEL "avx-vnni" -#elif defined(__AVX2__) -#define IDOT_KERNEL "avx2" -#elif defined(__ARM_NEON) && defined(__ARM_FEATURE_MATMUL_INT8) -#define IDOT_KERNEL "neon-i8mm" -#elif defined(__ARM_NEON) -#define IDOT_KERNEL "neon" -#elif defined(__VSX__) -#define IDOT_KERNEL "vsx" -#else -#define IDOT_KERNEL "scalar" -#endif -static int g_idot=1; -#if defined(__ARM_NEON) && defined(__ARM_FEATURE_DOTPROD) -static int g_i4s=1; /* SDOT presente: int4 IDOT conviene anche a S=1 (decode). Misurato - * su Apple M-series: +14%%, expert-matmul -16%%. EN: with SDOT, int4 - * IDOT pays even at S=1 (decode); measured on Apple M-series. */ -#elif defined(__VSX__) -static int g_i4s=1; /* POWER8 vec_msum: qui il fallback f32 e' SCALARE, quindi l'IDOT - * int4 conviene anche a S=1. Misurato su POWER8 S824 (vedi PR). - * EN: on VSX the f32 fallback is plain scalar C, so int4 IDOT - * pays even at S=1. Measured on a POWER8 S824 (see PR). */ -#else -static int g_i4s=2; /* senza SDOT / altrove: soglia originale (misura AVX2 dell'autore). - * EN: without SDOT / elsewhere: original threshold (author's AVX2). */ -#endif -static inline float qrow_i8(const float *x, int8_t *q, int I){ - float amax=0; for(int i=0;iamax)amax=a; } - float s=amax/127.f; if(s<1e-12f) s=1e-12f; float inv=1.f/s; - for(int i=0;i s32 directly, 64 bytes/iter, no 16-bit intermediate. - * AVX-512 has no vpsignb: |w| via abs, sign folded into x with a mask-negate - * (w==0 -> product 0 either way). |x|<=127 (qrow_i8), |w|<=128 as u8: each - * s32 lane adds <= 4*128*127, safe up to I=16384 like the AVX2 bound. */ - __m512i acc=_mm512_setzero_si512(); - for(;i+64<=I;i+=64){ - __m512i wv=_mm512_loadu_si512((const void*)(w+i)); - __m512i xv=_mm512_loadu_si512((const void*)(x+i)); - __mmask64 neg=_mm512_movepi8_mask(wv); - __m512i xs=_mm512_mask_sub_epi8(xv,neg,_mm512_setzero_si512(),xv); - acc=_mm512_dpbusd_epi32(acc,_mm512_abs_epi8(wv),xs); - } - sum=_mm512_reduce_add_epi32(acc); -#elif defined(__AVXVNNI__) && defined(__AVX2__) - /* AVX-VNNI 128-bit: vpdpbusd u8*s8 -> s32, 16 byte/iter. Stesso trucco del - * segno della variante 512-bit: |w| via abs, segno piegato in x con maschera - * (w==0 -> product 0). __AVX2__ serve per _mm_sign_epi8 / abs. */ - __m128i acc=_mm_setzero_si128(); - for(;i+16<=I;i+=16){ - __m128i wv=_mm_loadu_si128((const __m128i*)(w+i)); - __m128i xv=_mm_loadu_si128((const __m128i*)(x+i)); - __m128i xs=_mm_sign_epi8(xv,wv); /* x * sign(w); _mm_sign zona __AVX2__ */ - acc=_mm_dpbusd_epi32(acc,_mm_abs_epi8(wv),xs); - } - sum=hsum128_i32(acc); -#elif defined(__AVX2__) - __m256i acc=_mm256_setzero_si256(); const __m256i ones=_mm256_set1_epi16(1); - for(;i+32<=I;i+=32){ - __m256i wv=_mm256_loadu_si256((const __m256i*)(w+i)); - __m256i xv=_mm256_loadu_si256((const __m256i*)(x+i)); - __m256i p=_mm256_maddubs_epi16(_mm256_sign_epi8(wv,wv),_mm256_sign_epi8(xv,wv)); - acc=_mm256_add_epi32(acc,_mm256_madd_epi16(p,ones)); - } - sum=hsum256_i32(acc); -#elif defined(__ARM_NEON) - /* ARM: SDOT nativo se disponibile (Apple Silicon: sempre); altrimenti vmull/vpadal. - * Stesso bound anti-overflow del trucco AVX2: coppie <= 128*127*2 = 32512 < 32767. */ -#if defined(__ARM_FEATURE_DOTPROD) - /* 4 accumulatori indipendenti: SDOT ha latenza ~3-4 cicli, con un solo acc la - * catena seriale strozza il core a ~26 GB/s di pesi; con 4 lane indipendenti il - * dot diventa memory-bound (misurato su M4: 26 -> 63 GB/s per core, 2.4x). */ - int32x4_t a0=vdupq_n_s32(0),a1=vdupq_n_s32(0),a2=vdupq_n_s32(0),a3=vdupq_n_s32(0); - for(;i+64<=I;i+=64){ - a0=vdotq_s32(a0,vld1q_s8(w+i), vld1q_s8(x+i)); - a1=vdotq_s32(a1,vld1q_s8(w+i+16),vld1q_s8(x+i+16)); - a2=vdotq_s32(a2,vld1q_s8(w+i+32),vld1q_s8(x+i+32)); - a3=vdotq_s32(a3,vld1q_s8(w+i+48),vld1q_s8(x+i+48)); - } - int32x4_t acc=vaddq_s32(vaddq_s32(a0,a1),vaddq_s32(a2,a3)); - for(;i+16<=I;i+=16) acc=vdotq_s32(acc,vld1q_s8(w+i),vld1q_s8(x+i)); - sum=vaddvq_s32(acc); -#else - int32x4_t acc=vdupq_n_s32(0); - for(;i+16<=I;i+=16){ - int8x16_t wv=vld1q_s8(w+i), xv=vld1q_s8(x+i); - int16x8_t p=vmull_s8(vget_low_s8(wv),vget_low_s8(xv)); - p=vmlal_s8(p,vget_high_s8(wv),vget_high_s8(xv)); - acc=vpadalq_s16(acc,p); - } - sum=vaddvq_s32(acc); -#endif -#elif defined(__VSX__) - /* POWER8: vec_msum (s8 x u8 -> s32) somma i prodotti byte DIRETTAMENTE in lane - * s32, 16 byte/iter: il bound anti-saturazione a 16 bit di maddubs qui non serve. - * Stesso trucco del segno (|w| u8 per x*sign(w) s8), ma |w| via select+sub MODULO - * e non vec_abs: -128 deve diventare 128 u8, non saturare a 127. - * EN: vec_msum accumulates byte products straight into s32 lanes; |w| is built - * with a modulo subtract select instead of vec_abs so w=-128 wraps to 128 (u8) - * rather than saturating to 127. |x|<=127 from qrow_i8, so x negation is safe. */ - __vector signed int acc=vec_splats(0); - const __vector signed char vz=vec_splats((signed char)0); - for(;i+16<=I;i+=16){ - __vector signed char wv=vec_xl(0,(const signed char*)(w+i)); - __vector signed char xv=vec_xl(0,(const signed char*)(x+i)); - __vector __bool char neg=vec_cmplt(wv,vz); - __vector signed char xs=vec_sel(xv,vec_sub(vz,xv),neg); - __vector unsigned char wa=(__vector unsigned char)vec_sel(wv,vec_sub(vz,wv),neg); - acc=vec_msum(xs,wa,acc); - } - sum=vec_extract(acc,0)+vec_extract(acc,1)+vec_extract(acc,2)+vec_extract(acc,3); -#endif - for(;i int8 [-8,7] al volo, poi stesso trucco */ -static inline int32_t dot_i4i8(const uint8_t *w4, const int8_t *x, int I){ - int32_t sum=0; int i=0; -#if defined(__AVX512VNNI__) && defined(__AVX512BW__) - /* 32 bytes = 64 nibbles -> int8 in [-8,7], one vpdpbusd per 64 values. - * 256-bit unpack leaves values in per-128-lane order [0-15][32-47]/[16-31][48-63]; - * dot pairing is order-invariant, so permute x's 128-bit blocks to match - * instead of re-ordering w (one vpermq per iter, off the critical unpack path). */ - const __m256i m4v=_mm256_set1_epi8(0x0F); - const __m512i b8v=_mm512_set1_epi8(8); - const __m512i xidx=_mm512_setr_epi64(0,1,4,5,2,3,6,7); - __m512i acc=_mm512_setzero_si512(); - for(;i+64<=I;i+=64){ - __m256i by=_mm256_loadu_si256((const __m256i*)(w4+(i>>1))); - __m256i lo=_mm256_and_si256(by,m4v), hi=_mm256_and_si256(_mm256_srli_epi16(by,4),m4v); - __m256i z0=_mm256_unpacklo_epi8(lo,hi), z1=_mm256_unpackhi_epi8(lo,hi); - __m512i wv=_mm512_sub_epi8(_mm512_inserti64x4(_mm512_castsi256_si512(z0),z1,1),b8v); - __m512i xv=_mm512_permutexvar_epi64(xidx,_mm512_loadu_si512((const void*)(x+i))); - __mmask64 neg=_mm512_movepi8_mask(wv); - __m512i xs=_mm512_mask_sub_epi8(xv,neg,_mm512_setzero_si512(),xv); - acc=_mm512_dpbusd_epi32(acc,_mm512_abs_epi8(wv),xs); - } - sum=_mm512_reduce_add_epi32(acc); -#elif defined(__AVXVNNI__) && defined(__AVX2__) - /* AVX-VNNI 128-bit, int4: 16 byte = 32 nibble -> int8 [-8,7] in due half - * (n0/n1), ciascuno alimentato a un vpdpbusd da 16 byte. Stesso unpack - * 128-bit del ramo AVX2 sotto; 32 elementi/iter come li. */ - const __m128i m4=_mm_set1_epi8(0x0F); const __m128i b8=_mm_set1_epi8(8); - __m128i acc=_mm_setzero_si128(); - for(;i+32<=I;i+=32){ - __m128i by=_mm_loadu_si128((const __m128i*)(w4+(i>>1))); /* 16 byte = 32 nibble */ - __m128i lo=_mm_and_si128(by,m4), hi=_mm_and_si128(_mm_srli_epi16(by,4),m4); - __m128i n0=_mm_unpacklo_epi8(lo,hi), n1=_mm_unpackhi_epi8(lo,hi); /* nibble in ordine */ - __m128i w0=_mm_sub_epi8(n0,b8), w1=_mm_sub_epi8(n1,b8); - __m128i x0=_mm_loadu_si128((const __m128i*)(x+i)); - __m128i x1=_mm_loadu_si128((const __m128i*)(x+i+16)); - acc=_mm_dpbusd_epi32(acc,_mm_abs_epi8(w0),_mm_sign_epi8(x0,w0)); - acc=_mm_dpbusd_epi32(acc,_mm_abs_epi8(w1),_mm_sign_epi8(x1,w1)); - } - sum=hsum128_i32(acc); -#elif defined(__AVX2__) - const __m128i m4=_mm_set1_epi8(0x0F); const __m256i b8=_mm256_set1_epi8(8); - const __m256i ones=_mm256_set1_epi16(1); - __m256i acc=_mm256_setzero_si256(); - for(;i+32<=I;i+=32){ - __m128i by=_mm_loadu_si128((const __m128i*)(w4+(i>>1))); /* 16 byte = 32 nibble */ - __m128i lo=_mm_and_si128(by,m4), hi=_mm_and_si128(_mm_srli_epi16(by,4),m4); - __m128i n0=_mm_unpacklo_epi8(lo,hi), n1=_mm_unpackhi_epi8(lo,hi); /* in ordine */ - __m256i wv=_mm256_sub_epi8(_mm256_set_m128i(n1,n0),b8); - __m256i xv=_mm256_loadu_si256((const __m256i*)(x+i)); - __m256i p=_mm256_maddubs_epi16(_mm256_sign_epi8(wv,wv),_mm256_sign_epi8(xv,wv)); - acc=_mm256_add_epi32(acc,_mm256_madd_epi16(p,ones)); - } - sum=hsum256_i32(acc); -#elif defined(__ARM_NEON) - const uint8x16_t m4q=vdupq_n_u8(0x0F); const int8x16_t b8q=vdupq_n_s8(8); -#if defined(__ARM_FEATURE_DOTPROD) - /* 4 accumulatori indipendenti (vedi dot_i8i8): spezza la catena seriale su acc. - * Misurato su M4: 12.4 -> 29.9 GB/s di pesi per core (2.4x). */ - int32x4_t a0=vdupq_n_s32(0),a1=vdupq_n_s32(0),a2=vdupq_n_s32(0),a3=vdupq_n_s32(0); - for(;i+64<=I;i+=64){ - uint8x16_t byA=vld1q_u8(w4+(i>>1)), byB=vld1q_u8(w4+(i>>1)+16); - uint8x16x2_t zA=vzipq_u8(vandq_u8(byA,m4q), vshrq_n_u8(byA,4)); /* nibble in ordine */ - uint8x16x2_t zB=vzipq_u8(vandq_u8(byB,m4q), vshrq_n_u8(byB,4)); - a0=vdotq_s32(a0,vsubq_s8(vreinterpretq_s8_u8(zA.val[0]),b8q),vld1q_s8(x+i)); - a1=vdotq_s32(a1,vsubq_s8(vreinterpretq_s8_u8(zA.val[1]),b8q),vld1q_s8(x+i+16)); - a2=vdotq_s32(a2,vsubq_s8(vreinterpretq_s8_u8(zB.val[0]),b8q),vld1q_s8(x+i+32)); - a3=vdotq_s32(a3,vsubq_s8(vreinterpretq_s8_u8(zB.val[1]),b8q),vld1q_s8(x+i+48)); - } - int32x4_t acc=vaddq_s32(vaddq_s32(a0,a1),vaddq_s32(a2,a3)); - for(;i+32<=I;i+=32){ - uint8x16_t by=vld1q_u8(w4+(i>>1)); /* 16 byte = 32 nibble */ - uint8x16x2_t z=vzipq_u8(vandq_u8(by,m4q), vshrq_n_u8(by,4)); /* nibble in ordine */ - acc=vdotq_s32(acc,vsubq_s8(vreinterpretq_s8_u8(z.val[0]),b8q),vld1q_s8(x+i)); - acc=vdotq_s32(acc,vsubq_s8(vreinterpretq_s8_u8(z.val[1]),b8q),vld1q_s8(x+i+16)); - } - sum=vaddvq_s32(acc); -#else - int32x4_t acc=vdupq_n_s32(0); - for(;i+32<=I;i+=32){ - uint8x16_t by=vld1q_u8(w4+(i>>1)); /* 16 byte = 32 nibble */ - uint8x16x2_t z=vzipq_u8(vandq_u8(by,m4q), vshrq_n_u8(by,4)); /* nibble in ordine */ - int8x16_t w0=vsubq_s8(vreinterpretq_s8_u8(z.val[0]),b8q); - int8x16_t w1=vsubq_s8(vreinterpretq_s8_u8(z.val[1]),b8q); - int8x16_t x0=vld1q_s8(x+i), x1=vld1q_s8(x+i+16); - int16x8_t p=vmull_s8(vget_low_s8(w0),vget_low_s8(x0)); /* |w|<=8: nessun overflow */ - p=vmlal_s8(p,vget_high_s8(w0),vget_high_s8(x0)); - acc=vpadalq_s16(acc,p); - p=vmull_s8(vget_low_s8(w1),vget_low_s8(x1)); - p=vmlal_s8(p,vget_high_s8(w1),vget_high_s8(x1)); - acc=vpadalq_s16(acc,p); - } - sum=vaddvq_s32(acc); -#endif -#elif defined(__VSX__) - /* 16 byte = 32 nibble. vec_mergeh/vec_mergel su ppc64le (GCC) interallacciano come - * unpacklo/unpackhi x86 (verificato empiricamente su POWER8): i nibble escono in - * ordine di memoria. |w|<=8 dopo il -8, quindi stesso trucco del segno di dot_i8i8. - * EN: vec_mergeh/l on ppc64le interleave like x86 unpacklo/hi (verified on POWER8), - * so nibbles come out in memory order; then the same sign trick as dot_i8i8. */ - const __vector unsigned char m4v=vec_splats((unsigned char)0x0F); - const __vector unsigned char sh4=vec_splats((unsigned char)4); - const __vector signed char b8v=vec_splats((signed char)8); - const __vector signed char vz=vec_splats((signed char)0); - __vector signed int acc=vec_splats(0); - for(;i+32<=I;i+=32){ - __vector unsigned char by=vec_xl(0,w4+(i>>1)); /* 16 byte = 32 nibble */ - __vector unsigned char lo=vec_and(by,m4v), hi=vec_sr(by,sh4); - __vector signed char w0=vec_sub((__vector signed char)vec_mergeh(lo,hi),b8v); - __vector signed char w1=vec_sub((__vector signed char)vec_mergel(lo,hi),b8v); - __vector signed char x0=vec_xl(0,(const signed char*)(x+i)); - __vector signed char x1=vec_xl(0,(const signed char*)(x+i+16)); - __vector __bool char n0=vec_cmplt(w0,vz), n1=vec_cmplt(w1,vz); - acc=vec_msum(vec_sel(x0,vec_sub(vz,x0),n0), - (__vector unsigned char)vec_sel(w0,vec_sub(vz,w0),n0),acc); - acc=vec_msum(vec_sel(x1,vec_sub(vz,x1),n1), - (__vector unsigned char)vec_sel(w1,vec_sub(vz,w1),n1),acc); - } - sum=vec_extract(acc,0)+vec_extract(acc,1)+vec_extract(acc,2)+vec_extract(acc,3); -#endif - for(;i+1>1]; sum+=((int)(b&0xF)-8)*x[i]+((int)(b>>4)-8)*x[i+1]; } - if(i>1]; sum+=((int)(b&0xF)-8)*x[i]; } - return sum; -} -#if defined(__ARM_NEON) && defined(__ARM_FEATURE_MATMUL_INT8) -/* SMMLA (i8mm): vmmlaq_s32 vede ogni int8x16_t come matrice 2x8 row-major (byte 0-7 = - * riga 0, byte 8-15 = riga 1) e accumula C += A*B^T nel 2x2 int32: lane0=a0.b0, - * lane1=a0.b1, lane2=a1.b0, lane3=a1.b1. vcombine di due mezze-righe costruisce la - * matrice: A = due righe di peso (o,o+1), B = due righe di attivazione (s,s+1), quindi - * meta' traffico pesi e doppio lavoro per istruzione a S>=2. EN: vmmlaq_s32 treats each - * int8x16_t as a 2x8 row-major matrix and does C += A*B^T on a 2x2 int32 tile; vcombine - * of vget_low/high halves builds the 2-row register from two weight/activation rows. */ -static inline int32x4_t mm_tile16(int32x4_t acc, int8x16_t wo, int8x16_t wo1, - int8x16_t xs, int8x16_t xs1){ - acc=vmmlaq_s32(acc, vcombine_s8(vget_low_s8(wo), vget_low_s8(wo1)), - vcombine_s8(vget_low_s8(xs), vget_low_s8(xs1))); - return vmmlaq_s32(acc, vcombine_s8(vget_high_s8(wo), vget_high_s8(wo1)), - vcombine_s8(vget_high_s8(xs), vget_high_s8(xs1))); -} -static void matmul_q_idot_mm(float *y, const int8_t *xq, const float *sx, const int8_t *q, - const float *scale, int S, int I, int O){ - #pragma omp parallel for schedule(static) - for(int o=0;o<(O&~1);o+=2){ - const int8_t *wo=q+(int64_t)o*I, *wo1=q+(int64_t)(o+1)*I; - float sc0=scale[o], sc1=scale[o+1]; - for(int s=0;s<(S&~1);s+=2){ - const int8_t *xs=xq+(int64_t)s*I, *xs1=xq+(int64_t)(s+1)*I; - /* 4 accumulatori indipendenti: una sola catena vmmla e' latency-bound. - * EN: 4 independent accumulators; a single vmmla chain is latency-bound. */ - int32x4_t a0=vdupq_n_s32(0),a1=vdupq_n_s32(0),a2=vdupq_n_s32(0),a3=vdupq_n_s32(0); int i=0; - for(;i+64<=I;i+=64){ - a0=mm_tile16(a0,vld1q_s8(wo+i), vld1q_s8(wo1+i), vld1q_s8(xs+i), vld1q_s8(xs1+i)); - a1=mm_tile16(a1,vld1q_s8(wo+i+16),vld1q_s8(wo1+i+16),vld1q_s8(xs+i+16),vld1q_s8(xs1+i+16)); - a2=mm_tile16(a2,vld1q_s8(wo+i+32),vld1q_s8(wo1+i+32),vld1q_s8(xs+i+32),vld1q_s8(xs1+i+32)); - a3=mm_tile16(a3,vld1q_s8(wo+i+48),vld1q_s8(wo1+i+48),vld1q_s8(xs+i+48),vld1q_s8(xs1+i+48)); - } - for(;i+16<=I;i+=16) - a0=mm_tile16(a0,vld1q_s8(wo+i),vld1q_s8(wo1+i),vld1q_s8(xs+i),vld1q_s8(xs1+i)); - int32x4_t acc=vaddq_s32(vaddq_s32(a0,a1),vaddq_s32(a2,a3)); - int32_t d00=vgetq_lane_s32(acc,0), d01=vgetq_lane_s32(acc,1); - int32_t d10=vgetq_lane_s32(acc,2), d11=vgetq_lane_s32(acc,3); - for(;i>1)), byo1=vld1q_u8(wo1+(i>>1)); - uint8x16_t cyo=vld1q_u8(wo+(i>>1)+16), cyo1=vld1q_u8(wo1+(i>>1)+16); - uint8x16x2_t zo =vzipq_u8(vandq_u8(byo, m4q), vshrq_n_u8(byo, 4)); - uint8x16x2_t zo1=vzipq_u8(vandq_u8(byo1,m4q), vshrq_n_u8(byo1,4)); - uint8x16x2_t ko =vzipq_u8(vandq_u8(cyo, m4q), vshrq_n_u8(cyo, 4)); - uint8x16x2_t ko1=vzipq_u8(vandq_u8(cyo1,m4q), vshrq_n_u8(cyo1,4)); - a0=mm_tile16(a0, vsubq_s8(vreinterpretq_s8_u8(zo.val[0]),b8q), - vsubq_s8(vreinterpretq_s8_u8(zo1.val[0]),b8q), - vld1q_s8(xs+i), vld1q_s8(xs1+i)); - a1=mm_tile16(a1, vsubq_s8(vreinterpretq_s8_u8(zo.val[1]),b8q), - vsubq_s8(vreinterpretq_s8_u8(zo1.val[1]),b8q), - vld1q_s8(xs+i+16), vld1q_s8(xs1+i+16)); - a2=mm_tile16(a2, vsubq_s8(vreinterpretq_s8_u8(ko.val[0]),b8q), - vsubq_s8(vreinterpretq_s8_u8(ko1.val[0]),b8q), - vld1q_s8(xs+i+32), vld1q_s8(xs1+i+32)); - a3=mm_tile16(a3, vsubq_s8(vreinterpretq_s8_u8(ko.val[1]),b8q), - vsubq_s8(vreinterpretq_s8_u8(ko1.val[1]),b8q), - vld1q_s8(xs+i+48), vld1q_s8(xs1+i+48)); - } - for(;i+32<=I;i+=32){ - uint8x16_t byo=vld1q_u8(wo+(i>>1)), byo1=vld1q_u8(wo1+(i>>1)); - uint8x16x2_t zo =vzipq_u8(vandq_u8(byo, m4q), vshrq_n_u8(byo, 4)); - uint8x16x2_t zo1=vzipq_u8(vandq_u8(byo1,m4q), vshrq_n_u8(byo1,4)); - a0=mm_tile16(a0, vsubq_s8(vreinterpretq_s8_u8(zo.val[0]),b8q), - vsubq_s8(vreinterpretq_s8_u8(zo1.val[0]),b8q), - vld1q_s8(xs+i), vld1q_s8(xs1+i)); - a1=mm_tile16(a1, vsubq_s8(vreinterpretq_s8_u8(zo.val[1]),b8q), - vsubq_s8(vreinterpretq_s8_u8(zo1.val[1]),b8q), - vld1q_s8(xs+i+16), vld1q_s8(xs1+i+16)); - } - int32x4_t acc=vaddq_s32(vaddq_s32(a0,a1),vaddq_s32(a2,a3)); - int32_t d00=vgetq_lane_s32(acc,0), d01=vgetq_lane_s32(acc,1); - int32_t d10=vgetq_lane_s32(acc,2), d11=vgetq_lane_s32(acc,3); - for(;i+1>1], bo1=wo1[i>>1]; - int a0=(int)(bo&0xF)-8, a1=(int)(bo>>4)-8, b0=(int)(bo1&0xF)-8, b1=(int)(bo1>>4)-8; - int u0=xs[i],u1=xs[i+1],v0=xs1[i],v1=xs1[i+1]; - d00+=a0*u0+a1*u1; d01+=a0*v0+a1*v1; d10+=b0*u0+b1*u1; d11+=b0*v0+b1*v1; } - if(i>1], bo1=wo1[i>>1]; - int a0=(int)(bo&0xF)-8, b0=(int)(bo1&0xF)-8; - d00+=a0*xs[i]; d01+=a0*xs1[i]; d10+=b0*xs[i]; d11+=b0*xs1[i]; } - y[(int64_t)s*O+o] =(float)d00*sc0*sx[s]; - y[(int64_t)s*O+(o+1)] =(float)d10*sc1*sx[s]; - y[(int64_t)(s+1)*O+o] =(float)d01*sc0*sx[s+1]; - y[(int64_t)(s+1)*O+(o+1)]=(float)d11*sc1*sx[s+1]; - } - if(S&1){ int s=S-1; const int8_t *xs=xq+(int64_t)s*I; - y[(int64_t)s*O+o] =(float)dot_i4i8(wo, xs,I)*sc0*sx[s]; - y[(int64_t)s*O+(o+1)]=(float)dot_i4i8(wo1,xs,I)*sc1*sx[s]; } - } - if(O&1){ int o=O-1; const uint8_t *w=q4+(int64_t)o*rb; float sc=scale[o]; - #pragma omp parallel for schedule(static) - for(int s=0;s=2){ matmul_q_idot_mm(y,xq,sx,q,scale,S,I,O); return; } -#endif - #pragma omp parallel for schedule(static) - for(int o=0;o=2){ matmul_i4_idot_mm(y,xq,sx,q4,scale,S,I,O); return; } -#endif - #pragma omp parallel for schedule(static) - for(int o=0;og_qscratch.xq_cap){ - int8_t *p=realloc(g_qscratch.xq,xn); - if(!p){ fprintf(stderr,"OOM quant scratch\n"); exit(1); } - g_qscratch.xq=p; g_qscratch.xq_cap=xn; - } - if(sn>g_qscratch.sx_cap){ - float *p=realloc(g_qscratch.sx,sn*sizeof(float)); - if(!p){ fprintf(stderr,"OOM quant scales\n"); exit(1); } - g_qscratch.sx=p; g_qscratch.sx_cap=sn; - } - *xq=g_qscratch.xq; *sx=g_qscratch.sx; -} /* allow_idot=0: forza il kernel int4/int8 ESATTO (attivazioni f32). Serve alle proiezioni di * attenzione: sono sensibili alla quantizzazione int8 delle attivazioni dell'IDOT. Misurato su @@ -1153,8 +520,6 @@ static void quant_scratch(size_t xn, size_t sn, int8_t **xq, float **sx){ * EN: allow_idot=0 forces the EXACT int4/int8 kernel (f32 activations). The attention * projections need it: IDOT's int8 activation quantization costs +0.117 nats/token there * (~+12% perplexity), measured. Every other prefill matmul keeps IDOT as before. */ -static void matmul_qt_ex(float *y, const float *x, QT *w, int S, int allow_idot); -static void matmul_qt(float *y, const float *x, QT *w, int S){ matmul_qt_ex(y,x,w,S,1); } static void matmul_qt_ex(float *y, const float *x, QT *w, int S, int allow_idot){ #ifdef COLI_METAL if(g_metal_enabled && S>=g_metal_gemm_min && !spec_pinned() && (w->fmt==1||w->fmt==2) && !omp_in_parallel()){