From 04fff34f4680c17fceeb8490403ef1471d2925b8 Mon Sep 17 00:00:00 2001 From: Benjamin Schwendinger Date: Mon, 26 Oct 2026 22:54:06 +0100 Subject: [PATCH 1/2] add per thread buffers and per-group parallization to gmedian and gsd/gvar --- inst/tests/optimize.Rraw | 7 +++++++ src/gsumm.c | 30 +++++++++++++++++++++--------- 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/inst/tests/optimize.Rraw b/inst/tests/optimize.Rraw index 9d246cd4fc..aac03a41bb 100644 --- a/inst/tests/optimize.Rraw +++ b/inst/tests/optimize.Rraw @@ -508,8 +508,15 @@ if (test_bit64) { test(2284.6, optimize=c(0L, Inf), DT[, .(sum(i64), sum(i64, na.rm=TRUE)), by=g1]) test(2284.7, optimize=c(0L, Inf), DT[, .(sum(i64), sum(i64, na.rm=TRUE)), by=g2]) } +setDTthreads(old) # GForce skips the gforce() setup when j is only .N DT = data.table(a=c(2L,1L,2L,3L,2L)) test(2285.1, optimize=2L, DT[, .N, by=a, verbose=TRUE], output="GForce optimized j to '.N'", notOutput="gforce initial population") test(2285.2, optimize=0:2, DT[, list(".N"), by=a], data.table(a=c(2L,1L,3L), V1=".N")) + +# gmedian/gvar/gsd run in parallel over groups using per-thread buffers of gx +set.seed(3L) +DT = data.table(g=sample(c(rep(0L, 600L), rep(1:200, 2L))), i=sample(c(1:3, NA), 1000L, TRUE), d=sample(c(rnorm(3), NA), 1000L, TRUE)) +test(2286.1, optimize=c(0L, Inf), DT[, .(var(i), var(i, na.rm=TRUE), sd(i), var(d), sd(d, na.rm=TRUE), median(d), median(d, na.rm=TRUE)), by=g]) # group 0 > nrow/2 gives each group its own buffer with >=2 threads +test(2286.2, optimize=c(0L, Inf), DT[g > 0L, .(var(i), var(i, na.rm=TRUE), sd(i), var(d), sd(d, na.rm=TRUE), median(d), median(d, na.rm=TRUE)), by=g]) # per-thread buffers diff --git a/src/gsumm.c b/src/gsumm.c index d2adb81654..33f56cc4ce 100644 --- a/src/gsumm.c +++ b/src/gsumm.c @@ -930,6 +930,12 @@ SEXP gmax(SEXP x, SEXP narm) return gminmax(x, narm, false); } +// offset of group i's scratch space in gx for gmedian/gvarsd1: each thread reuses its own cache-hot buffer of maxgrpn, +// unless nth*maxgrpn exceeds nrow due to a huge group; then each group gets its own buffer +static inline size_t my_buffer(int i, bool pergroup) { + return pergroup ? ff[i]-1 : (size_t)omp_get_thread_num()*maxgrpn; +} + // gmedian, always returns numeric type (to avoid as.numeric() wrap..) SEXP gmedian(SEXP x, SEXP narmArg) { if (!IS_TRUE_OR_FALSE(narmArg)) @@ -943,12 +949,15 @@ SEXP gmedian(SEXP x, SEXP narmArg) { SEXP ans = PROTECT(allocVector(REALSXP, ngrp)); double *ansd = REAL(ans); const bool nosubset = irowslen==-1; + const int nth = getDTthreads(ngrp, false); + const bool pergroup = (int64_t)nth*maxgrpn > nrow; switch(TYPEOF(x)) { case REALSXP: { - double *subd = REAL(PROTECT(allocVector(REALSXP, maxgrpn))); // allocate once upfront and reuse const int64_t *xi64 = (const int64_t *)REAL_RO(x); const double *xd = REAL_RO(x); + #pragma omp parallel for schedule(dynamic) num_threads(nth) for (int i=0; i nrow; switch(TYPEOF(x)) { case LGLSXP: case INTSXP: { - sub = PROTECT(allocVector(INTSXP, maxgrpn)); // allocate once upfront - int *subd = INTEGER(sub); const int *xd = INTEGER_RO(x); + #pragma omp parallel for schedule(dynamic) num_threads(nth) for (int i=0; i Date: Wed, 28 Oct 2026 13:49:45 +0100 Subject: [PATCH 2/2] use thread gate --- src/gsumm.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gsumm.c b/src/gsumm.c index 33f56cc4ce..23b9d23733 100644 --- a/src/gsumm.c +++ b/src/gsumm.c @@ -949,7 +949,7 @@ SEXP gmedian(SEXP x, SEXP narmArg) { SEXP ans = PROTECT(allocVector(REALSXP, ngrp)); double *ansd = REAL(ans); const bool nosubset = irowslen==-1; - const int nth = getDTthreads(ngrp, false); + const int nth = getDTthreads(nrow, true); const bool pergroup = (int64_t)nth*maxgrpn > nrow; switch(TYPEOF(x)) { case REALSXP: { @@ -1108,7 +1108,7 @@ static SEXP gvarsd1(SEXP x, SEXP narmArg, bool isSD) double *ansd = REAL(ans); const bool nosubset = irowslen==-1; const bool narm = LOGICAL(narmArg)[0]; - const int nth = getDTthreads(ngrp, false); + const int nth = getDTthreads(nrow, true); const bool pergroup = (int64_t)nth*maxgrpn > nrow; switch(TYPEOF(x)) { case LGLSXP: case INTSXP: {