diff --git a/NEWS.md b/NEWS.md
index cb4a6231fb..fc28800a1d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -92,6 +92,8 @@
 
 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).
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;
   highSize = ((ngrp-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<nrow; i++) if (xp[i]==NA_INTEGER) { *anyNA=true; break; }
+    } break;
+    case REALSXP: {
+      if (INHERITS(x, char_integer64)) {
+        const int64_t *xp = (const int64_t *)REAL_RO(x);
+        for (int i=0; i<nrow; i++) if (xp[i]==NA_INTEGER64) { *anyNA=true; break; }
+      } else {
+        const double *xp = REAL_RO(x);
+        for (int i=0; i<nrow; i++) if (ISNAN(xp[i])) { *anyNA=true; break; }
+      }
+    } break;
+    case CPLXSXP: {
+      const Rcomplex *xp = COMPLEX_RO(x);
+      for (int i=0; i<nrow; i++) if (ISNAN(xp[i].r) && ISNAN(xp[i].i)) { *anyNA=true; break; }
+    } break;
+    default : // # nocov
+      error(_("gather implemented for INTSXP, REALSXP, and CPLXSXP but not '%s'"), type2char(TYPEOF(x)));   // # nocov
+    }
+    if (verbose) { Rprintf(_("gather took %.3fs\n"), wallclock()-started); }
+    return (void *)DATAPTR_RO(x);
+  }
   switch (TYPEOF(x)) {
   case LGLSXP: case INTSXP: {
     const int *restrict thisx = INTEGER_RO(x);