diff --git a/build/src/build.bqn b/build/src/build.bqn index a6306709..d39b4124 100755 --- a/build/src/build.bqn +++ b/build/src/build.bqn @@ -576,7 +576,7 @@ cachedBin‿linkerCache ← { "xa"‿"src/builtins/arithd.c"‿"dyarith", "xa"‿"src/builtins/cmp.c"‿"cmp", "xa"‿"src/builtins/squeeze.c"‿"squeeze" "x."‿"src/builtins/select.c"‿"select", "x."‿"src/builtins/fold.c"‿"fold", "x."‿"src/builtins/scan.c"‿"scan" "x."‿"src/builtins/scan.c"‿"neq", "x."‿"src/builtins/slash.c"‿"slash", "x."‿"src/builtins/slash.c"‿"constrep" - "x."‿"src/builtins/transpose.c"‿"transpose" + "xa"‿"src/builtins/transpose.c"‿"transpose" ⟩ objs ← ⟨⟩ diff --git a/src/builtins/md1.c b/src/builtins/md1.c index 4d8ff0ad..76b23a71 100644 --- a/src/builtins/md1.c +++ b/src/builtins/md1.c @@ -322,6 +322,12 @@ static NOINLINE B match_cells(bool ne, B w, B x, ur wr, ur xr, usz len) { return r; } +B transp_c2(B, B, B); +static B transp_cells(ur ax, B x) { + i8* wp; B w=m_i8arrv(&wp, 2); wp[0]=0; wp[1]=ax; + return C2(transp, w, x); +} + B shape_c1(B, B); B fold_rows(Md1D* d, B x); // From fold.c B cell_c1(Md1D* d, B x) { B f = d->f; @@ -362,6 +368,7 @@ B cell_c1(Md1D* d, B x) { B f = d->f; B xf = getFillR(x); if (!noFill(xf)) return shift_cells(xf, x, TI(x,elType), rtid); } + if (rtid==n_transp) return xr<=2? x : transp_cells(xr-1, x); if (TY(f) == t_md1D) { Md1D* fd = c(Md1D,f); u8 rtid = fd->m1->flags-1; @@ -425,6 +432,7 @@ B cell_c2(Md1D* d, B w, B x) { B f = d->f; } if (rtid==n_take && xr>1 && isF64(w)) return takedrop_highrank(1, m_hVec2(m_f64(SH(x)[0]), w), x); if (rtid==n_drop && xr>1 && isF64(w)) return takedrop_highrank(0, m_hVec2(m_f64(0), w), x); + if (rtid==n_transp && q_usz(w)) { usz a=o2sG(w); if (a @@ -10,25 +46,83 @@ #endif #endif -#define TRANSPOSE_LOOP( DST, SRC, W, H) PLAINLOOP for(usz y=0;y< H;y++) NOVECTORIZE for(usz x=0;x< W;x++) DST[x*H+y] = SRC[xi++] -#define TRANSPOSE_BLOCK(DST, SRC, BW, BH, W, H) PLAINLOOP for(usz y=0;ya; + for (usz i=0; ia, xsh+1, xr-1); @@ -47,116 +140,300 @@ B transp_c1(B t, B x) { return taga(r); } usz w = xsh[1] * shProd(xsh, 2, xr); - if (w==1) goto no_reorder; - Arr* r; - usz xi = 0; - u8 xe = TI(x,elType); - bool toBit = false; - if (h==2) { - if (xe==el_B) { - B* xp = TO_BPTR(x); - B* x0 = xp; B* x1 = x0+w; - HArr_p rp = m_harrUp(ia); - for (usz i=0; ia; - for (usz i=0; i1) thrM("⍉: 𝕨 must have rank at most 1"); + wia = IA(w); + if (wia==0) { decG(w); return isArr(x)? x : m_atomUnit(x); } + } + ur xr; + if (isAtm(x) || (xr=RNK(x))=xr) thrF("⍉: Axis %s does not exist (%i≡=𝕩)", a, xr); + if (a==xr-1) { TFREE(alloc); return C1(transp, x); } + p[0] = a; + } else { + SGetU(w) + for (usz i=0; i=xr) thrF("⍉: Axis %s does not exist (%i≡=𝕩)", a, xr); + p[i] = a; + } + decG(w); + } + + B r; + + // Compute shape for the given axes + usz* xsh = SH(x); + usz* rsh = (usz*)(p + xr); // Length xr + usz dup = 0, max = 0, id = 0; + usz no_sh = -(usz)1; + for (usz j=0; jmax? j : max; + if (xl= rr) thrF("⍉: Skipped result axis"); + if (wia 1)) { // Not all duplicates + sh = m_shArr(rr); + shcpy(sh->a, rsh, rr); + } + + // Empty result + if (IA(x) == 0) { + Arr* ra = m_fillarrpEmpty(getFillQ(x)); + shSet(ra, rr, sh); + decG(x); + r = taga(ra); goto ret; + } + + // Add up stride for each axis + ur na = max + 1; // Number of result axes that moved + usz* st = rsh + xr; // Length na + for (usz j=0; j=1 */ \ + ur a = a0; \ + usz str = st[a]; \ + for (usz k=0; k= (32*8) >> xlw) { // cell >= 32 bytes + usz ria = csz * shProd(rsh, 0, na); + MAKE_MUT_INIT(rm, ria, xe); MUTG_INIT(rm); + AXIS_LOOP(na, csz, mut_copyG(rm, i, x, j, csz)); + Arr* ra = mut_fp(rm); + shSet(ra, rr, sh); + r = withFill(taga(ra), getFillQ(x)); + decG(x); goto ret; + } + #undef AXIS_LOOP + if ((csz & (csz-1))==0 && csz<=64>>xlw && csz<=8 // CPU-sized cells + && xe!=el_B && na>=2) { + // If some result axis has stride 1 (guaranteed if dup==0), then it + // corresponds to the last argument axis and we have a strided + // transpose swapping that with the last result axis + usz rai = na-1; + usz xai=rai; while (st[--xai]!=1) if (xai==0) goto skip_2d; + if (rsh[xai]*rsh[rai] < (256*8) >> xlw) goto skip_2d; + TranspFn tran = transposeFns[CTZ(csz<>= 3-xlw; + else if (xlw>3) for (usz i=0; ia[0] = IA(x)/csz; + zsh->a[1] = csz; + Arr* z = TI(x,slice)(x, 0, IA(x)); + arr_shSetU(z, 2, zsh); + x = taga(z); + // (+⌜´st×⟜↕¨rsh)⊏⥊𝕩 + B ind = bi_N; + for (ur k=na; k--; ) { + B v = C2(mul, m_usz(st[k]/csz), C1(ud, m_f64(rsh[k]))); + if (q_N(ind)) ind = v; + else ind = M1C2(tbl, add, v, ind); + } + r = C2(select, ind, x); + Arr* ra = cpyWithShape(r); r = taga(ra); + if (rr>1) arr_shReplace(ra, rr, sh); + else { decSh((Value*)ra); arr_shVec(ra); } + + ret:; + TFREE(alloc); + return r; +} + B transp_im(B t, B x) { if (isAtm(x)) thrM("⍉⁼: 𝕩 must not be an atom"); - if (RNK(x)<=2) return transp_c1(t, x); - return def_fn_im(bi_transp, x); + ur xr = RNK(x); + if (xr<=1) return x; + + usz ia = IA(x); + usz* xsh = SH(x); + usz w = xsh[xr-1]; + if (ia==0 || w==1 || w==ia /*h==1*/) { + Arr* r = cpyWithShape(x); + ShArr* sh = m_shArr(xr); + sh->a[0] = w; + shcpy(sh->a+1, xsh, xr-1); + arr_shReplace(r, xr, sh); + return taga(r); + } + usz h = xsh[0] * shProd(xsh, 1, xr-1); + + Arr* r = transpose_noshape(&x, ia, w, h); + + usz* rsh = arr_shAlloc(r, xr); + rsh[0] = w; + if (xr==2) rsh[1] = h; else shcpy(rsh+1, SH(x), xr-1); + decG(x); return taga(r); } -B transp_uc1(B t, B o, B x) { return transp_im(m_f64(0), c1(o, transp_c1(t, x))); } +B transp_uc1(B t, B o, B x) { + return transp_im(m_f64(0), c1(o, transp_c1(t, x))); +} + +// Consumes w; return bi_N if w contained duplicates +static B invert_transp_w(B w, ur xr) { + if (isAtm(w)) { + if (xr<1) thrM("⍉⁼: Length of 𝕨 must be at most rank of 𝕩"); + usz a=o2s(w); + if (a>=xr) thrF("⍉⁼: Axis %s does not exist (%i≡=𝕩)", a, xr); + i32* wp; w = m_i32arrv(&wp, a); + PLAINLOOP for (usz i=0; i1) thrM("⍉⁼: 𝕨 must have rank at most 1"); + usz wia = IA(w); + if (wia==0) return w; + if (xr=xr) thrF("⍉⁼: Axis %s does not exist (%i≡=𝕩)", a, xr); + if (p[a]!=xr) { TFREE(p); decG(w); return bi_N; } // Handled by caller + max = a>max? a : max; + p[a] = i; + } + decG(w); + usz n = max+1; + i32* wp; w = m_i32arrv(&wp, n); + for (usz i=0, j=wia; iuc1 = transp_uc1; c(BFn,bi_transp)->im = transp_im; + c(BFn,bi_transp)->ix = transp_ix; + c(BFn,bi_transp)->uc1 = transp_uc1; + c(BFn,bi_transp)->ucw = transp_ucw; } diff --git a/src/load.c b/src/load.c index f6485267..7fbdf0be 100644 --- a/src/load.c +++ b/src/load.c @@ -107,7 +107,7 @@ B comp_currSrc; B comp_currRe; // ⟨REPL mode ⋄ scope ⋄ compiler ⋄ runtime ⋄ glyphs ⋄ sysval names ⋄ sysval values⟩ B rt_undo, rt_select, rt_slash, rt_insert, rt_depth, - rt_group, rt_under, rt_find, rt_transp; + rt_group, rt_under, rt_find; Block* load_compObj(B x, B src, B path, Scope* sc) { // consumes x,src SGet(x) usz xia = IA(x); @@ -442,7 +442,6 @@ void load_init() { // very last init function rt_group = Get(rtObjRaw, n_group ); rt_under = Get(rtObjRaw, n_under ); rt_find = Get(rtObjRaw, n_find ); - rt_transp = Get(rtObjRaw, n_transp ); rt_depth = Get(rtObjRaw, n_depth ); rt_insert = Get(rtObjRaw, n_insert ); @@ -485,7 +484,7 @@ void load_init() { // very last init function } load_rtObj = frtObj; load_compArg = m_hVec2(load_rtObj, incG(bi_sys)); - rt_select=rt_slash=rt_group=rt_find=rt_transp=rt_invFnReg=rt_invFnSwap = incByG(bi_invalidFn, 7); + rt_select=rt_slash=rt_group=rt_find=rt_invFnReg=rt_invFnSwap = incByG(bi_invalidFn, 7); rt_undo=rt_insert = incByG(bi_invalidMd1, 2); rt_under=rt_depth = incByG(bi_invalidMd2, 2); rt_invFnRegFn=rt_invFnSwapFn = invalidFn_c1; @@ -497,7 +496,6 @@ void load_init() { // very last init function gc_add(rt_group); gc_add(rt_under); gc_add(rt_find); - gc_add(rt_transp); gc_add(rt_depth); gc_add(rt_insert); diff --git a/src/singeli/src/avx.singeli b/src/singeli/src/avx.singeli index 4a580745..b5d874ac 100644 --- a/src/singeli/src/avx.singeli +++ b/src/singeli/src/avx.singeli @@ -1,12 +1,3 @@ -def w256{T} = 0 -def w256{T & isvec{T}} = width{T}==256 -def w256{T,w} = 0 -def w256{T,w & w256{T}} = elwidth{T}==w -def w256i = genchk{w256, {T} => isint{T}} -def w256s = genchk{w256, {T} => issigned{T}} -def w256u = genchk{w256, {T} => isunsigned{T}} -def w256f = genchk{w256, {T} => isfloat{T}} - def v2i{x:T & w256{T}} = [32]u8 ~~ x # for compact casting for the annoying intrinsic type system def v2f{x:T & w256{T}} = [8]f32 ~~ x def v2d{x:T & w256{T}} = [4]f64 ~~ x diff --git a/src/singeli/src/base.singeli b/src/singeli/src/base.singeli index 2eb6e795..4bd83f18 100644 --- a/src/singeli/src/base.singeli +++ b/src/singeli/src/base.singeli @@ -11,15 +11,6 @@ def exportT{name, fs} = { v:*type{tupsel{0,fs}} = fs; export{name, v} } def elwidth{T} = width{eltype{T}} -def genchk{B, F} = { - def r{T} = 0 - def r{T & B{T}} = F{eltype{T}} - def r{T,w} = 0 - def r{T,w & B{T}} = F{eltype{T}} & (elwidth{T}==w) - def r{T & ~isvec{T}} = 0 - r -} - # ceiling divide def cdiv{a,b} = (a+b-1)/b @@ -61,6 +52,34 @@ def anyInt{x} = 0 def anyInt{x & knum{x}} = (x>>0) == x def anyInt{x & isreg{x}|isconst{x}} = isint{x} +# vector width/type checks +def w128{T} = 0 +def w128{T & isvec{T}} = width{T}==128 +def w128{T,w} = 0 +def w128{T,w & w128{T}} = elwidth{T}==w +def w256{T} = 0 +def w256{T & isvec{T}} = width{T}==256 +def w256{T,w} = 0 +def w256{T,w & w256{T}} = elwidth{T}==w + +# width+type checks +def genchk{B, F} = { + def r{T} = 0 + def r{T & B{T}} = F{eltype{T}} + def r{T,w} = 0 + def r{T,w & B{T}} = F{eltype{T}} & (elwidth{T}==w) + def r{T & ~isvec{T}} = 0 + r +} +def w128i = genchk{w128, {T} => isint{T}} +def w128s = genchk{w128, {T} => issigned{T}} +def w128u = genchk{w128, {T} => isunsigned{T}} +def w128f = genchk{w128, {T} => isfloat{T}} +def w256i = genchk{w256, {T} => isint{T}} +def w256s = genchk{w256, {T} => issigned{T}} +def w256u = genchk{w256, {T} => isunsigned{T}} +def w256f = genchk{w256, {T} => isfloat{T}} + def trunc{T, x:U & isint{T} & isint{U} & T<=U} = emit{T, '', x} def trunc{T, x & knum{x}} = cast{T, x} diff --git a/src/singeli/src/sse3.singeli b/src/singeli/src/sse3.singeli index 744919ef..47a1401d 100644 --- a/src/singeli/src/sse3.singeli +++ b/src/singeli/src/sse3.singeli @@ -1,12 +1,3 @@ -def w128{T} = 0 -def w128{T & isvec{T}} = width{T}==128 -def w128{T,w} = 0 -def w128{T,w & w128{T}} = elwidth{T}==w -def w128i = genchk{w128, {T} => isint{T}} -def w128s = genchk{w128, {T} => issigned{T}} -def w128u = genchk{w128, {T} => isunsigned{T}} -def w128f = genchk{w128, {T} => isfloat{T}} - def v2i{x:T & w128{T}} = [16]u8 ~~ x # for compact casting for the annoying intrinsic type system def v2f{x:T & w128{T}} = [4]f32 ~~ x def v2d{x:T & w128{T}} = [2]f64 ~~ x diff --git a/src/singeli/src/transpose.singeli b/src/singeli/src/transpose.singeli index bbdb9057..72050c9e 100644 --- a/src/singeli/src/transpose.singeli +++ b/src/singeli/src/transpose.singeli @@ -65,16 +65,8 @@ def for_mult_max{k, m}{vars,begin,end,block} = { } } -fn transpose{T, k, kh}(r0:*void, x0:*void, w:u64, h:u64) : void = { - # Scalar transpose defined in C - def ts = if (T==i8) 'i8' else if (T==i16) 'i16' else if (T==i32) 'i32' else 'i64' - def call_base{...a} = emit{void, merge{'base_transpose_',ts}, ...a, w, h} - - rp:*T = *T~~r0 - xp:*T = *T~~x0 - if (w 2*k or h&(line_elts-1) != 0) { + if (line_elts > 2*k or h&(line_elts-1) != 0 or h != hs) { ho := h%k # Effective height, like we for w he := h; if (use_overlap{ho}) he += k - ho @@ -96,7 +88,7 @@ fn transpose{T, k, kh}(r0:*void, x0:*void, w:u64, h:u64) : void = { # Main transpose @for_mult_max{kh, h-kh} (y to he) { @for_mult_max{k, wm} (x to we) { - kernel{...at{x,y}, k, kh, w, h} + kernel{...at{x,y}, k, kh, ws, hs} } } # Half-row(s) for non-square i16 case @@ -106,12 +98,12 @@ fn transpose{T, k, kh}(r0:*void, x0:*void, w:u64, h:u64) : void = { @for (yi to n) { y:u64 = 0; if (yi == n-1) y = h - e @for_mult_max{k, wm} (x to we) { - kernel{...at{x,y}, k, k, w, h} + kernel{...at{x,y}, k, k, ws, hs} } } } # Base transpose used if overlap wasn't - if (ho!=0 and he==h) { hs := h-ho; call_base{rp+hs, xp+w*hs, w, ho} } + if (ho!=0 and he==h) { hd := h-ho; call_base{rp+hd, xp+ws*hd, w, ho} } } else { # Result rows are aligned with each other so it's possible to # write a full cache line at a time @@ -125,49 +117,67 @@ fn transpose{T, k, kh}(r0:*void, x0:*void, w:u64, h:u64) : void = { def vt{i} = transpose_square{VT, k, each{loadx, k*i + iota{k}}} each{tup, ...each{vt, iota{line_vecs}}} } - ro := tail{6, -u64~~r0} / (width{T}/8) # Offset to align within cache line; assume elt-aligned - wh := w*h + ro := tail{6, -u64~~rp} / (width{T}/8) # Offset to align within cache line; assume elt-aligned + wh := ws*h yn := h if (ro != 0) { ra := line_elts - ro y := h - ra - rpe := rp + y + (w-1)*h # Cache aligned + rpe := rp + y + (w-1)*hs # Cache aligned # Part of first and last result row aren't covered by the split loop - def trtail{dst, src, len} = @for (i to len) store{dst, i, load{src, w*i}} + def trtail{dst, src, len} = @for (i to len) store{dst, i, load{src, ws*i}} trtail{rp, xp, ro} - trtail{rpe, xp + y*w + w-1, ra} + trtail{rpe, xp + y*ws + w-1, ra} # Transpose first few rows and last few rows together @for_mult_max{k, wm} (x to we) { {xpo,rpo} := at{x, y} - o := w*y + x + o := ws*y + x def loadx{_} = { l:=load{*VT~~(xp+o)} - o+=w; if (o>wh-k) o -= wh-1 # Jump from last source row to first, shifting right 1 + o+=ws; if (o>wh-k) o -= wh-1 # Jump from last source row to first, shifting right 1 l } def rls = get_lines{loadx} # 4 rows of 2 vectors each - each{{i,v} => {p:=rpo+i*h; if (i<3 or p {p:=rpo+i*hs; if (i<3 or p load{*VT~~(xpo+i*w), 0}} - each{{i,v} => store_line{*VT~~(rpo+i*h), v}, iota{k}, rls} + def rls = get_lines{{i} => load{*VT~~(xpo+i*ws), 0}} + each{{i,v} => store_line{*VT~~(rpo+i*hs), v}, iota{k}, rls} } } } - if (we==w) @for(ws from w-wo to w) { - xpo:=xp+ws; rpo:=rp+h*ws - @for (i to h) store{rpo, i, load{xpo, w*i}} + if (we==w) @for(wd from w-wo to w) { + xpo:=xp+wd; rpo:=rp+hs*wd + @for (i to h) store{rpo, i, load{xpo, ws*i}} + } +} + +fn transpose{T, k, kh}(r0:*void, x0:*void, w:u64, h:u64, ws:u64, hs:u64) : void = { + # Scalar transpose defined in C + def ts = if (T==i8) 'i8' else if (T==i16) 'i16' else if (T==i32) 'i32' else 'i64' + def call_base{...a} = emit{void, merge{'transpose_',ts}, ...a, ws, hs} + + rp:*T = *T~~r0 + xp:*T = *T~~x0 + if (hasarch{'X86_64'} and w>=k and h>=k) { + transpose_with_kernel{T, k, kh, call_base, rp, xp, w, h, ws, hs} + } else { + if (h==2 and h==hs) @for (x0 in xp, x1 in xp+ws over i to w) { store{rp, i*2, x0}; store{rp, i*2+1, x1} } + else if (w==2 and w==ws) @for (r0 in rp, r1 in rp+hs over i to h) { r0 = load{xp, i*2}; r1 = load{xp, i*2+1} } + else call_base{rp, xp, w, h} } } def transpose{T, k} = transpose{T, k, k} -export{'simd_transpose_i8', transpose{i8 , 16}} -export{'simd_transpose_i16', transpose{i16, 8, 16}} -export{'simd_transpose_i32', transpose{i32, 8}} -export{'simd_transpose_i64', transpose{i64, 4}} +exportT{'simd_transpose', tup{ + transpose{i8 , 16}, + transpose{i16, 8, 16}, + transpose{i32, 8}, + transpose{i64, 4} +}}