From 9bc4dd94f3fbe8efa1d9258944e097c9315fdaf9 Mon Sep 17 00:00:00 2001 From: Benjamin Schwendinger Date: Sun, 18 Oct 2026 22:00:01 +0200 Subject: [PATCH 1/4] speed up single-threaded GForce sum/mean by skipping batching (#3330) With one thread, use a single batch and put as much of the group id into low as fits in uint16_t. For up to 65536 groups there is then a single high bucket, gather() returns x itself instead of copying it, and gsum and gmean reduce to direct loops. --- NEWS.md | 2 ++ inst/tests/tests.Rraw | 19 +++++++++++++++++++ src/gsumm.c | 33 +++++++++++++++++++++++++++++++-- 3 files changed, 52 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index cb4a6231fb..f1ea56e2f8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -114,6 +114,8 @@ 10. `fwrite()` returns a clearer error message when `na = data.frame()` is used, [#7866](https://github.com/Rdatatable/data.table/issues/7866). Thanks @mcol for the report and the fix. +11. GForce grouping is faster when running single-threaded, e.g. after `setDTthreads(1)`, [#3330](https://github.com/Rdatatable/data.table/issues/3330). The batching used to parallelize `sum()` and `mean()` is skipped when only one thread is available, and with up to 65,536 groups those functions now read the column directly instead of first copying it. Thanks @jangorecki for the benchmarks, @mattdowle for the report, and @ben-schwen for the fix. + ## data.table [v1.18.6.1](https://github.com/Rdatatable/data.table/milestone/47) (22 August 2026) diff --git a/inst/tests/tests.Rraw b/inst/tests/tests.Rraw index 78123b73fe..1262120bc7 100644 --- a/inst/tests/tests.Rraw +++ b/inst/tests/tests.Rraw @@ -22011,3 +22011,22 @@ test(2388.21, nafill(c(1, NA, NA, 5), type="nocb", limit=1, nan=NaN), c(1, NA, 5 test(2388.22, nafill(c(1L, NA, NA, 5L), type="nocb", limit=1), c(1L, NA, 5L, 5L)) test(2388.23, nafill(c("a", NA, NA), type="const", fill="x", limit=1), c("a", "x", NA)) test(2388.24, nafill(c("a", NA, NA, "z"), type="nocb", limit=1), c("a", NA, "z", "z")) + +# single-threaded GForce uses one batch and, with up to 65536 groups, gsum/gmean read x directly instead of gathering it, #3330 +old = setDTthreads(1L) +set.seed(3330L) +N = 2e5L +DT = data.table(g1=sample(100L, N, TRUE), g2=sample(70000L, N, TRUE), i=sample(c(1:10, NA), N, TRUE), k=sample(5L, N, TRUE), d=rnorm(N), c=complex(real=rnorm(N), imaginary=rnorm(N))) +DT[c(3L, 7L), d := NA][c(5L, 9L), c := NA] +test(2389.1, optimize=c(0L, Inf), DT[, .(sum(i), sum(i, na.rm=TRUE), mean(i), mean(i, na.rm=TRUE), sum(k), mean(k), sum(d), sum(d, na.rm=TRUE), mean(d), mean(d, na.rm=TRUE), sum(c), sum(c, na.rm=TRUE), mean(c), mean(c, na.rm=TRUE)), keyby=g1]) +test(2389.2, optimize=c(0L, Inf), DT[, .(sum(i), sum(i, na.rm=TRUE), mean(i), mean(i, na.rm=TRUE), sum(k), mean(k), sum(d), sum(d, na.rm=TRUE), mean(d), mean(d, na.rm=TRUE), sum(c), sum(c, na.rm=TRUE), mean(c), mean(c, na.rm=TRUE)), by=g2]) +test(2389.3, optimize=c(0L, Inf), DT[d > 0, .(sum(i), mean(i, na.rm=TRUE), sum(k), mean(d), sum(c)), by=g1]) +setkey(DT, g1) +test(2389.4, optimize=c(0L, Inf), DT[, .(sum(i), sum(i, na.rm=TRUE), mean(i), mean(i, na.rm=TRUE), sum(k), mean(k), sum(d), sum(d, na.rm=TRUE), mean(d), mean(d, na.rm=TRUE), sum(c), sum(c, na.rm=TRUE), mean(c), mean(c, na.rm=TRUE)), by=g1]) +test(2389.5, options=list(datatable.optimize=Inf), data.table(g=c(1L, 1L, 2L), x=c(.Machine$integer.max, 1L, 1L))[, sum(x), by=g], data.table(g=1:2, V1=c(2147483648, 1)), warning="coerced to 'numeric'") +if (test_bit64) { + DT[, i64 := as.integer64(i)] + test(2389.6, optimize=c(0L, Inf), DT[, .(sum(i64), sum(i64, na.rm=TRUE)), by=g1]) + test(2389.7, optimize=c(0L, Inf), DT[, .(sum(i64), sum(i64, na.rm=TRUE)), by=g2]) +} +setDTthreads(old) diff --git a/src/gsumm.c b/src/gsumm.c index 5051d72ca8..42b5164c88 100644 --- a/src/gsumm.c +++ b/src/gsumm.c @@ -75,8 +75,11 @@ SEXP gforce(SEXP env, SEXP jsub, SEXP o, SEXP f, SEXP l, SEXP irowsArg) { if (length(tt)==1 && INTEGER(tt)[0]!=maxgrpn) internal_error(__func__, "o's maxgrpn attribute mismatches recalculated maxgrpn"); // # nocov } + // Single-threaded, the batching below is pure overhead (#3330). Use one batch and put as much of g into low as fits in + // uint16_t, so that for ngrp<=65536 there is one high bucket and gather/gsum/gmean reduce to direct loops. + const bool single = getDTthreads(nrow, false)==1; int nb = nbit(ngrp-1); - bitshift = nb/2; // /2 so that high and low can be uint16_t, and no limit (even for nb=4) to stress-test. + bitshift = single ? MIN(nb, 16) : nb/2; // /2 so that high and low can be uint16_t, and no limit (even for nb=4) to stress-test. // bitshift=MAX(nb-8,0); if (bitshift>16) bitshift=nb/2; // TODO: when we have stress-test off mode, do this mask = (1<>bitshift) + 1; @@ -85,7 +88,7 @@ SEXP gforce(SEXP env, SEXP jsub, SEXP o, SEXP f, SEXP l, SEXP irowsArg) { // maybe better to malloc to avoid R's heap. This grp isn't global, so it doesn't need to be R_alloc const int *restrict fp = INTEGER_RO(f); - nBatch = MIN((nrow+1)/2, getDTthreads(nrow, true)*2); // *2 to reduce last-thread-home. TODO: experiment. The higher this is though, the bigger is counts[] + nBatch = single ? 1 : MIN((nrow+1)/2, getDTthreads(nrow, true)*2); // *2 to reduce last-thread-home. TODO: experiment. The higher this is though, the bigger is counts[] batchSize = MAX(1, (nrow-1)/nBatch); lastBatchSize = nrow - (nBatch-1)*batchSize; // We deliberate use, for example, 40 batches of just 14 rows, to stress-test tests. This strategy proved to be a good one as #3204 immediately came to light. @@ -222,6 +225,32 @@ void *gather(SEXP x, bool *anyNA) { double started=wallclock(); const bool verbose = GetVerbose(); + if (nBatch==1 && highSize==1 && irowslen==-1) { + // gather would be an identical copy of x, so return x itself after checking for NA (#3330). Callers only read gx. + switch (TYPEOF(x)) { + case LGLSXP: case INTSXP: { + const int *xp = INTEGER_RO(x); + for (int i=0; i Date: Sun, 18 Oct 2026 22:00:01 +0200 Subject: [PATCH 2/4] speed up single-threaded GForce sum/mean by skipping batching (#3330) With one thread, use a single batch and put as much of the group id into low as fits in uint16_t. For up to 65536 groups there is then a single high bucket, gather() returns x itself instead of copying it, and gsum and gmean reduce to direct loops. --- NEWS.md | 4 +++- inst/tests/optimize.Rraw | 23 +++++++++++++++++++++-- src/gsumm.c | 33 +++++++++++++++++++++++++++++++-- 3 files changed, 55 insertions(+), 5 deletions(-) diff --git a/NEWS.md b/NEWS.md index cb4a6231fb..6045c590cc 100644 --- a/NEWS.md +++ b/NEWS.md @@ -92,6 +92,9 @@ 20. `setnafill()` with `type = "locf"` or `type = "nocb"` on a 0-row table no longer segfaults, [#7910](https://github.com/Rdatatable/data.table/issues/7910). Thanks @MichaelChirico for the report and fix. +21. GForce grouping is again faster when running single-threaded, e.g. after `setDTthreads(1)`, [#3330](https://github.com/Rdatatable/data.table/issues/3330). The batching used to parallelize `gsum()` and `gmean()` is skipped when only one thread is available. Thanks @jangorecki for the benchmarks, @mattdowle for the report, and @ben-schwen for the fix. + + ### Notes 1. {data.table} now depends on R 3.5.0 (2018). @@ -114,7 +117,6 @@ 10. `fwrite()` returns a clearer error message when `na = data.frame()` is used, [#7866](https://github.com/Rdatatable/data.table/issues/7866). Thanks @mcol for the report and the fix. - ## data.table [v1.18.6.1](https://github.com/Rdatatable/data.table/milestone/47) (22 August 2026) ### NEW FEATURES diff --git a/inst/tests/optimize.Rraw b/inst/tests/optimize.Rraw index 4c25baee8b..6afe56b65b 100644 --- a/inst/tests/optimize.Rraw +++ b/inst/tests/optimize.Rraw @@ -110,7 +110,7 @@ test(1184.2,optimize=1L, DT[, mean(v), by=x, verbose=TRUE], output="(GForce FALS test(1185.1,optimize=c(0L, 1L, 2L), DT[, list(sum(y), sum(v), sum(y,na.rm=TRUE), sum(v,na.rm=TRUE)), by=x], data.table(x=c("a","b","c","d"), V1=c(NA,10L,NA,NA), V2=c(6,NA,NA,NA), V3=c(4L,10L,7L,0L), V4=c(6,10,15,0))) # test 1185.1 subsumes 1185.1 and 1187.1 for testing all levels -test(1185.2,optimize=c(0L,1L,Inf), DT[, list(mean(y), mean(v), mean(y,na.rm=TRUE), mean(v,na.rm=TRUE)), by=x, verbose=TRUE], +test(1185.2,optimize=c(0L,1L,Inf), DT[, list(mean(y), mean(v), mean(y,na.rm=TRUE), mean(v,na.rm=TRUE)), by=x, verbose=TRUE], output=c("All optimizations.*off", "Old mean.*changed j", "GForce optimized j to")) # test 1185.2 subsumes 1185.2, 1185.3, and 1185.4 for testing all levels MyVar = TRUE @@ -334,7 +334,7 @@ test(1594.08,optimize=opt, DT[(idx), var(z1, na.rm=TRUE), by=x], ans) test(1594.09,optimize=opt,DT[, lapply(.SD, sd, na.rm=FALSE), by=x]) test(1594.10,optimize=opt, DT[, lapply(.SD, sd, na.rm=TRUE), by=x], DT[, lapply(.SD, stats::sd, na.rm=TRUE), by=x]) test(1594.11,optimize=opt, DT[, lapply(.SD, sd, na.rm=TRUE), by=x, verbose=TRUE], output=out) - + test(1594.12,optimize=opt, DT[, lapply(.SD, prod, na.rm=FALSE), by=x]) test(1594.13,optimize=opt, DT[, lapply(.SD, prod, na.rm=TRUE), by=x]) test(1594.14,optimize=opt, DT[, lapply(.SD, prod, na.rm=TRUE), by=x, verbose=TRUE], output=out) @@ -490,3 +490,22 @@ test(2283 + 0.19, optimize=opts, names(M[, c(mpg[1], list(mpg, b=hp), c(lapply(.SD, mean))), by="cyl", .SDcols=c("vs", "am")]), c("cyl", "V1", "V2", "b", "vs", "am"), context=sprintf("optimize=%s [XIX]", format(opt))) + +# single-threaded GForce uses one batch and gsum/gmean read x directly instead of gathering it, #3330 +old = setDTthreads(1L) +set.seed(3330L) +N = 2e5L +DT = data.table(g1=sample(100L, N, TRUE), g2=sample(70000L, N, TRUE), i=sample(c(1:10, NA), N, TRUE), k=sample(5L, N, TRUE), d=rnorm(N), c=complex(real=rnorm(N), imaginary=rnorm(N))) +DT[c(3L, 7L), d := NA][c(5L, 9L), c := NA] +test(2284.1, optimize=c(0L, Inf), DT[, .(sum(i), sum(i, na.rm=TRUE), mean(i), mean(i, na.rm=TRUE), sum(k), mean(k), sum(d), sum(d, na.rm=TRUE), mean(d), mean(d, na.rm=TRUE), sum(c), sum(c, na.rm=TRUE), mean(c), mean(c, na.rm=TRUE)), keyby=g1]) +test(2284.2, optimize=c(0L, Inf), DT[, .(sum(i), sum(i, na.rm=TRUE), mean(i), mean(i, na.rm=TRUE), sum(k), mean(k), sum(d), sum(d, na.rm=TRUE), mean(d), mean(d, na.rm=TRUE), sum(c), sum(c, na.rm=TRUE), mean(c), mean(c, na.rm=TRUE)), by=g2]) +test(2284.3, optimize=c(0L, Inf), DT[d > 0, .(sum(i), mean(i, na.rm=TRUE), sum(k), mean(d), sum(c)), by=g1]) +setkey(DT, g1) +test(2284.4, optimize=c(0L, Inf), DT[, .(sum(i), sum(i, na.rm=TRUE), mean(i), mean(i, na.rm=TRUE), sum(k), mean(k), sum(d), sum(d, na.rm=TRUE), mean(d), mean(d, na.rm=TRUE), sum(c), sum(c, na.rm=TRUE), mean(c), mean(c, na.rm=TRUE)), by=g1]) +test(2284.5, options=list(datatable.optimize=Inf), data.table(g=c(1L, 1L, 2L), x=c(.Machine$integer.max, 1L, 1L))[, sum(x), by=g], data.table(g=1:2, V1=c(2147483648, 1)), warning="coerced to 'numeric'") +if (test_bit64) { + DT[, i64 := as.integer64(i)] + 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) diff --git a/src/gsumm.c b/src/gsumm.c index 5051d72ca8..88d3cdb7de 100644 --- a/src/gsumm.c +++ b/src/gsumm.c @@ -75,8 +75,11 @@ SEXP gforce(SEXP env, SEXP jsub, SEXP o, SEXP f, SEXP l, SEXP irowsArg) { if (length(tt)==1 && INTEGER(tt)[0]!=maxgrpn) internal_error(__func__, "o's maxgrpn attribute mismatches recalculated maxgrpn"); // # nocov } + // Single-threaded, the batching below is pure overhead (#3330). Use one batch and put as much of g into low as fits in + // uint16_t, so that for ngrp<=65536 there is one high bucket and gather/gsum/gmean reduce to direct loops. + const bool single = getDTthreads(nrow, false)==1; int nb = nbit(ngrp-1); - bitshift = nb/2; // /2 so that high and low can be uint16_t, and no limit (even for nb=4) to stress-test. + bitshift = single ? MIN(nb, 16) : nb/2; // /2 so that high and low can be uint16_t, and no limit (even for nb=4) to stress-test. // bitshift=MAX(nb-8,0); if (bitshift>16) bitshift=nb/2; // TODO: when we have stress-test off mode, do this mask = (1<>bitshift) + 1; @@ -85,7 +88,7 @@ SEXP gforce(SEXP env, SEXP jsub, SEXP o, SEXP f, SEXP l, SEXP irowsArg) { // maybe better to malloc to avoid R's heap. This grp isn't global, so it doesn't need to be R_alloc const int *restrict fp = INTEGER_RO(f); - nBatch = MIN((nrow+1)/2, getDTthreads(nrow, true)*2); // *2 to reduce last-thread-home. TODO: experiment. The higher this is though, the bigger is counts[] + nBatch = single ? 1 : MIN((nrow+1)/2, getDTthreads(nrow, true)*2); // *2 to reduce last-thread-home. TODO: experiment. The higher this is though, the bigger is counts[] batchSize = MAX(1, (nrow-1)/nBatch); lastBatchSize = nrow - (nBatch-1)*batchSize; // We deliberate use, for example, 40 batches of just 14 rows, to stress-test tests. This strategy proved to be a good one as #3204 immediately came to light. @@ -222,6 +225,32 @@ void *gather(SEXP x, bool *anyNA) { double started=wallclock(); const bool verbose = GetVerbose(); + if (nBatch==1 && highSize==1 && irowslen==-1) { // hot path for single threads and ngrp<=65536 + // gather would be an identical copy of x, so return x itself after checking for NA (#3330). Callers only read gx. + switch (TYPEOF(x)) { + case LGLSXP: case INTSXP: { + const int *xp = INTEGER_RO(x); + for (int i=0; i Date: Sun, 18 Oct 2026 23:28:43 +0200 Subject: [PATCH 3/4] fix merge --- NEWS.md | 3 +-- inst/tests/tests.Rraw | 19 ------------------- src/gsumm.c | 4 ---- 3 files changed, 1 insertion(+), 25 deletions(-) diff --git a/NEWS.md b/NEWS.md index 0650e835d0..fc28800a1d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -94,8 +94,6 @@ 21. GForce grouping is again faster when running single-threaded, e.g. after `setDTthreads(1)`, [#3330](https://github.com/Rdatatable/data.table/issues/3330). The batching used to parallelize `gsum()` and `gmean()` is skipped when only one thread is available. Thanks @jangorecki for the benchmarks, @mattdowle for the report, and @ben-schwen for the fix. -22. GForce grouping is again faster when running single-threaded, e.g. after `setDTthreads(1)`, [#3330](https://github.com/Rdatatable/data.table/issues/3330). The batching used to parallelize `gsum()` and `gmean()` is skipped when only one thread is available. Thanks @jangorecki for the benchmarks, @mattdowle for the report, and @ben-schwen for the fix. - ### Notes 1. {data.table} now depends on R 3.5.0 (2018). @@ -118,6 +116,7 @@ 10. `fwrite()` returns a clearer error message when `na = data.frame()` is used, [#7866](https://github.com/Rdatatable/data.table/issues/7866). Thanks @mcol for the report and the fix. + ## data.table [v1.18.6.1](https://github.com/Rdatatable/data.table/milestone/47) (22 August 2026) ### NEW FEATURES diff --git a/inst/tests/tests.Rraw b/inst/tests/tests.Rraw index 1262120bc7..78123b73fe 100644 --- a/inst/tests/tests.Rraw +++ b/inst/tests/tests.Rraw @@ -22011,22 +22011,3 @@ test(2388.21, nafill(c(1, NA, NA, 5), type="nocb", limit=1, nan=NaN), c(1, NA, 5 test(2388.22, nafill(c(1L, NA, NA, 5L), type="nocb", limit=1), c(1L, NA, 5L, 5L)) test(2388.23, nafill(c("a", NA, NA), type="const", fill="x", limit=1), c("a", "x", NA)) test(2388.24, nafill(c("a", NA, NA, "z"), type="nocb", limit=1), c("a", NA, "z", "z")) - -# single-threaded GForce uses one batch and, with up to 65536 groups, gsum/gmean read x directly instead of gathering it, #3330 -old = setDTthreads(1L) -set.seed(3330L) -N = 2e5L -DT = data.table(g1=sample(100L, N, TRUE), g2=sample(70000L, N, TRUE), i=sample(c(1:10, NA), N, TRUE), k=sample(5L, N, TRUE), d=rnorm(N), c=complex(real=rnorm(N), imaginary=rnorm(N))) -DT[c(3L, 7L), d := NA][c(5L, 9L), c := NA] -test(2389.1, optimize=c(0L, Inf), DT[, .(sum(i), sum(i, na.rm=TRUE), mean(i), mean(i, na.rm=TRUE), sum(k), mean(k), sum(d), sum(d, na.rm=TRUE), mean(d), mean(d, na.rm=TRUE), sum(c), sum(c, na.rm=TRUE), mean(c), mean(c, na.rm=TRUE)), keyby=g1]) -test(2389.2, optimize=c(0L, Inf), DT[, .(sum(i), sum(i, na.rm=TRUE), mean(i), mean(i, na.rm=TRUE), sum(k), mean(k), sum(d), sum(d, na.rm=TRUE), mean(d), mean(d, na.rm=TRUE), sum(c), sum(c, na.rm=TRUE), mean(c), mean(c, na.rm=TRUE)), by=g2]) -test(2389.3, optimize=c(0L, Inf), DT[d > 0, .(sum(i), mean(i, na.rm=TRUE), sum(k), mean(d), sum(c)), by=g1]) -setkey(DT, g1) -test(2389.4, optimize=c(0L, Inf), DT[, .(sum(i), sum(i, na.rm=TRUE), mean(i), mean(i, na.rm=TRUE), sum(k), mean(k), sum(d), sum(d, na.rm=TRUE), mean(d), mean(d, na.rm=TRUE), sum(c), sum(c, na.rm=TRUE), mean(c), mean(c, na.rm=TRUE)), by=g1]) -test(2389.5, options=list(datatable.optimize=Inf), data.table(g=c(1L, 1L, 2L), x=c(.Machine$integer.max, 1L, 1L))[, sum(x), by=g], data.table(g=1:2, V1=c(2147483648, 1)), warning="coerced to 'numeric'") -if (test_bit64) { - DT[, i64 := as.integer64(i)] - test(2389.6, optimize=c(0L, Inf), DT[, .(sum(i64), sum(i64, na.rm=TRUE)), by=g1]) - test(2389.7, optimize=c(0L, Inf), DT[, .(sum(i64), sum(i64, na.rm=TRUE)), by=g2]) -} -setDTthreads(old) diff --git a/src/gsumm.c b/src/gsumm.c index a181089921..88d3cdb7de 100644 --- a/src/gsumm.c +++ b/src/gsumm.c @@ -225,11 +225,7 @@ void *gather(SEXP x, bool *anyNA) { double started=wallclock(); const bool verbose = GetVerbose(); -<<<<<<< HEAD if (nBatch==1 && highSize==1 && irowslen==-1) { // hot path for single threads and ngrp<=65536 -======= - if (nBatch==1 && highSize==1 && irowslen==-1) { ->>>>>>> 9bc4dd94f3fbe8efa1d9258944e097c9315fdaf9 // gather would be an identical copy of x, so return x itself after checking for NA (#3330). Callers only read gx. switch (TYPEOF(x)) { case LGLSXP: case INTSXP: { From fc67973e553ec6a42ff30d267a672c29a8720175 Mon Sep 17 00:00:00 2001 From: Benjamin Schwendinger Date: Mon, 19 Oct 2026 11:12:40 +0200 Subject: [PATCH 4/4] remove coercion from gmean --- NEWS.md | 2 +- inst/tests/optimize.Rraw | 1 - src/gsumm.c | 46 +++++++++++++++++++++++++++++++++++----- 3 files changed, 42 insertions(+), 7 deletions(-) diff --git a/NEWS.md b/NEWS.md index fc28800a1d..c0a0d1613d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -92,7 +92,7 @@ 20. `setnafill()` with `type = "locf"` or `type = "nocb"` on a 0-row table no longer segfaults, [#7910](https://github.com/Rdatatable/data.table/issues/7910). Thanks @MichaelChirico for the report and fix. -21. GForce grouping is again faster when running single-threaded, e.g. after `setDTthreads(1)`, [#3330](https://github.com/Rdatatable/data.table/issues/3330). The batching used to parallelize `gsum()` and `gmean()` is skipped when only one thread is available. Thanks @jangorecki for the benchmarks, @mattdowle for the report, and @ben-schwen for the fix. +21. GForce grouping is again faster when running single-threaded, e.g. after `setDTthreads(1)`, [#3330](https://github.com/Rdatatable/data.table/issues/3330). The batching used to parallelize `gsum()` and `gmean()` is skipped when only one thread is available. Also, `mean()` of an `integer` or `logical` column by group no longer coerces the whole column to `double` first, which made it slower than in v1.12.0 regardless of the number of threads. Thanks @jangorecki for the benchmarks, @mattdowle for the report, and @ben-schwen for the fix. ### Notes diff --git a/inst/tests/optimize.Rraw b/inst/tests/optimize.Rraw index 6afe56b65b..d345a06edb 100644 --- a/inst/tests/optimize.Rraw +++ b/inst/tests/optimize.Rraw @@ -508,4 +508,3 @@ 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) diff --git a/src/gsumm.c b/src/gsumm.c index 88d3cdb7de..d2adb81654 100644 --- a/src/gsumm.c +++ b/src/gsumm.c @@ -75,11 +75,12 @@ SEXP gforce(SEXP env, SEXP jsub, SEXP o, SEXP f, SEXP l, SEXP irowsArg) { if (length(tt)==1 && INTEGER(tt)[0]!=maxgrpn) internal_error(__func__, "o's maxgrpn attribute mismatches recalculated maxgrpn"); // # nocov } - // Single-threaded, the batching below is pure overhead (#3330). Use one batch and put as much of g into low as fits in - // uint16_t, so that for ngrp<=65536 there is one high bucket and gather/gsum/gmean reduce to direct loops. + // Single-threaded, the batching below is pure overhead (#3330). Use one batch and, when all of g fits in uint16_t low, + // a single high bucket so that gather/gsum/gmean reduce to direct loops. With more groups keep nb/2: the small blocks + // of ans per high bucket stay in cache and are faster than one large block. const bool single = getDTthreads(nrow, false)==1; int nb = nbit(ngrp-1); - bitshift = single ? MIN(nb, 16) : nb/2; // /2 so that high and low can be uint16_t, and no limit (even for nb=4) to stress-test. + bitshift = single && nb<=16 ? nb : nb/2; // /2 so that high and low can be uint16_t, and no limit (even for nb=4) to stress-test. // bitshift=MAX(nb-8,0); if (bitshift>16) bitshift=nb/2; // TODO: when we have stress-test off mode, do this mask = (1<>bitshift) + 1; @@ -624,8 +625,43 @@ SEXP gmean(SEXP x, SEXP narmArg) SEXP ans=R_NilValue; int protecti=0; switch(TYPEOF(x)) { - case LGLSXP: case INTSXP: - x = PROTECT(coerceVector(x, REALSXP)); protecti++; + case LGLSXP: case INTSXP: { + // accumulate the gathered ints into double directly; coerceVector(x, REALSXP) first would copy all of x (#3330) + const int *restrict gx = gather(x, &anyNA); + ans = PROTECT(allocVector(REALSXP, ngrp)); protecti++; + double *restrict ansp = REAL(ans); + memset(ansp, 0, ngrp*sizeof(double)); + const bool count = narm && anyNA; + int *restrict nna_counts = count ? calloc(ngrp, sizeof(*nna_counts)) : NULL; + if (count && !nna_counts) + error(_("Unable to allocate %d * %zu bytes for non-NA counts in gmean na.rm=TRUE"), ngrp, sizeof(*nna_counts)); // # nocov + #pragma omp parallel for num_threads(getDTthreads(highSize, false)) + for (int h=0; h