Use Singeli max/min scan for Indices and Sort

This commit is contained in:
Marshall Lochbaum 2022-09-21 17:29:10 -04:00
parent bc4f5167fb
commit 8e1d8bb42c
2 changed files with 62 additions and 17 deletions

View File

@ -30,21 +30,47 @@
rp[j] = xi; \
}
#if SINGELI
extern void (*const avx2_scan_max8)(int8_t* v0,int8_t* v1,uint64_t v2);
extern void (*const avx2_scan_min8)(int8_t* v0,int8_t* v1,uint64_t v2);
extern void (*const avx2_scan_max16)(int16_t* v0,int16_t* v1,uint64_t v2);
extern void (*const avx2_scan_min16)(int16_t* v0,int16_t* v1,uint64_t v2);
#define COUNT_THRESHOLD 32
#define WRITE_SPARSE_i8 \
for (usz i=0; i<n; i++) rp[i]=j; \
while (ij<n) { rp[ij]=GRADE_UD(++j,--j); ij+=c0o[j]; } \
GRADE_UD(avx2_scan_max8,avx2_scan_min8)(rp,rp,n);
#define WRITE_SPARSE_i16 \
usz b = 1<<10; \
for (usz k=0; ; ) { \
usz e = b<n-k? k+b : n; \
for (usz i=k; i<e; i++) rp[i]=j; \
while (ij<e) { rp[ij]=GRADE_UD(++j,--j); ij+=c0o[j]; } \
GRADE_UD(avx2_scan_max16,avx2_scan_min16)(rp+k,rp+k,e-k); \
if (e==n) {break;} k=e; \
}
#define WRITE_SPARSE(T) WRITE_SPARSE_##T
#else
#define COUNT_THRESHOLD 16
#define WRITE_SPARSE(T) \
for (usz i=0; i<n; i++) rp[i]=0; \
usz js = j; \
while (ij<n) { rp[ij]GRADE_UD(++,--); ij+=c0o[GRADE_UD(++j,--j)]; } \
for (usz i=0; i<n; i++) js=rp[i]+=js;
#endif
#define COUNTING_SORT(T) \
usz C=1<<(8*sizeof(T)); \
TALLOC(usz, c0, C); usz *c0o=c0+C/2; \
for (usz j=0; j<C; j++) c0[j]=0; \
for (usz i=0; i<n; i++) c0o[xp[i]]++; \
if (n/16 <= C) { /* Sum-based */ \
for (usz i=0; i<n; i++) rp[i]=0; \
usz j=GRADE_UD(0,C-1), i; \
while ((i=c0[j])==0) GRADE_UD(j++,j--); \
usz js = j - C/2; \
while (i<n) { rp[i]++; i+=c0[GRADE_UD(++j,--j)]; } \
for (usz i=0; i<n; i++) js=rp[i]+=js; \
} else { /* Branchy */ \
FOR(j,C) for (usz c=c0[j]; c--; ) *rp++ = j-C/2; \
} \
usz C=1<<(8*sizeof(T)); \
TALLOC(usz, c0, C); usz *c0o=c0+C/2; \
for (usz j=0; j<C; j++) c0[j]=0; \
for (usz i=0; i<n; i++) c0o[xp[i]]++; \
if (n/(COUNT_THRESHOLD*sizeof(T)) <= C) { /* Scan-based */ \
T j=GRADE_UD(-C/2,C/2-1); \
usz ij; while ((ij=c0o[j])==0) GRADE_UD(j++,j--); \
WRITE_SPARSE(T) \
} else { /* Branchy */ \
FOR(j,C) for (usz c=c0[j]; c--; ) *rp++ = j-C/2; \
} \
TFREE(c0)
// Radix sorting
@ -181,6 +207,10 @@ B SORT_C1(B t, B x) {
#undef SORT_C1
#undef INSERTION_SORT
#undef COUNTING_SORT
#if SINGELI
#undef WRITE_SPARSE_i8
#undef WRITE_SPARSE_i16
#endif
#define GRADE_CHR GRADE_UD("⍋","⍒")

View File

@ -110,6 +110,7 @@ extern void (*const avx2_scan_pluswrap_u16)(uint16_t* v0,uint16_t* v1,uint64_t v
extern void (*const avx2_scan_pluswrap_u32)(uint32_t* v0,uint32_t* v1,uint64_t v2,uint32_t v3);
#define avx2_scan_pluswrap_u64(V0,V1,V2,V3) for (usz i=k; i<e; i++) js=rp[i]+=js;
#define PLUS_SCAN(T) avx2_scan_pluswrap_##T(rp+k,rp+k,e-k,js); js=rp[e-1];
extern void (*const avx2_scan_max32)(int32_t* v0,int32_t* v1,uint64_t v2);
#else
#define PLUS_SCAN(T) for (usz i=k; i<e; i++) js=rp[i]+=js;
#endif
@ -475,7 +476,20 @@ B slash_c1(B t, B x) {
for (u64 j = 0; j < c; j++) *rp++ = i;
}
} else {
if (s/16 <= xia) { // Sparse case: type of x matters
if (s/32 <= xia) { // Sparse case: type of x matters
#if SINGELI
#define SPARSE_IND(T) \
T* xp = T##any_ptr(x); \
usz b = 1<<10; \
for (usz k=0, j=0, ij=xp[0]; ; ) { \
usz e = b<s-k? k+b : s; \
for (usz i=k; i<e; i++) rp[i]=0; \
rp[k]=j; \
while (ij<e) { rp[ij]=++j; ij+=xp[j]; } \
avx2_scan_max32(rp+k,rp+k,e-k); \
if (e==s) {break;} k=e; \
}
#else
#define SPARSE_IND(T) \
T* xp = T##any_ptr(x); \
usz b = 1<<10; \
@ -483,10 +497,11 @@ B slash_c1(B t, B x) {
usz e = b<s-k? k+b : s; \
for (usz i=k; i<e; i++) rp[i]=0; \
while (ij<e) { rp[ij]++; ij+=xp[++j]; } \
PLUS_SCAN(u32) \
PLUS_SCAN(i32) \
if (e==s) {break;} k=e; \
}
u32* rp; r = m_i32arrv((i32**)&rp, s);
#endif
i32* rp; r = m_i32arrv(&rp, s);
if (xe == el_i8 ) { SPARSE_IND(i8 ); }
else if (xe == el_i16) { SPARSE_IND(i16); }
else { SPARSE_IND(i32); }