From 432b4eaaa65da38ddf5d68185d624aa255d45623 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Thu, 23 Mar 2023 13:41:06 -0400 Subject: [PATCH 01/26] Transpose implementation comments --- src/builtins/transpose.c | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 73b39f1c..3cadb5be 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -1,3 +1,26 @@ +// Transpose and Reorder Axes (⍉) + +// Transpose +// One length-2 axis: dedicated code +// Boolean: pdep for height 2; pext for width 2 +// SHOULD use a generic implementation if BMI2 not present +// SHOULD optimize other short lengths with pdep/pext and shuffles +// Boolean 𝕩: convert to integer +// SHOULD have bit matrix transpose kernel +// CPU sizes: native or SIMD code +// Large SIMD kernels used when they fit, overlapping for odd sizes +// i8: 16Γ—16; i16: 16Γ—8; i32: 8Γ—8; f64: 4Γ—4 +// COULD use half-width or smaller kernels to improve odd sizes +// Scalar transpose or loop used for overhang of 1 + +// Reorder Axes: self-hosted runtime (based on +⌜ and ⊏, not that slow) + +// Transpose inverse ⍉⁼ +// Same as ⍉ for a rank ≀2 argument +// SHOULD share data movement with ⍉ for other sizes +// COULD implement fast β‰βŸn +// SHOULD convert ⍉ with rank to a Reorder Axes call + #include "../core.h" #include "../utils/each.h" #include "../utils/talloc.h" @@ -154,7 +177,9 @@ B transp_im(B t, B x) { return def_fn_im(bi_transp, x); } -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))); +} void transp_init(void) { c(BFn,bi_transp)->uc1 = transp_uc1; From b555e3c0350aa1aee7eebfc2c3b3ccdbc04426a6 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Thu, 23 Mar 2023 20:56:25 -0400 Subject: [PATCH 02/26] Reorder Axes translated from runtime version --- src/builtins/transpose.c | 115 ++++++++++++++++++++++++++++++++++++++- src/load.c | 6 +- 2 files changed, 114 insertions(+), 7 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 3cadb5be..f6f8ce50 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -13,7 +13,7 @@ // COULD use half-width or smaller kernels to improve odd sizes // Scalar transpose or loop used for overhang of 1 -// Reorder Axes: self-hosted runtime (based on +⌜ and ⊏, not that slow) +// Reorder Axes: generate indices and select with +⌜ and ⊏ // Transpose inverse ⍉⁼ // Same as ⍉ for a rank ≀2 argument @@ -25,6 +25,7 @@ #include "../utils/each.h" #include "../utils/talloc.h" #include "../builtins.h" +#include "../utils/calls.h" #ifdef __BMI2__ #include @@ -51,7 +52,6 @@ #endif -extern B rt_transp; B transp_c1(B t, B x) { if (RARE(isAtm(x))) return m_atomUnit(x); ur xr = RNK(x); @@ -169,7 +169,116 @@ B transp_c1(B t, B x) { } decG(x); return taga(toBit? (Arr*)cpyBitArr(taga(r)) : r); } -B transp_c2(B t, B w, B x) { return c2rt(transp, w, x); } + +B mul_c2(B,B,B); +B ud_c1(B,B); +B tbl_c2(Md1D*,B,B); +B select_c2(B,B,B); + +B transp_c2(B t, B w, B x) { + usz wia=1; + if (isArr(w)) { + if (RNK(w)>1) 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); + 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); + } + + // compute shape for the given axes + usz* xsh = SH(x); + TALLOC(usz, rsh, xr); + usz dup = 0, max = 0; + usz no_sh = -(usz)1; + for (usz j=0; jmax? j : max; + if (xl= rr) thrF("⍉: Skipped result axis"); + if (wiaa, rsh, rr); + arr_shSetU(ra, rr, sh); + } + decG(x); + r = taga(ra); goto ret; + } + + // Number of axes that move + ur ar = max+1+dup; + if (ar == 1) { r = x; goto ret; } + // Add up stride for each axis + TALLOC(u64, st, rr); + for (usz j=0; j1) { + zsh = m_shArr(zr); + zsh->a[0] = c; + shcpy(zsh->a+1, xsh+ar, xr-ar); + } + Arr* z = TI(x,slice)(x, 0, IA(x)); + if (zr>1) arr_shSetU(z, zr, zsh); + else arr_shVec(z); + x = taga(z); + } + // (+⌜´stΓ—βŸœβ†•Β¨rsh)⊏β₯Šπ•© + B ind = bi_N; + for (ur k=ar-dup; k--; ) { + B v = C2(mul, m_f64(st[k]), C1(ud, m_f64(rsh[k]))); + if (q_N(ind)) ind = v; + else ind = M1C2(tbl, add, v, ind); + } + TFREE(st); + r = C2(select, ind, x); + + ret:; + TFREE(rsh); + TFREE(p); + return r; +} + B transp_im(B t, B x) { if (isAtm(x)) thrM("⍉⁼: 𝕩 must not be an atom"); 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); From dea063cb7b271228f81e6789df577177e644d566 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Thu, 23 Mar 2023 22:02:00 -0400 Subject: [PATCH 03/26] =?UTF-8?q?Easy=20optimizations:=20(1-=CB=9C=3D?= =?UTF-8?q?=F0=9D=95=A9)=E2=8D=89=F0=9D=95=A9=20to=20=E2=8D=89=F0=9D=95=A9?= =?UTF-8?q?=20and=20detect=20trailing=20unchanged=20axes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins/transpose.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index f6f8ce50..e1e04e81 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -190,6 +190,7 @@ B transp_c2(B t, B w, B x) { if (isAtm(w)) { usz a=o2s(w); if (a>=xr) thrF("⍉: Axis %s does not exist (%i≑=𝕩)", a, xr); + if (a==xr-1) { TFREE(p); return C1(transp, x); } p[0] = a; } else { SGetU(w) @@ -242,7 +243,8 @@ B transp_c2(B t, B w, B x) { // Number of axes that move ur ar = max+1+dup; - if (ar == 1) { r = x; goto ret; } + if (!dup) while (ar>1 && p[ar-1]==ar-1) ar--; // Unmoved trailing + if (ar <= 1) { r = x; goto ret; } // Add up stride for each axis TALLOC(u64, st, rr); for (usz j=0; j Date: Fri, 24 Mar 2023 16:17:33 -0400 Subject: [PATCH 04/26] =?UTF-8?q?Separate=20transpose=20el=5FB=20case=20an?= =?UTF-8?q?d=20fix=20dropped=20fill=20for=20shape=202=E2=80=BFn?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins/transpose.c | 95 +++++++++++++++++++--------------------- 1 file changed, 44 insertions(+), 51 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index e1e04e81..c2613303 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -60,8 +60,7 @@ B transp_c1(B t, B x) { usz ia = IA(x); usz* xsh = SH(x); usz h = xsh[0]; - if (ia==0 || h==1) { - no_reorder:; + if (ia==0 || h==1 || h==ia /*w==1*/) { Arr* r = cpyWithShape(x); ShArr* sh = m_shArr(xr); shcpy(sh->a, xsh+1, xr-1); @@ -70,45 +69,57 @@ 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); + if (xe==el_B) { + B xf = getFillR(x); + B* xp = TO_BPTR(x); + + HArr_p p = m_harrUp(ia); + if (h==2) { B* x0 = xp; B* x1 = x0+w; - HArr_p rp = m_harrUp(ia); - for (usz i=0; ia; - for (usz i=0; ia; + for (usz i=0; i Date: Sat, 25 Mar 2023 10:18:19 -0400 Subject: [PATCH 05/26] Separate out transpose_move function and use for el_B --- src/builtins/transpose.c | 124 +++++++++++++++++++-------------------- 1 file changed, 61 insertions(+), 63 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index c2613303..255a2d3b 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -34,8 +34,8 @@ #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; - for (usz i=0; i Date: Sat, 25 Mar 2023 10:43:35 -0400 Subject: [PATCH 06/26] =?UTF-8?q?Separate=20out=20transpose=20data=20movem?= =?UTF-8?q?ent=20function=20and=20share=20with=20monadic=20=E2=8D=89?= =?UTF-8?q?=E2=81=BC?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins/transpose.c | 99 +++++++++++++++++++++++----------------- 1 file changed, 56 insertions(+), 43 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 255a2d3b..50e574d0 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -15,9 +15,7 @@ // Reorder Axes: generate indices and select with +⌜ and ⊏ -// Transpose inverse ⍉⁼ -// Same as ⍉ for a rank ≀2 argument -// SHOULD share data movement with ⍉ for other sizes +// Transpose inverse ⍉⁼𝕩: data movement of ⍉ with different shape logic // COULD implement fast β‰βŸn // SHOULD convert ⍉ with rank to a Reorder Axes call @@ -77,45 +75,22 @@ static void transpose_move(void* rv, void* xv, u8 xe, usz w, usz h) { } } } - -B transp_c1(B t, B x) { - if (RARE(isAtm(x))) return m_atomUnit(x); - ur xr = RNK(x); - if (xr<=1) return x; - - usz ia = IA(x); - usz* xsh = SH(x); - usz h = xsh[0]; - if (ia==0 || h==1 || h==ia /*w==1*/) { - Arr* r = cpyWithShape(x); - ShArr* sh = m_shArr(xr); - shcpy(sh->a, xsh+1, xr-1); - sh->a[xr-1] = h; - arr_shReplace(r, xr, sh); - return taga(r); - } - usz w = xsh[1] * shProd(xsh, 2, xr); - - Arr* r; +// Return an array with data from x transposed as though it's shape h,w +// Shape of result needs to be set afterwards! +static Arr* transpose_noshape(B* px, usz ia, usz w, usz h) { + B x = *px; u8 xe = TI(x,elType); + Arr* r; if (xe==el_B) { B xf = getFillR(x); B* xp = TO_BPTR(x); - HArr_p p = m_harrUp(ia); + HArr_p p = m_harrUv(ia); // Debug build complains with harrUp transpose_move(p.a, xp, el_f64, w, h); for (usz xi=0; xia, xsh+1, xr-1); + sh->a[xr-1] = h; + arr_shReplace(r, xr, sh); + return taga(r); } + usz w = xsh[1] * shProd(xsh, 2, xr); + + Arr* r = transpose_noshape(&x, ia, w, h); + + usz* rsh = arr_shAlloc(r, xr); + if (xr==2) rsh[0] = w; else shcpy(rsh, SH(x)+1, xr-1); + rsh[xr-1] = h; decG(x); return taga(r); } @@ -275,8 +268,28 @@ B transp_c2(B t, B w, B x) { 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) { From ebbddf9dcf6813dbf9389a5ef54ea2fef4220fd2 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Sat, 25 Mar 2023 15:18:50 -0400 Subject: [PATCH 07/26] Move height/width 2 transposes to Singeli --- src/builtins/transpose.c | 26 +++++--------------------- src/singeli/src/transpose.singeli | 7 ++++++- 2 files changed, 11 insertions(+), 22 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 50e574d0..0af9e620 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -52,27 +52,11 @@ static void transpose_move(void* rv, void* xv, u8 xe, usz w, usz h) { assert(xe!=el_bit); assert(xe!=el_B); - if (h==2) { - switch(xe) { default: UD; - case el_i8: case el_c8: { u8* x0=xv; u8* x1=x0+w; u8* rp=rv; for (usz i=0; i Date: Sun, 26 Mar 2023 20:33:24 -0400 Subject: [PATCH 08/26] Move w128i and similar from sse3 and avx.singeli to base --- src/singeli/src/avx.singeli | 9 --------- src/singeli/src/base.singeli | 37 +++++++++++++++++++++++++++--------- src/singeli/src/sse3.singeli | 9 --------- 3 files changed, 28 insertions(+), 27 deletions(-) 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 From 410d51501d154215eccb7553f61d6c095d7facd3 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Sat, 25 Mar 2023 15:27:25 -0400 Subject: [PATCH 09/26] Check for AVX2 in Singeli transpose so it can build on non-x86 architectures --- build/src/build.bqn | 2 +- src/builtins/transpose.c | 2 +- src/singeli/src/transpose.singeli | 33 +++++++++++++++++-------------- 3 files changed, 20 insertions(+), 17 deletions(-) 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/transpose.c b/src/builtins/transpose.c index 0af9e620..0bd23027 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -35,7 +35,7 @@ #define TRANSPOSE_LOOP( DST, SRC, W, H) PLAINLOOP for(usz y=0,xi=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 ;y=k and h>=k) { + transpose_with_kernel{T, k, kh, call_base, rp, xp, w, h} + } else { + if (h==2) @for (x0 in xp, x1 in xp+w over i to w) { store{rp, i*2, x0}; store{rp, i*2+1, x1} } + else if (w==2) @for (r0 in rp, r1 in rp+h 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}} From b82907627120b74a7eb3b8d43711b0173e192519 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Sat, 25 Mar 2023 15:48:53 -0400 Subject: [PATCH 10/26] Put CPU-size transpose functions in a table --- src/builtins/transpose.c | 25 +++++++++++++------------ src/singeli/src/transpose.singeli | 10 ++++++---- 2 files changed, 19 insertions(+), 16 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 0bd23027..4e1f9985 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -32,32 +32,33 @@ #endif #endif -#define TRANSPOSE_LOOP( DST, SRC, W, H) PLAINLOOP for(usz y=0,xi=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 ;y Date: Sun, 26 Mar 2023 22:28:53 -0400 Subject: [PATCH 11/26] Transpose large cells with mut_copyG --- src/builtins/transpose.c | 51 +++++++++++++++++++++++++++++++++------- 1 file changed, 42 insertions(+), 9 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 4e1f9985..35786e71 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -144,6 +144,16 @@ B ud_c1(B,B); B tbl_c2(Md1D*,B,B); B select_c2(B,B,B); +static void shSet_alloc(Arr* ra, ur rr, usz* rsh) { + if (RARE(rr <= 1)) { + arr_shVec(ra); + } else { + ShArr* sh=m_shArr(rr); + shcpy(sh->a, rsh, rr); + arr_shSetU(ra, rr, sh); + } +} + B transp_c2(B t, B w, B x) { usz wia=1; if (isArr(w)) { @@ -199,13 +209,7 @@ B transp_c2(B t, B w, B x) { // Empty result if (IA(x) == 0) { Arr* ra = m_fillarrpEmpty(getFillQ(x)); - if (RARE(rr <= 1)) { - arr_shVec(ra); - } else { - ShArr* sh=m_shArr(rr); - shcpy(sh->a, rsh, rr); - arr_shSetU(ra, rr, sh); - } + shSet_alloc(ra, rr, rsh); decG(x); r = taga(ra); goto ret; } @@ -214,12 +218,40 @@ B transp_c2(B t, B w, B x) { ur ar = max+1+dup; if (!dup) while (ar>1 && p[ar-1]==ar-1) ar--; // Unmoved trailing if (ar <= 1) { r = x; goto ret; } + ur na = ar - dup; // Add up stride for each axis TALLOC(u64, st, rr); for (usz j=0; j= (32*8) >> elWidthLogBits(xe)) { // cell >= 32 bytes + usz ria = csz * shProd(rsh, 0, na); + MAKE_MUT_INIT(rm, ria, xe); MUTG_INIT(rm); + for (usz i=0; i Date: Mon, 27 Mar 2023 20:44:31 -0400 Subject: [PATCH 12/26] Fast bottom-level loop for Reorder Axes loop structure --- src/builtins/transpose.c | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 35786e71..15584b90 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -231,18 +231,24 @@ B transp_c2(B t, B w, B x) { usz ria = csz * shProd(rsh, 0, na); MAKE_MUT_INIT(rm, ria, xe); MUTG_INIT(rm); for (usz i=0; i=0 + ur a = na - 1; + usz str = st[a]; + usz l = rsh[a]; + for (usz k=0; k Date: Mon, 27 Mar 2023 22:12:12 -0400 Subject: [PATCH 13/26] SIMD kernel for Reorder Axes transposing last two axes --- src/builtins/transpose.c | 80 +++++++++++++++++++++++++++------------- 1 file changed, 55 insertions(+), 25 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 15584b90..57cb8f45 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -32,6 +32,7 @@ #endif #endif +typedef void (*TranspFn)(void*,void*,u64,u64); #if SINGELI #define transposeFns simd_transpose #define DECL_BASE(T) \ @@ -50,7 +51,7 @@ } DECL_BASE(i8) DECL_BASE(i16) DECL_BASE(i32) DECL_BASE(i64) #undef DECL_BASE - static void (*transposeFns[])(void*,void*,u64,u64) = { + static TranspFn transposeFns[] = { transpose_i8, transpose_i16, transpose_i32, transpose_i64 }; #endif @@ -227,36 +228,65 @@ B transp_c2(B t, B w, B x) { u8 xe = TI(x,elType); usz csz = shProd(xsh, ar, xr); - if (csz >= (32*8) >> elWidthLogBits(xe)) { // cell >= 32 bytes - usz ria = csz * shProd(rsh, 0, na); - MAKE_MUT_INIT(rm, ria, xe); MUTG_INIT(rm); - for (usz i=0; i=0 - ur a = na - 1; - usz str = st[a]; - usz l = rsh[a]; - for (usz k=0; k=1 */ \ + ur a = a0; \ + usz str = st[a]; \ + usz l = rsh[a]; \ + for (usz k=0; k= (32*8) >> xlw) { // cell >= 32 bytes + usz ria = csz * shProd(rsh, 0, na); + AXIS_LOOP(MAKE_MUT_INIT(rm, ria, xe); MUTG_INIT(rm); , + na, csz, mut_copyG(rm, i, x, j, csz)); Arr* ra = mut_fp(rm); shSet_alloc(ra, rr, rsh); r = withFill(taga(ra), getFillQ(x)); decG(x); goto ret_decst; } + if ((csz & (csz-1))==0 && csz<=64>>xlw && csz<=8 // CPU-sized cells + && !dup && na>=2 && p[na-1]==ar-2 && p[na-2]==ar-1 // Last two axes transposed + && rsh[na-2]*rsh[na-1] >= (256*8) >> xlw // And large-ish + && xe!=el_B) { + TranspFn tran = transposeFns[CTZ(csz< Date: Tue, 28 Mar 2023 16:32:41 -0400 Subject: [PATCH 14/26] Combine transp_c2 temp allocations into one TALLOC --- src/builtins/transpose.c | 29 ++++++++++++----------------- 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 57cb8f45..0390b2ce 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -166,7 +166,7 @@ B transp_c2(B t, B w, B x) { if (isAtm(x) || (xr=RNK(x))=xr) thrF("⍉: Axis %s does not exist (%i≑=𝕩)", a, xr); @@ -184,7 +184,7 @@ B transp_c2(B t, B w, B x) { // compute shape for the given axes usz* xsh = SH(x); - TALLOC(usz, rsh, xr); + usz *rsh = (usz*)(p + xr); // Length xr usz dup = 0, max = 0; usz no_sh = -(usz)1; for (usz j=0; j=1 */ \ ur a = a0; \ usz str = st[a]; \ - usz l = rsh[a]; \ for (usz k=0; k= (32*8) >> xlw) { // cell >= 32 bytes usz ria = csz * shProd(rsh, 0, na); - AXIS_LOOP(MAKE_MUT_INIT(rm, ria, xe); MUTG_INIT(rm); , - na, csz, mut_copyG(rm, i, x, j, csz)); + 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_alloc(ra, rr, rsh); r = withFill(taga(ra), getFillQ(x)); - decG(x); goto ret_decst; + decG(x); goto ret; } if ((csz & (csz-1))==0 && csz<=64>>xlw && csz<=8 // CPU-sized cells && !dup && na>=2 && p[na-1]==ar-2 && p[na-2]==ar-1 // Last two axes transposed @@ -280,11 +278,11 @@ B transp_c2(B t, B w, B x) { csz = (csz< Date: Tue, 28 Mar 2023 20:17:04 -0400 Subject: [PATCH 15/26] Save a shape object in Reorder Axes --- src/builtins/transpose.c | 60 +++++++++++++++++++--------------------- 1 file changed, 28 insertions(+), 32 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 0390b2ce..8bb00681 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -145,14 +145,9 @@ B ud_c1(B,B); B tbl_c2(Md1D*,B,B); B select_c2(B,B,B); -static void shSet_alloc(Arr* ra, ur rr, usz* rsh) { - if (RARE(rr <= 1)) { - arr_shVec(ra); - } else { - ShArr* sh=m_shArr(rr); - shcpy(sh->a, rsh, rr); - arr_shSetU(ra, rr, sh); - } +static void shSet(Arr* ra, ur rr, ShArr* sh) { + if (RARE(rr <= 1)) arr_shVec(ra); + else arr_shSetU(ra, rr, sh); } B transp_c2(B t, B w, B x) { @@ -182,7 +177,7 @@ B transp_c2(B t, B w, B x) { decG(w); } - // compute shape for the given axes + // Compute shape for the given axes usz* xsh = SH(x); usz *rsh = (usz*)(p + xr); // Length xr usz dup = 0, max = 0; @@ -196,7 +191,7 @@ B transp_c2(B t, B w, B x) { if (xl= rr) thrF("⍉: Skipped result axis"); if (wia 1)) { // Not all duplicates + sh = m_shArr(rr); + shcpy(sh->a, rsh, rr); + } B r; // Empty result if (IA(x) == 0) { Arr* ra = m_fillarrpEmpty(getFillQ(x)); - shSet_alloc(ra, rr, rsh); + shSet(ra, rr, sh); decG(x); r = taga(ra); goto ret; } @@ -218,11 +219,11 @@ B transp_c2(B t, B w, B x) { // Number of axes that move ur ar = max+1+dup; if (!dup) while (ar>1 && p[ar-1]==ar-1) ar--; // Unmoved trailing - if (ar <= 1) { r = x; goto ret; } + if (ar <= 1) { if (rr>1) ptr_dec(sh); r = x; goto ret; } ur na = ar - dup; // Add up stride for each axis - usz* st = rsh + xr; // Length rr - for (usz j=0; j=1 */ \ @@ -257,12 +258,12 @@ B transp_c2(B t, B w, B x) { 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_alloc(ra, rr, rsh); + shSet(ra, rr, sh); r = withFill(taga(ra), getFillQ(x)); decG(x); goto ret; } if ((csz & (csz-1))==0 && csz<=64>>xlw && csz<=8 // CPU-sized cells - && !dup && na>=2 && p[na-1]==ar-2 && p[na-2]==ar-1 // Last two axes transposed + && na>=2 && st[na-1]==rsh[na-2] && st[na-2]==1 // Last two axes transposed && rsh[na-2]*rsh[na-1] >= (256*8) >> xlw // And large-ish && xe!=el_B) { TranspFn tran = transposeFns[CTZ(csz<1) { - zsh = m_shArr(zr); - zsh->a[0] = c; - shcpy(zsh->a+1, xsh+ar, xr-ar); - } - Arr* z = TI(x,slice)(x, 0, IA(x)); - if (zr>1) arr_shSetU(z, zr, zsh); - else arr_shVec(z); - x = taga(z); - } + // Reshape x for selection + ShArr* zsh = m_shArr(2); + zsh->a[0] = c; + 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--; ) { @@ -308,6 +302,8 @@ B transp_c2(B t, B w, B x) { else ind = M1C2(tbl, add, v, ind); } r = C2(select, ind, x); + if (rr>1) arr_shReplace(a(r), rr, sh); + else { decSh(v(r)); arr_shVec(a(r)); } ret:; TFREE(p); From 814a6776764cfaab2a58e8c347fdb7043e440b88 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Tue, 28 Mar 2023 21:52:13 -0400 Subject: [PATCH 16/26] Make transpose movement functions strided --- src/builtins/transpose.c | 29 ++++++++----------- src/singeli/src/transpose.singeli | 46 +++++++++++++++---------------- 2 files changed, 35 insertions(+), 40 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 8bb00681..900b08a7 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -32,25 +32,20 @@ #endif #endif -typedef void (*TranspFn)(void*,void*,u64,u64); +#define DECL_BASE(T) \ + static NOINLINE void transpose_##T(void* rv, void* xv, u64 bw, u64 bh, u64 w, u64 h) { \ + T* rp=rv; T* xp=xv; \ + PLAINLOOP for(usz y=0;y 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 @@ -88,7 +88,7 @@ def transpose_with_kernel{T, k, kh, call_base, rp:*T, xp:*T, w, h} = { # 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 @@ -98,12 +98,12 @@ def transpose_with_kernel{T, k, kh, call_base, rp:*T, xp:*T, w, h} = { @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 @@ -118,57 +118,57 @@ def transpose_with_kernel{T, k, kh, call_base, rp:*T, xp:*T, w, h} = { each{tup, ...each{vt, iota{line_vecs}}} } ro := tail{6, -u64~~rp} / (width{T}/8) # Offset to align within cache line; assume elt-aligned - wh := w*h + 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) : void = { +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{'base_transpose_',ts}, ...a, w, h} + 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} + transpose_with_kernel{T, k, kh, call_base, rp, xp, w, h, ws, hs} } else { - if (h==2) @for (x0 in xp, x1 in xp+w over i to w) { store{rp, i*2, x0}; store{rp, i*2+1, x1} } - else if (w==2) @for (r0 in rp, r1 in rp+h over i to h) { r0 = load{xp, i*2}; r1 = load{xp, i*2+1} } + 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} } } From e923a718811d6b3ec0e6bbbd76a7605e0726e91d Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Wed, 29 Mar 2023 08:22:15 -0400 Subject: [PATCH 17/26] Use strided width in transpose-based Reorder Axes --- src/builtins/transpose.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 900b08a7..da9eaba2 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -258,23 +258,23 @@ B transp_c2(B t, B w, B x) { decG(x); goto ret; } if ((csz & (csz-1))==0 && csz<=64>>xlw && csz<=8 // CPU-sized cells - && na>=2 && st[na-1]==rsh[na-2] && st[na-2]==1 // Last two axes transposed - && rsh[na-2]*rsh[na-1] >= (256*8) >> xlw // And large-ish + && na>=2 && st[na-2]==1 // Last axis ends up second-to-last + && rsh[na-2]*rsh[na-1] >= (256*8) >> xlw // Large-ish axes && xe!=el_B) { TranspFn tran = transposeFns[CTZ(csz< Date: Wed, 29 Mar 2023 11:31:52 -0400 Subject: [PATCH 18/26] And strided height, so all orders are fast if axes are long --- src/builtins/transpose.c | 40 ++++++++++++++++++++++++++++++---------- 1 file changed, 30 insertions(+), 10 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index da9eaba2..5c75aadb 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -257,30 +257,50 @@ B transp_c2(B t, B w, B x) { 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 - && na>=2 && st[na-2]==1 // Last axis ends up second-to-last - && rsh[na-2]*rsh[na-1] >= (256*8) >> xlw // Large-ish axes - && xe!=el_B) { + && 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< Date: Wed, 29 Mar 2023 20:08:36 -0400 Subject: [PATCH 19/26] Full simplification pass for Reorder Axes --- src/builtins/transpose.c | 51 +++++++++++++++++++++++++--------------- 1 file changed, 32 insertions(+), 19 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 5c75aadb..c3f3d533 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -172,19 +172,23 @@ B transp_c2(B t, B w, B x) { 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; + usz dup = 0, max = 0, id = 0; usz no_sh = -(usz)1; for (usz j=0; jmax? j : max; if (xla, rsh, rr); } - B r; - // Empty result if (IA(x) == 0) { Arr* ra = m_fillarrpEmpty(getFillQ(x)); @@ -211,23 +213,33 @@ B transp_c2(B t, B w, B x) { r = taga(ra); goto ret; } - // Number of axes that move - ur ar = max+1+dup; - if (!dup) while (ar>1 && p[ar-1]==ar-1) ar--; // Unmoved trailing - if (ar <= 1) { if (rr>1) ptr_dec(sh); r = x; goto ret; } - ur na = ar - dup; // Add up stride for each axis - usz* st = rsh + xr; // Length ar - for (usz j=0; j=1 */ \ @@ -280,8 +292,9 @@ B transp_c2(B t, B w, B x) { usz i_skip = (w-1)*hs*csz; usz end = rf*csz - i_skip; ur a0 = na - 1; - for (usz i=0; i>= 3-xlw; + else if (xlw>3) for (usz i=0; ia[0] = c; + zsh->a[0] = IA(x)/csz; zsh->a[1] = csz; Arr* z = TI(x,slice)(x, 0, IA(x)); arr_shSetU(z, 2, zsh); @@ -312,7 +325,7 @@ B transp_c2(B t, B w, B x) { // (+⌜´stΓ—βŸœβ†•Β¨rsh)⊏β₯Šπ•© B ind = bi_N; for (ur k=na; k--; ) { - B v = C2(mul, m_f64(st[k]), C1(ud, m_f64(rsh[k]))); + B v = C2(mul, m_f64(st[k]/csz), C1(ud, m_f64(rsh[k]))); if (q_N(ind)) ind = v; else ind = M1C2(tbl, add, v, ind); } From 89016fc88dbe7b16b2d81e4eda316486753850af Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Wed, 29 Mar 2023 20:23:16 -0400 Subject: [PATCH 20/26] Update Reorder Axes implementation comments --- src/builtins/transpose.c | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index c3f3d533..fc4d3c5c 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -13,7 +13,17 @@ // COULD use half-width or smaller kernels to improve odd sizes // Scalar transpose or loop used for overhang of 1 -// Reorder Axes: generate indices and select with +⌜ and ⊏ +// Reorder Axes +// If 𝕨 indicates the identity permutation, return 𝕩 +// Simplify: remove length-1 axes; coalesce adjacent and trailing axes +// Empty result or trivial reordering: reshape 𝕩 +// Large cells: slow outer loop plus mut_copy +// CPU-sized cells, large last 𝕩 and result axes: strided 2D transposes +// Otherwise, generate indices and select with +⌜ and ⊏ +// SHOULD generate for a cell and virtualize the rest to save space +// COULD decompose axis permutations to use 2D transpose when possible +// COULD convert boolean to integer for some axis reorderings +// SHOULD have a small-subarray transposer using one or a few shuffles // Transpose inverse ⍉⁼𝕩: data movement of ⍉ with different shape logic // COULD implement fast β‰βŸn From d53f3dbd7208015ea99b2341133f8092c66b4e9a Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Thu, 30 Mar 2023 18:07:31 -0400 Subject: [PATCH 21/26] =?UTF-8?q?Native=20dyadic=20=E2=8D=89=E2=81=BC?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins/transpose.c | 42 +++++++++++++++++++++++++++++++++++++--- 1 file changed, 39 insertions(+), 3 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index fc4d3c5c..6855e6a2 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -379,7 +379,43 @@ B transp_uc1(B t, B o, B x) { return transp_im(m_f64(0), c1(o, transp_c1(t, x))); } -void transp_init(void) { - c(BFn,bi_transp)->uc1 = transp_uc1; - c(BFn,bi_transp)->im = transp_im; +B transp_ix(B t, B w, B x) { + if (isAtm(x)) thrM("⍉⁼: 𝕩 must not be an atom"); + ur xr=RNK(x); + + 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) { decG(w); return x; } + if (xr=xr) thrF("⍉⁼: Axis %s does not exist (%i≑=𝕩)", a, xr); + if (p[i]!=xr) thrM("⍉⁼: Duplicate axes"); + 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; iim = transp_im; + c(BFn,bi_transp)->ix = transp_ix; + c(BFn,bi_transp)->uc1 = transp_uc1; } From 9fe6a2e7b223b0c57aa00223ad75bfd0b817054b Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Thu, 30 Mar 2023 19:03:24 -0400 Subject: [PATCH 22/26] Under for invertible Reorder Axes --- src/builtins/transpose.c | 29 ++++++++++++++++++++++------- 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 6855e6a2..2a724507 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -25,7 +25,10 @@ // COULD convert boolean to integer for some axis reorderings // SHOULD have a small-subarray transposer using one or a few shuffles -// Transpose inverse ⍉⁼𝕩: data movement of ⍉ with different shape logic +// ⍉⁼𝕩: data movement of ⍉ with different shape logic +// 𝕨⍉⁼𝕩: compute inverse 𝕨, length 1+βŒˆΒ΄π•¨ +// Under Transpose supports invertible cases +// SHOULD implement Under with duplicate axes, maybe as Under Select // COULD implement fast β‰βŸn // SHOULD convert ⍉ with rank to a Reorder Axes call @@ -379,10 +382,8 @@ B transp_uc1(B t, B o, B x) { return transp_im(m_f64(0), c1(o, transp_c1(t, x))); } -B transp_ix(B t, B w, B x) { - if (isAtm(x)) thrM("⍉⁼: 𝕩 must not be an atom"); - ur xr=RNK(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); @@ -392,7 +393,7 @@ B transp_ix(B t, B w, B x) { } else { if (RNK(w)>1) thrM("⍉⁼: 𝕨 must have rank at most 1"); usz wia = IA(w); - if (wia==0) { decG(w); return x; } + if (wia==0) return w; if (xr=xr) thrF("⍉⁼: Axis %s does not exist (%i≑=𝕩)", a, xr); - if (p[i]!=xr) thrM("⍉⁼: Duplicate axes"); + if (p[a]!=xr) { TFREE(p); decG(w); return bi_N; } // Handled by caller max = a>max? a : max; p[a] = i; } @@ -411,11 +412,25 @@ B transp_ix(B t, B w, B x) { for (usz i=0, j=wia; iim = transp_im; c(BFn,bi_transp)->ix = transp_ix; c(BFn,bi_transp)->uc1 = transp_uc1; + c(BFn,bi_transp)->ucw = transp_ucw; } From 69e2f6c69044e13e8aaad9aaefef601f6c570978 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Thu, 30 Mar 2023 21:10:29 -0400 Subject: [PATCH 23/26] =?UTF-8?q?Convert=20=E2=8D=89=CB=98=20with=20atom?= =?UTF-8?q?=20or=20no=20left=20argument=20to=20dyadic=20=E2=8D=89?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins/md1.c | 8 ++++++++ src/builtins/transpose.c | 3 ++- 2 files changed, 10 insertions(+), 1 deletion(-) 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 Date: Fri, 31 Mar 2023 13:44:32 +0300 Subject: [PATCH 24/26] don't assume sizeof(ur)==1 --- src/builtins/transpose.c | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 92f7f5a7..f71fcb10 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -170,11 +170,12 @@ B transp_c2(B t, B w, B x) { if (isAtm(x) || (xr=RNK(x))=xr) thrF("⍉: Axis %s does not exist (%i≑=𝕩)", a, xr); - if (a==xr-1) { TFREE(p); return C1(transp, x); } + if (a==xr-1) { TFREE(alloc); return C1(transp, x); } p[0] = a; } else { SGetU(w) @@ -190,7 +191,7 @@ B transp_c2(B t, B w, B x) { // Compute shape for the given axes usz* xsh = SH(x); - usz *rsh = (usz*)(p + xr); // Length xr + usz* rsh = (usz*)(p + xr); // Length xr usz dup = 0, max = 0, id = 0; usz no_sh = -(usz)1; for (usz j=0; j Date: Fri, 31 Mar 2023 14:14:01 +0300 Subject: [PATCH 25/26] don't assume select returns refcount 1 array that probably would be a safe assumption right now, but who knows about the future; cpyWithShape should return back its argument directly if it already has refcount 1 anyway --- src/builtins/transpose.c | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index f71fcb10..3ef94226 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -345,8 +345,9 @@ B transp_c2(B t, B w, B x) { else ind = M1C2(tbl, add, v, ind); } r = C2(select, ind, x); - if (rr>1) arr_shReplace(a(r), rr, sh); - else { decSh(v(r)); arr_shVec(a(r)); } + 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); From c513724fd4de57868bb2cda3afa6655058bf6f3d Mon Sep 17 00:00:00 2001 From: dzaima Date: Fri, 31 Mar 2023 14:26:25 +0300 Subject: [PATCH 26/26] hide clangd warning and remove whitespace which was previously used for alignment --- src/builtins/transpose.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/builtins/transpose.c b/src/builtins/transpose.c index 3ef94226..0cc0d4be 100644 --- a/src/builtins/transpose.c +++ b/src/builtins/transpose.c @@ -340,7 +340,7 @@ B transp_c2(B t, B w, B x) { // (+⌜´stΓ—βŸœβ†•Β¨rsh)⊏β₯Šπ•© B ind = bi_N; for (ur k=na; k--; ) { - B v = C2(mul, m_f64(st[k]/csz), C1(ud, m_f64(rsh[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); } @@ -382,7 +382,7 @@ B transp_im(B t, B x) { } B transp_uc1(B t, B o, B x) { - return transp_im(m_f64(0), c1(o, transp_c1(t, x))); + return transp_im(m_f64(0), c1(o, transp_c1(t, x))); } // Consumes w; return bi_N if w contained duplicates