diff --git a/DESCRIPTION b/DESCRIPTION
index 2bb30c59..0971acfc 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -38,7 +38,9 @@ Imports:
     methods,
     statmod,
     parallel,
-    rlang
+    rlang,
+    matter,
+    BiocParallel
 Suggests: 
     BiocStyle,
     knitr,
diff --git a/NAMESPACE b/NAMESPACE
index e674f5c6..ba6c1354 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -22,6 +22,7 @@ export(MSstatsSelectFeatures)
 export(MSstatsSummarizationOutput)
 export(MSstatsSummarizeSingleLinear)
 export(MSstatsSummarizeSingleTMP)
+export(MSstatsSummarizeWithMultipleCores)
 export(MSstatsSummarizeWithSingleCore)
 export(MZMinetoMSstatsFormat)
 export(MaxQtoMSstatsFormat)
@@ -54,6 +55,12 @@ import(data.table)
 import(ggplot2)
 import(limma)
 import(lme4)
+importFrom(BiocParallel,bpisup)
+importFrom(BiocParallel,bplapply)
+importFrom(BiocParallel,bpnworkers)
+importFrom(BiocParallel,bpprogressbar)
+importFrom(BiocParallel,bpstart)
+importFrom(BiocParallel,bpstop)
 importFrom(MASS,rlm)
 importFrom(MSstatsConvert,DIANNtoMSstatsFormat)
 importFrom(MSstatsConvert,DIAUmpiretoMSstatsFormat)
@@ -74,9 +81,12 @@ importFrom(MSstatsConvert,SkylinetoMSstatsFormat)
 importFrom(MSstatsConvert,SpectronauttoMSstatsFormat)
 importFrom(Rcpp,sourceCpp)
 importFrom(data.table,as.data.table)
+importFrom(data.table,data.table)
+importFrom(data.table,fifelse)
 importFrom(data.table,melt)
 importFrom(data.table,rbindlist)
 importFrom(data.table,setDT)
+importFrom(data.table,setDTthreads)
 importFrom(data.table,uniqueN)
 importFrom(ggrepel,geom_text_repel)
 importFrom(gplots,heatmap.2)
@@ -97,6 +107,7 @@ importFrom(htmltools,tagList)
 importFrom(limma,squeezeVar)
 importFrom(lme4,lmer)
 importFrom(marray,maPalette)
+importFrom(matter,SnowfastParam)
 importFrom(methods,is)
 importFrom(parallel,clusterExport)
 importFrom(parallel,makeCluster)
diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R
new file mode 100644
index 00000000..1d4678fa
--- /dev/null
+++ b/R/MSstatsSummarizeWithMultipleCores.R
@@ -0,0 +1,476 @@
+#' Maximum RAM (in MB) this R process has ever used
+#'
+#' There are two paths:
+#'
+#' \enumerate{
+#'   \item \strong{Linux} (non-Windows and \code{/proc/self/status} exists):
+#'     read the number Linux itself keeps. \code{/proc/self/status} is a
+#'     plain-text file of \code{Key: value} lines describing this process; the
+#'     line \code{VmHWM:} holds the maximum RAM ever used, in kB. Preferred
+#'     because \code{/proc/self/clear_refs} (see \code{\link{.reset_max_rss}})
+#'     can reset it, so the reading reflects only recent work. The function 
+#'     pulls the number out of a line that looks like 
+#'     \code{"VmHWM:   123456 kB"}.
+#'   \item \strong{Everything else} (macOS, Windows, or Linux without
+#'     \code{/proc}): fall back to the compiled \code{peak_rss_mb()} in
+#'     \code{src/peak_rss.cpp}, which asks the operating system directly.
+#' }
+#'
+#' @return maximum RAM used, in MB
+#' @keywords internal
+#' @noRd
+.max_rss_mb <- function() {
+  if (.Platform$OS.type != "windows" && file.exists("/proc/self/status")) {
+    ln <- grep("^VmHWM:", readLines("/proc/self/status"), value = TRUE)
+    if (length(ln)) return(as.numeric(sub("\\D+(\\d+).*", "\\1", ln)) / 1024)
+  }
+  peak_rss_mb()
+}
+
+#' Reset the process maximum-RSS high-water mark to the current RSS
+#'
+#' The kernel-tracked maximum RSS (\code{VmHWM} on Linux, \code{ru_maxrss} on
+#' macOS/BSD, \code{PeakWorkingSetSize} on Windows) is a running maximum
+#' since process start. Without a reset, an earlier one-off allocation in the
+#' same R session (e.g. reading an 8 GB CSV) inflates every maximum-memory
+#' measurement taken afterwards even if that memory has since been freed,
+#' making it useless as a baseline for benchmarking this function's own
+#' memory use.
+#'
+#' Linux only: writing \code{"5"} to \code{/proc/self/clear_refs} resets
+#' \code{VmHWM} to the current RSS (kernel >= 4.0). macOS and Windows have no
+#' equivalent public facility to reset their maximum-RSS counters, so this is a
+#' no-op there, and \code{track_memory} reports on those platforms may still
+#' include memory used before this function was entered.
+#'
+#' @return invisible \code{TRUE} if the high-water mark was reset,
+#'   \code{FALSE} otherwise
+#' @keywords internal
+#' @noRd
+.reset_max_rss <- function() {
+    if (Sys.info()[["sysname"]] == "Linux" && file.exists("/proc/self/clear_refs")) {
+        reset_ok <- tryCatch({
+            cat("5", file = "/proc/self/clear_refs")
+            TRUE
+        }, error = function(e) FALSE, warning = function(w) FALSE)
+        return(invisible(reset_ok))
+    }
+    invisible(FALSE)
+}
+
+.print_memory_report <- function(function_name, checkpoints, worker_max_rss_mb = NULL,
+                                 elapsed = NULL) {
+    rule_width <- 65L
+    rule <- strrep("─", rule_width)
+    format_mb    <- function(x) if (is.na(x)) "    n/a" else sprintf("%7.1f", x)
+    format_delta <- function(now, previous_value) {
+        if (is.na(now) || is.na(previous_value)) return("")
+        sprintf("  (%+.1f MB)", now - previous_value)
+    }
+    lines <- c(
+        rule,
+        sprintf(" MSstats Memory Report — %s", function_name),
+        rule,
+        sprintf("  %-36s  %7s  %s", "Checkpoint", "RSS MB", "Delta")
+    )
+    previous_value <- NA_real_
+    for (checkpoint_name in names(checkpoints)) {
+        checkpoint_value <- checkpoints[[checkpoint_name]]
+        lines <- c(lines,
+                   sprintf("  %-36s  %s%s", checkpoint_name,
+                           format_mb(checkpoint_value),
+                           format_delta(checkpoint_value, previous_value)))
+        previous_value <- checkpoint_value
+    }
+    if (!is.null(worker_max_rss_mb)) {
+        observed_max_memory <- worker_max_rss_mb[!is.na(worker_max_rss_mb)]
+        if (length(observed_max_memory)) {
+            lines <- c(lines, "",
+                       sprintf("  Worker RSS  min / mean / max :  %.1f / %.1f / %.1f MB",
+                               min(observed_max_memory), mean(observed_max_memory),
+                               max(observed_max_memory)))
+        }
+    }
+    if (!is.null(elapsed))
+        lines <- c(lines, sprintf("  Total elapsed: %.1f s", as.numeric(elapsed)))
+    lines <- c(lines, rule)
+    message(paste(lines, collapse = "\n"))
+}
+
+#' Pack one protein slot into a double vector plus metadata
+#'
+#' @param protein_dt data.table rows for one protein (or protein x label) slot
+#' @param slot_index integer position of this slot in the global protein list
+#' @param all_runs character vector of all run names in global order
+#' @return list with elements \code{packed} (double vector) and \code{meta} (list)
+#' @keywords internal
+#' @noRd
+.pack_protein_slot <- function(protein_dt, slot_index, all_runs) {
+
+    n_runs <- length(all_runs)
+
+    has_peptide <- "PEPTIDE" %in% colnames(protein_dt)
+    feature_label_dt <- unique(protein_dt[, .(
+        FEATURE = as.character(FEATURE),
+        LABEL   = as.character(LABEL),
+        PEPTIDE = if (has_peptide) as.character(PEPTIDE) else as.character(FEATURE)
+    )])
+    data.table::setorder(feature_label_dt, FEATURE, LABEL)
+    n_feature_labels <- nrow(feature_label_dt)
+
+    feature_row_idx <- match(
+        paste(as.character(protein_dt$FEATURE), as.character(protein_dt$LABEL), sep = "\t"),
+        paste(feature_label_dt$FEATURE, feature_label_dt$LABEL, sep = "\t"))
+    run_col_idx_lookup <- seq_len(n_runs); names(run_col_idx_lookup) <- all_runs
+    run_col_idx  <- run_col_idx_lookup[as.character(protein_dt$RUN)]
+    row_is_valid <- !is.na(feature_row_idx) & !is.na(run_col_idx)
+
+    scatter_into_matrix <- function(col_vec, fill = NA_real_) {
+        m <- matrix(fill, nrow = n_feature_labels, ncol = n_runs)
+        m[cbind(feature_row_idx[row_is_valid], run_col_idx[row_is_valid])] <-
+            as.double(col_vec[row_is_valid])
+        m
+    }
+
+    new_abundance_mat <- scatter_into_matrix(protein_dt$newABUNDANCE)
+
+    has_ABUNDANCE <- "ABUNDANCE" %in% colnames(protein_dt)
+    abundance_mat <- if (has_ABUNDANCE) scatter_into_matrix(protein_dt$ABUNDANCE) else
+        matrix(NA_real_, n_feature_labels, n_runs)
+
+    has_censored <- "censored" %in% colnames(protein_dt)
+    censored_mat <- if (has_censored) scatter_into_matrix(as.double(protein_dt$censored)) else
+        matrix(0.0, n_feature_labels, n_runs)
+
+    has_cen <- "cen" %in% colnames(protein_dt)
+    event_mat <- if (has_cen) scatter_into_matrix(protein_dt$cen) else
+        matrix(NA_real_, n_feature_labels, n_runs)
+
+    has_anom <- "ANOMALYSCORES" %in% colnames(protein_dt) &&
+        !all(is.na(protein_dt$ANOMALYSCORES))
+    anomaly_scores_mat <- if (has_anom) scatter_into_matrix(protein_dt$ANOMALYSCORES) else
+        matrix(NA_real_, n_feature_labels, n_runs)
+
+    n_obs_by_feature_label <- protein_dt[, .(n_obs = as.double(n_obs[1L])),
+                      by = .(FEATURE = as.character(FEATURE),
+                             LABEL   = as.character(LABEL))]
+    feature_label_with_nobs <- n_obs_by_feature_label[feature_label_dt, on = c("FEATURE", "LABEL")]
+    data.table::setorder(feature_label_with_nobs, FEATURE, LABEL)
+    n_obs_vec   <- feature_label_with_nobs$n_obs
+
+    run_level_scalars <- protein_dt[, .(n_obs_run    = as.double(n_obs_run[1L]),
+                       prop_features = as.double(prop_features[1L])),
+                   by = .(RUN = as.character(RUN))]
+    run_scalars_all_runs <- run_level_scalars[data.table::data.table(RUN = all_runs), on = "RUN"]
+    n_obs_run_vec     <- run_scalars_all_runs$n_obs_run
+    prop_features_vec <- run_scalars_all_runs$prop_features
+
+    packed <- c(
+        as.double(slot_index),
+        as.vector(new_abundance_mat),
+        as.vector(abundance_mat),
+        as.vector(censored_mat),
+        as.vector(event_mat),
+        as.vector(anomaly_scores_mat),
+        n_obs_vec,
+        n_obs_run_vec,
+        prop_features_vec
+    )
+
+    meta <- list(
+        PROTEIN           = as.character(protein_dt$PROTEIN[1L]),
+        feature_label_dt  = as.data.frame(feature_label_dt),
+        runs              = all_runs,
+        n_feature_labels  = n_feature_labels,
+        n_runs            = n_runs,
+        is_labeled_ref    = "is_labeled_ref" %in% colnames(protein_dt) &&
+            isTRUE(any(protein_dt$is_labeled_ref, na.rm = TRUE)),
+        has_ABUNDANCE     = has_ABUNDANCE,
+        has_censored      = has_censored,
+        has_cen           = has_cen,
+        has_anom          = has_anom,
+        add_ref_covariate = "ref_covariate" %in% colnames(protein_dt)
+    )
+
+    list(packed = packed, meta = meta)
+}
+
+
+#' Reconstruct a per-protein data.table from a packed double vector
+#'
+#' @param packed double vector produced by \code{.pack_protein_slot}
+#' @param meta metadata list from \code{.pack_protein_slot}
+#' @return data.table compatible with \code{MSstatsSummarizeSingleTMP} /
+#'   \code{MSstatsSummarizeSingleLinear}
+#' @keywords internal
+#' @noRd
+.unpack_protein_slot <- function(packed, meta) {
+
+    n_feature_labels <- meta$n_feature_labels
+    n_runs           <- meta$n_runs
+    matrix_len       <- n_feature_labels * n_runs
+
+    cursor <- 2L
+    read_next_matrix <- function() {
+        m <- matrix(packed[cursor:(cursor + matrix_len - 1L)], nrow = n_feature_labels, ncol = n_runs)
+        cursor <<- cursor + matrix_len
+        m
+    }
+    new_abundance_mat  <- read_next_matrix()
+    abundance_mat      <- read_next_matrix()
+    censored_mat       <- read_next_matrix()
+    event_mat          <- read_next_matrix()
+    anomaly_scores_mat <- read_next_matrix()
+
+    n_obs_vec         <- packed[cursor:(cursor + n_feature_labels - 1L)]; cursor <- cursor + n_feature_labels
+    n_obs_run_vec     <- packed[cursor:(cursor + n_runs - 1L)]; cursor <- cursor + n_runs
+    prop_features_vec <- packed[cursor:(cursor + n_runs - 1L)]
+
+    feature_label_dt <- meta$feature_label_dt
+    runs   <- meta$runs
+    n_rows <- n_feature_labels * n_runs
+
+    protein_dt <- data.table::data.table(
+        PROTEIN       = rep(meta$PROTEIN, n_rows),
+        FEATURE       = rep(feature_label_dt$FEATURE,  times = n_runs),
+        LABEL         = rep(feature_label_dt$LABEL,    times = n_runs),
+        PEPTIDE       = rep(feature_label_dt$PEPTIDE,  times = n_runs),
+        RUN           = rep(runs,         each  = n_feature_labels),
+        newABUNDANCE  = as.vector(new_abundance_mat),
+        n_obs         = as.integer(rep(n_obs_vec,     times = n_runs)),
+        n_obs_run     = as.integer(rep(n_obs_run_vec, each  = n_feature_labels)),
+        prop_features = rep(prop_features_vec, each = n_feature_labels)
+    )
+
+    if (meta$has_ABUNDANCE) {
+        protein_dt[, ABUNDANCE := as.vector(abundance_mat)]
+    }
+
+    protein_dt[, censored := {
+        v <- as.vector(censored_mat)
+        if (meta$has_censored) as.logical(v > 0.5) else rep(FALSE, n_rows)
+    }]
+
+    if (meta$has_cen) {
+        protein_dt[, cen := as.vector(event_mat)]
+    }
+
+    protein_dt[, ANOMALYSCORES := as.vector(anomaly_scores_mat)]
+
+    if (meta$is_labeled_ref) {
+        protein_dt[, is_labeled_ref := (LABEL == "H")]
+        if (meta$add_ref_covariate) {
+            protein_dt[, ref_covariate := factor(
+                data.table::fifelse(LABEL == "L", as.character(RUN), "0"))]
+        }
+    }
+
+    protein_dt
+}
+
+
+#' Build the per-record worker closure for \code{MSstatsSummarizeWithMultipleCores}
+#'
+#' Defined at package top level so the closure only captures the scalar run
+#' parameters, not the caller's run-scale objects.
+#'
+#' @keywords internal
+#' @noRd
+.build_summarize_worker <- function(
+        use_TMP, impute, censored_symbol, remove50missing,
+        aft_iterations, equal_variance
+) {
+    unpack_fn        <- .unpack_protein_slot
+    use_TMP_         <- use_TMP
+    impute_          <- impute
+    censored_symbol_ <- censored_symbol
+    remove50missing_ <- remove50missing
+    aft_iterations_  <- aft_iterations
+    equal_variance_  <- equal_variance
+
+    function(record) {
+        meta   <- record$meta
+        protein_dt <- unpack_fn(record$packed, meta)
+        result <- if (use_TMP_) {
+            MSstatsSummarizeSingleTMP(
+                protein_dt, impute_, censored_symbol_,
+                remove50missing_, aft_iterations_)
+        } else {
+            MSstatsSummarizeSingleLinear(
+                protein_dt, impute_, censored_symbol_,
+                remove50missing_, aft_iterations_,
+                equal_variances = equal_variance_)
+        }
+        result
+    }
+}
+
+#' Per-worker maximum-RSS query task for \code{MSstatsSummarizeWithMultipleCores}
+#'
+#' @keywords internal
+#' @noRd
+.report_worker_max_rss <- function(i) {
+    list(worker = i, pid = Sys.getpid(), max_rss_mb = .max_rss_mb())
+}
+
+.warmup_worker <- function(i) {
+    library(MSstats, quietly = TRUE, warn.conflicts = FALSE)
+    data.table::setDTthreads(1)
+    NULL
+}
+
+
+#' Feature-level data summarization via socket-dispatched protein records
+#'
+#' @param input feature-level data processed by dataProcess subfunctions
+#' @param method summarization method: "linear" or "TMP"
+#' @param impute only for method = "TMP"; imputes censored values via AFT model when TRUE
+#' @param censored_symbol how censored values are encoded: 'NA', '0', or NULL for none
+#' @param remove50missing only for method = "TMP"; drops proteins missing >=50\% per peptide in every run
+#' @param equal_variance only for method = "linear"; assume equal variance among feature intensities
+#' @param numberOfCores number of cores for parallel processing (Linux/Mac only)
+#' @param aft_iterations number of AFT model iterations
+#' @param verbose whether to print verbose output
+#' @param BPPARAM optional \code{BiocParallelParam} instance
+#' @param track_memory whether to report per-worker maximum RSS memory usage.
+#'   On Linux, the process maximum RSS is reset on entry (via
+#'   \code{/proc/self/clear_refs}) so the reported baseline reflects memory
+#'   at function entry rather than a historical high-water mark inflated by
+#'   earlier operations in the session (e.g. reading a large input file).
+#'   This reset is Linux-only; on macOS and Windows there is no supported
+#'   facility to reset the maximum-RSS counter, so reports there may still
+#'   include memory used before this function was entered.
+#' @param max_proteins_per_worker caps protein records per \code{bplapply} task;
+#'   0 uses BiocParallel's default split, default is 50.
+#'
+#' @return A named list with one element per protein slot, keyed by protein
+#'   (or protein \eqn{\times} label) identifier.
+#'
+#' @importFrom matter SnowfastParam
+#' @importFrom BiocParallel bplapply bpstart bpstop bpisup bpnworkers bpprogressbar
+#' @importFrom data.table data.table fifelse setDTthreads
+#' @importFrom stats median
+#'
+#' @export
+MSstatsSummarizeWithMultipleCores <- function(
+        input,
+        method,
+        impute,
+        censored_symbol,
+        remove50missing,
+        equal_variance,
+        numberOfCores  = 1L,
+        aft_iterations = 90L,
+        verbose        = FALSE,
+        BPPARAM        = NULL,
+        track_memory   = FALSE,
+        max_proteins_per_worker = 50L
+) {
+    if (numberOfCores <= 1L && is.null(BPPARAM)) {
+        return(MSstatsSummarizeWithSingleCore(
+            input, method, impute, censored_symbol,
+            remove50missing, equal_variance, aft_iterations))
+    }
+
+    start_time <- proc.time()[["elapsed"]]
+    memory_checkpoints <- list()
+
+    if (track_memory) {
+        max_rss_reset <- .reset_max_rss()
+        if (!max_rss_reset) {
+            getOption("MSstatsLog")("INFO",
+                paste0("Maximum RSS reset on entry is only supported on Linux ",
+                       "(via /proc/self/clear_refs); on this platform, the ",
+                       "maximum-memory report below may include usage from ",
+                       "before this function was entered."))
+        }
+        memory_checkpoints[["baseline (max RSS reset on entry)"]] <- .max_rss_mb()
+    }
+
+    is_labeled_reference <- "is_labeled_ref" %in% colnames(input) &&
+        any(input$is_labeled_ref, na.rm = TRUE)
+    split_keys <- if (is_labeled_reference) list(input$PROTEIN) else
+        list(input$PROTEIN, input$LABEL)
+    protein_indices <- split(seq_len(nrow(input)), split_keys)
+    protein_ids     <- names(protein_indices)
+    num_proteins    <- length(protein_indices)
+
+    all_runs <- if (is.factor(input$RUN)) levels(input$RUN) else
+        as.character(sort(unique(input$RUN)))
+
+    getOption("MSstatsLog")("INFO",
+                            paste0("Packing ", num_proteins, " proteins × ",
+                                   length(all_runs), " runs into per-protein records"))
+
+    protein_records <- vector("list", num_proteins)
+    for (slot_index in seq_len(num_proteins)) {
+        protein_records[[slot_index]] <- .pack_protein_slot(
+            input[protein_indices[[slot_index]], ], slot_index, all_runs)
+    }
+
+    payload_mb <- sum(vapply(protein_records,
+                             function(r) length(r$packed), integer(1))) * 8 / 1024^2
+    getOption("MSstatsLog")("INFO",
+                            paste0("Dispatching via sockets (",
+                                   format(round(payload_mb, 1)),
+                                   " MB total packed payload; metadata sharded per task)"))
+
+    use_TMP <- identical(method, "TMP")
+
+    worker_fn <- .build_summarize_worker(
+        use_TMP, impute, censored_symbol, remove50missing,
+        aft_iterations, equal_variance)
+
+    if (is.null(BPPARAM)) {
+        tasks <- if (max_proteins_per_worker > 0L) {
+            as.integer(ceiling(num_proteins / max_proteins_per_worker))
+        } else {
+            0L
+        }
+        getOption("MSstatsLog")("INFO",
+                                paste0("Dispatching as ",
+                                       if (tasks > 0L) tasks else numberOfCores,
+                                       " task(s) (max_proteins_per_worker = ",
+                                       max_proteins_per_worker, ")"))
+        BPPARAM <- matter::SnowfastParam(
+            workers       = numberOfCores,
+            tasks         = tasks,
+            progressbar   = TRUE,
+            force.GC      = TRUE,
+            stop.on.error = FALSE)
+    }
+
+    started_here <- !BiocParallel::bpisup(BPPARAM)
+    if (started_here) {
+        BiocParallel::bpstart(BPPARAM)
+        on.exit(BiocParallel::bpstop(BPPARAM), add = TRUE)
+    }
+    show_progress <- BiocParallel::bpprogressbar(BPPARAM)
+
+    BiocParallel::bpprogressbar(BPPARAM) <- FALSE
+    BiocParallel::bplapply(
+        seq_len(BiocParallel::bpnworkers(BPPARAM)),
+        .warmup_worker, BPPARAM = BPPARAM)
+    BiocParallel::bpprogressbar(BPPARAM) <- show_progress
+
+    results <- BiocParallel::bplapply(protein_records, worker_fn, BPPARAM = BPPARAM)
+    names(results) <- protein_ids
+
+    worker_max_rss <- NULL
+    if (track_memory) {
+        BiocParallel::bpprogressbar(BPPARAM) <- FALSE
+        worker_max_rss <- BiocParallel::bplapply(
+            seq_len(BiocParallel::bpnworkers(BPPARAM)),
+            .report_worker_max_rss, BPPARAM = BPPARAM)
+        BiocParallel::bpprogressbar(BPPARAM) <- show_progress
+        memory_checkpoints[["parent maximum (main)"]] <- .max_rss_mb()
+        worker_max_rss_mb <- vapply(worker_max_rss, function(x) x$max_rss_mb, numeric(1L))
+        .print_memory_report(
+            "MSstatsSummarizeWithMultipleCores",
+            memory_checkpoints, worker_max_rss_mb,
+            elapsed = proc.time()[["elapsed"]] - start_time)
+    }
+
+    getOption("MSstatsLog")("INFO", "Summarization complete.")
+    results
+}
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 65f2c4fb..26dbc834 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -17,3 +17,7 @@ median_polish_summary <- function(x, eps = 0.01, maxiter = 10L) {
     .Call(`_MSstats_median_polish_summary`, x, eps, maxiter)
 }
 
+peak_rss_mb <- function() {
+    .Call(`_MSstats_peak_rss_mb`)
+}
+
diff --git a/R/dataProcess.R b/R/dataProcess.R
index 073d72b0..ef23f518 100755
--- a/R/dataProcess.R
+++ b/R/dataProcess.R
@@ -181,99 +181,6 @@ dataProcess = function(
     output
 }
 
-#' Feature-level data summarization with multiple cores
-#' 
-#' @param input feature-level data processed by dataProcess subfunctions
-#' @param method summarization method: "linear" or "TMP" 
-#' @param equal_variance only for summaryMethod = "linear". Default is TRUE. 
-#' Logical variable for whether the model should account for heterogeneous variation 
-#' among intensities from different features. Default is TRUE, which assume equal
-#' variance among intensities from features. FALSE means that we cannot assume 
-#' equal variance among intensities from features, then we will account for
-#' heterogeneous variation from different features.
-#' @param censored_symbol Indicates how censored missing values are encoded in
-#' the 'Intensity' column. 'NA' (default) treats all NA intensities as
-#' left-censored (i.e., below the limit of detection). '0' treats zero
-#' intensities as left-censored; in this case NA intensities are assumed to be
-#' missing at random and are not censored. Skyline output should use '0'. NULL
-#' assumes that all missing values are missing at random — no values are treated
-#' as censored, and imputation is disabled (impute is ignored).
-#' @param remove50missing only for summaryMethod = "TMP". TRUE removes the proteins
-#' where every run has at least 50\% missing values for each peptide. FALSE is default.
-#' @param impute only for summaryMethod = "TMP" and censored_symbol = 'NA' or '0'.
-#' TRUE (default) imputes censored missing values using an Accelerated Failure
-#' Time model. FALSE excludes censored observations from summarization entirely,
-#' treating them as missing at random; no imputed values are introduced.
-#' Has no effect when censored_symbol = NULL, since no values are considered censored.
-#' @param numberOfCores Number of cores for parallel processing. When > 1, 
-#' a logfile named `MSstats_dataProcess_log_progress.log` is created to 
-#' track progress. Only works for Linux & Mac OS. Default is 1.
-#' @param aft_iterations Number of iterations for AFT model fitting. Default is 90.
-#' 
-#' @importFrom parallel makeCluster parLapply stopCluster clusterExport
-#' 
-#' @return list of length one with run-level data.
-#' 
-MSstatsSummarizeWithMultipleCores = function(input, method, impute, censored_symbol,
-                              remove50missing, equal_variance, numberOfCores = 1,
-                              aft_iterations = 90) {
-    if (numberOfCores > 1) {
-        is_labeled_reference = "is_labeled_ref" %in% colnames(input) && any(input$is_labeled_ref, na.rm = TRUE)
-        if (is_labeled_reference) {
-            protein_indices = split(seq_len(nrow(input)), list(input$PROTEIN))
-        } else {
-            protein_indices = split(seq_len(nrow(input)), list(input$PROTEIN, input$LABEL))
-        }
-        num_proteins = length(protein_indices)
-        function_environment = environment()
-        cl = parallel::makeCluster(numberOfCores)
-        getOption("MSstatsLog")("INFO",
-                                "Starting the cluster setup for summarization")
-        parallel::clusterExport(cl, c("MSstatsSummarizeSingleTMP", 
-                                      "MSstatsSummarizeSingleLinear",
-                                      "input", "impute", "censored_symbol",
-                                      "remove50missing", "protein_indices", 
-                                      "equal_variance", "aft_iterations"), 
-                                envir = function_environment)
-        cat(paste0("Number of proteins to process: ", num_proteins), 
-            sep = "\n", file = "MSstats_dataProcess_log_progress.log")
-        if (method == "TMP") {
-            summarized_results = parallel::parLapply(cl, seq_len(num_proteins), function(i) {
-                if (i %% 100 == 0) {
-                    cat("Finished processing an additional 100 proteins", 
-                        sep = "\n", file = "MSstats_dataProcess_log_progress.log", append = TRUE)
-                }
-                single_protein = input[protein_indices[[i]],]
-                MSstatsSummarizeSingleTMP(
-                    single_protein, impute, censored_symbol, remove50missing,
-                    aft_iterations)
-            })
-        } else {
-            summarized_results = parallel::parLapply(cl, seq_len(num_proteins), function(i) {
-                if (i %% 100 == 0) {
-                    cat("Finished processing an additional 100 proteins", 
-                        sep = "\n", file = "MSstats_dataProcess_log_progress.log", append = TRUE)
-                }
-                single_protein = input[protein_indices[[i]],]
-                MSstatsSummarizeSingleLinear(
-                    single_protein,
-                    impute, 
-                    censored_symbol, 
-                    remove50missing,
-                    aft_iterations)
-            })
-        }
-        parallel::stopCluster(cl)
-        return(summarized_results)
-    } else {
-        return(MSstatsSummarizeWithSingleCore(input, method, impute, 
-                                              censored_symbol, 
-                                              remove50missing, 
-                                              equal_variance,
-                                              aft_iterations))
-    }
-}
-
 #' Feature-level data summarization with 1 core
 #' 
 #' @inheritParams MSstatsSummarizeWithMultipleCores
diff --git a/R/utils_censored.R b/R/utils_censored.R
index 73b81cb2..252708fe 100644
--- a/R/utils_censored.R
+++ b/R/utils_censored.R
@@ -84,7 +84,7 @@ MSstatsHandleMissing = function(input, summary_method, impute,
 #' possible values for left-censored data as the `time` input to the Surv function.
 #' @param input `data.table` in MSstats format
 #' @param censored_symbol censoredInt parameter to `dataProcess`
-#' @param remove50missing if TRUE, features with at least 50% missing values
+#' @param remove50missing if TRUE, features with at least 50\% missing values
 #' will be removed
 #' @keywords internal
 .setCensoredByThreshold = function(input, censored_symbol, remove50missing) {
diff --git a/R/utils_summarization.R b/R/utils_summarization.R
index 47bbef21..26225316 100644
--- a/R/utils_summarization.R
+++ b/R/utils_summarization.R
@@ -1,6 +1,6 @@
 #' Check if a protein can be summarized with TMP
 #' @param input data.table
-#' @param remove50missing if TRUE, proteins with more than 50% missing values
+#' @param remove50missing if TRUE, proteins with more than 50\% missing values
 #' in all runs will not be summarized
 #' @return data.table
 #' @keywords internal 
diff --git a/benchmark/benchmark_summarize_perf_selevsek.R b/benchmark/benchmark_summarize_perf_selevsek.R
new file mode 100644
index 00000000..5ea301aa
--- /dev/null
+++ b/benchmark/benchmark_summarize_perf_selevsek.R
@@ -0,0 +1,91 @@
+library(MSstats)
+
+# Performance expectations for MSstatsSummarizeWithMultipleCores on this dataset
+# (not enforced automatically here -- read the job output and compare by hand):
+#   - 4 cores should reduce wall time by at least 25% vs. 1 core
+#   - peak RSS per worker should stay under ~1GB (see the memory report below)
+
+northeastern_high_performance_cluster_file_path <-
+  "/projects/VitekLab/Data/MS/selevsek/before_summarization.csv"
+
+input <- data.table::fread(northeastern_high_performance_cluster_file_path)
+
+cat("=== MSstatsSummarizeWithMultipleCores performance check ===\n")
+cat("Expectation: >=25% reduction in wall time on 4 cores vs. 1 core.\n")
+cat("Expectation: peak RSS per worker should stay under ~1GB (see memory report below).\n\n")
+
+time_1core <- system.time(
+  result_1core <- MSstatsSummarizeWithSingleCore(
+    input, "TMP", TRUE, "NA", FALSE, TRUE
+  )
+)
+
+time_4core <- system.time(
+  result_4core <- MSstatsSummarizeWithMultipleCores(
+    input, "TMP", TRUE, "NA", FALSE, TRUE, 4, track_memory = TRUE,
+    max_proteins_per_worker = 200
+  )
+)
+
+elapsed_1core <- time_1core[["elapsed"]]
+elapsed_4core <- time_4core[["elapsed"]]
+speedup_pct <- (elapsed_1core - elapsed_4core) / elapsed_1core * 100
+
+cat(sprintf("1 core wall time:  %.2f s\n", elapsed_1core))
+cat(sprintf("4 core wall time:  %.2f s\n", elapsed_4core))
+cat(sprintf("Observed speedup:  %.1f%% (expectation: >= 25%%)\n", speedup_pct))
+cat("Compare the 'Worker RSS min / mean / max' line above against the ~1GB expectation.\n")
+
+cat("\n=== Verifying 1-core and 4-core results are identical ===\n")
+
+compare_df <- function(df1, df2) {
+  df1 <- as.data.frame(df1)
+  df2 <- as.data.frame(df2)
+
+  if (!setequal(names(df1), names(df2))) return(FALSE)
+
+  df2 <- df2[, names(df1), drop = FALSE]
+
+  # Convert factors to character so ordering/comparison is based on
+  # content, not on (possibly differing) factor level order
+  df1[] <- lapply(df1, function(x) if (is.factor(x)) as.character(x) else x)
+  df2[] <- lapply(df2, function(x) if (is.factor(x)) as.character(x) else x)
+
+  ord1 <- do.call(order, as.list(df1))
+  ord2 <- do.call(order, as.list(df2))
+
+  df1 <- df1[ord1, , drop = FALSE]
+  df2 <- df2[ord2, , drop = FALSE]
+
+  rownames(df1) <- NULL
+  rownames(df2) <- NULL
+
+  isTRUE(all.equal(df1, df2, tolerance = 1e-14))
+}
+
+n_proteins <- length(result_1core)
+
+protein_level_matches <- logical(n_proteins)
+for (i in seq_len(n_proteins)) {
+  protein_level_matches[i] <- compare_df(result_1core[[i]][[1]], result_4core[[i]][[1]])
+}
+
+feature_level_data_matches <- logical(n_proteins)
+for (i in seq_len(n_proteins)) {
+  feature_level_data_matches[i] <- compare_df(result_1core[[i]][[2]], result_4core[[i]][[2]])
+}
+
+n_protein_level_match <- sum(protein_level_matches)
+n_feature_level_data_match <- sum(feature_level_data_matches)
+
+cat(sprintf("Protein-level results matching: %d / %d\n", n_protein_level_match, n_proteins))
+cat(sprintf("Feature_level data matching:         %d / %d\n", n_feature_level_data_match, n_proteins))
+
+stopifnot(
+  "Protein-level results differ between 1 core and 4 cores" =
+    n_protein_level_match == n_proteins,
+  "Feature_level data differs between 1 core and 4 cores" =
+    n_feature_level_data_match == n_proteins
+)
+
+cat("1-core and 4-core results are identical.\n")
diff --git a/benchmark/config.slurm b/benchmark/config.slurm
index 28332e1a..28355762 100644
--- a/benchmark/config.slurm
+++ b/benchmark/config.slurm
@@ -27,7 +27,7 @@ remotes::install_github('Vitek-Lab/MSstats', ref = 'devel', lib = Sys.getenv('R_
 remotes::install_github('Vitek-Lab/MSstatsConvert', ref = 'master', lib = Sys.getenv('R_LIBS_USER')); \
 install.packages(c('dplyr', 'stringr', 'ggplot2'), lib = Sys.getenv('R_LIBS_USER'), repos = 'https://cloud.r-project.org')"
 
-R_SCRIPTS=("benchmark_Dowell2021-HEqe408_LFQ.R" "benchmark_Puyvelde2022-HYE5600735_LFQ.R")
+R_SCRIPTS=("benchmark_Dowell2021-HEqe408_LFQ.R" "benchmark_Puyvelde2022-HYE5600735_LFQ.R" "benchmark_summarize_perf_selevsek.R")
 
 for script in "${R_SCRIPTS[@]}"; do
     echo "Executing script: $script" >> job_output.txt  
diff --git a/man/MSstatsSummarizeSingleTMP.Rd b/man/MSstatsSummarizeSingleTMP.Rd
index 041ef76e..cd115723 100644
--- a/man/MSstatsSummarizeSingleTMP.Rd
+++ b/man/MSstatsSummarizeSingleTMP.Rd
@@ -15,22 +15,11 @@ MSstatsSummarizeSingleTMP(
 \arguments{
 \item{single_protein}{feature-level data for a single protein}
 
-\item{impute}{only for summaryMethod = "TMP" and censored_symbol = 'NA' or '0'.
-TRUE (default) imputes censored missing values using an Accelerated Failure
-Time model. FALSE excludes censored observations from summarization entirely,
-treating them as missing at random; no imputed values are introduced.
-Has no effect when censored_symbol = NULL, since no values are considered censored.}
+\item{impute}{only for method = "TMP"; imputes censored values via AFT model when TRUE}
 
-\item{censored_symbol}{Indicates how censored missing values are encoded in
-the 'Intensity' column. 'NA' (default) treats all NA intensities as
-left-censored (i.e., below the limit of detection). '0' treats zero
-intensities as left-censored; in this case NA intensities are assumed to be
-missing at random and are not censored. Skyline output should use '0'. NULL
-assumes that all missing values are missing at random — no values are treated
-as censored, and imputation is disabled (impute is ignored).}
+\item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none}
 
-\item{remove50missing}{only for summaryMethod = "TMP". TRUE removes the proteins
-where every run has at least 50\% missing values for each peptide. FALSE is default.}
+\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50\% per peptide in every run}
 
 \item{aft_iterations}{number of iterations for AFT model fitting}
 }
diff --git a/man/MSstatsSummarizeWithMultipleCores.Rd b/man/MSstatsSummarizeWithMultipleCores.Rd
index dfd7749d..b95b1e81 100644
--- a/man/MSstatsSummarizeWithMultipleCores.Rd
+++ b/man/MSstatsSummarizeWithMultipleCores.Rd
@@ -1,8 +1,8 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/dataProcess.R
+% Please edit documentation in R/MSstatsSummarizeWithMultipleCores.R
 \name{MSstatsSummarizeWithMultipleCores}
 \alias{MSstatsSummarizeWithMultipleCores}
-\title{Feature-level data summarization with multiple cores}
+\title{Feature-level data summarization via socket-dispatched protein records}
 \usage{
 MSstatsSummarizeWithMultipleCores(
   input,
@@ -11,8 +11,12 @@ MSstatsSummarizeWithMultipleCores(
   censored_symbol,
   remove50missing,
   equal_variance,
-  numberOfCores = 1,
-  aft_iterations = 90
+  numberOfCores = 1L,
+  aft_iterations = 90L,
+  verbose = FALSE,
+  BPPARAM = NULL,
+  track_memory = FALSE,
+  max_proteins_per_worker = 50L
 )
 }
 \arguments{
@@ -20,39 +24,38 @@ MSstatsSummarizeWithMultipleCores(
 
 \item{method}{summarization method: "linear" or "TMP"}
 
-\item{impute}{only for summaryMethod = "TMP" and censored_symbol = 'NA' or '0'.
-TRUE (default) imputes censored missing values using an Accelerated Failure
-Time model. FALSE excludes censored observations from summarization entirely,
-treating them as missing at random; no imputed values are introduced.
-Has no effect when censored_symbol = NULL, since no values are considered censored.}
-
-\item{censored_symbol}{Indicates how censored missing values are encoded in
-the 'Intensity' column. 'NA' (default) treats all NA intensities as
-left-censored (i.e., below the limit of detection). '0' treats zero
-intensities as left-censored; in this case NA intensities are assumed to be
-missing at random and are not censored. Skyline output should use '0'. NULL
-assumes that all missing values are missing at random — no values are treated
-as censored, and imputation is disabled (impute is ignored).}
-
-\item{remove50missing}{only for summaryMethod = "TMP". TRUE removes the proteins
-where every run has at least 50\% missing values for each peptide. FALSE is default.}
-
-\item{equal_variance}{only for summaryMethod = "linear". Default is TRUE. 
-Logical variable for whether the model should account for heterogeneous variation 
-among intensities from different features. Default is TRUE, which assume equal
-variance among intensities from features. FALSE means that we cannot assume 
-equal variance among intensities from features, then we will account for
-heterogeneous variation from different features.}
-
-\item{numberOfCores}{Number of cores for parallel processing. When > 1, 
-a logfile named `MSstats_dataProcess_log_progress.log` is created to 
-track progress. Only works for Linux & Mac OS. Default is 1.}
-
-\item{aft_iterations}{Number of iterations for AFT model fitting. Default is 90.}
+\item{impute}{only for method = "TMP"; imputes censored values via AFT model when TRUE}
+
+\item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none}
+
+\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50\% per peptide in every run}
+
+\item{equal_variance}{only for method = "linear"; assume equal variance among feature intensities}
+
+\item{numberOfCores}{number of cores for parallel processing (Linux/Mac only)}
+
+\item{aft_iterations}{number of AFT model iterations}
+
+\item{verbose}{whether to print verbose output}
+
+\item{BPPARAM}{optional \code{BiocParallelParam} instance}
+
+\item{track_memory}{whether to report per-worker maximum RSS memory usage.
+On Linux, the process maximum RSS is reset on entry (via
+\code{/proc/self/clear_refs}) so the reported baseline reflects memory
+at function entry rather than a historical high-water mark inflated by
+earlier operations in the session (e.g. reading a large input file).
+This reset is Linux-only; on macOS and Windows there is no supported
+facility to reset the maximum-RSS counter, so reports there may still
+include memory used before this function was entered.}
+
+\item{max_proteins_per_worker}{caps protein records per \code{bplapply} task;
+0 uses BiocParallel's default split, default is 50.}
 }
 \value{
-list of length one with run-level data.
+A named list with one element per protein slot, keyed by protein
+  (or protein \eqn{\times} label) identifier.
 }
 \description{
-Feature-level data summarization with multiple cores
+Feature-level data summarization via socket-dispatched protein records
 }
diff --git a/man/MSstatsSummarizeWithSingleCore.Rd b/man/MSstatsSummarizeWithSingleCore.Rd
index 4892dccb..fc711d42 100644
--- a/man/MSstatsSummarizeWithSingleCore.Rd
+++ b/man/MSstatsSummarizeWithSingleCore.Rd
@@ -19,29 +19,13 @@ MSstatsSummarizeWithSingleCore(
 
 \item{method}{summarization method: "linear" or "TMP"}
 
-\item{impute}{only for summaryMethod = "TMP" and censored_symbol = 'NA' or '0'.
-TRUE (default) imputes censored missing values using an Accelerated Failure
-Time model. FALSE excludes censored observations from summarization entirely,
-treating them as missing at random; no imputed values are introduced.
-Has no effect when censored_symbol = NULL, since no values are considered censored.}
+\item{impute}{only for method = "TMP"; imputes censored values via AFT model when TRUE}
 
-\item{censored_symbol}{Indicates how censored missing values are encoded in
-the 'Intensity' column. 'NA' (default) treats all NA intensities as
-left-censored (i.e., below the limit of detection). '0' treats zero
-intensities as left-censored; in this case NA intensities are assumed to be
-missing at random and are not censored. Skyline output should use '0'. NULL
-assumes that all missing values are missing at random — no values are treated
-as censored, and imputation is disabled (impute is ignored).}
+\item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none}
 
-\item{remove50missing}{only for summaryMethod = "TMP". TRUE removes the proteins
-where every run has at least 50\% missing values for each peptide. FALSE is default.}
+\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50\% per peptide in every run}
 
-\item{equal_variance}{only for summaryMethod = "linear". Default is TRUE. 
-Logical variable for whether the model should account for heterogeneous variation 
-among intensities from different features. Default is TRUE, which assume equal
-variance among intensities from features. FALSE means that we cannot assume 
-equal variance among intensities from features, then we will account for
-heterogeneous variation from different features.}
+\item{equal_variance}{only for method = "linear"; assume equal variance among feature intensities}
 
 \item{aft_iterations}{Number of iterations for AFT model fitting. Default is 90.}
 }
diff --git a/man/dot-getNonMissingFilterStats.Rd b/man/dot-getNonMissingFilterStats.Rd
index 7f8b8970..b4ca502f 100644
--- a/man/dot-getNonMissingFilterStats.Rd
+++ b/man/dot-getNonMissingFilterStats.Rd
@@ -9,13 +9,7 @@
 \arguments{
 \item{input}{data.table with data for a single protein}
 
-\item{censored_symbol}{Indicates how censored missing values are encoded in
-the 'Intensity' column. 'NA' (default) treats all NA intensities as
-left-censored (i.e., below the limit of detection). '0' treats zero
-intensities as left-censored; in this case NA intensities are assumed to be
-missing at random and are not censored. Skyline output should use '0'. NULL
-assumes that all missing values are missing at random — no values are treated
-as censored, and imputation is disabled (impute is ignored).}
+\item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none}
 }
 \value{
 data.table
diff --git a/man/dot-isSummarizable.Rd b/man/dot-isSummarizable.Rd
index b30b49cf..b13e68fb 100644
--- a/man/dot-isSummarizable.Rd
+++ b/man/dot-isSummarizable.Rd
@@ -9,7 +9,7 @@
 \arguments{
 \item{input}{data.table}
 
-\item{remove50missing}{if TRUE, proteins with more than 50% missing values
+\item{remove50missing}{if TRUE, proteins with more than 50\% missing values
 in all runs will not be summarized}
 }
 \value{
diff --git a/man/dot-runTukey.Rd b/man/dot-runTukey.Rd
index b6425a04..c1049973 100644
--- a/man/dot-runTukey.Rd
+++ b/man/dot-runTukey.Rd
@@ -15,16 +15,9 @@ subtracting the H value and adding back the H median, and only L results
 are returned. If FALSE (e.g. protein turnover), each label is summarized
 independently and results for all labels are returned.}
 
-\item{censored_symbol}{Indicates how censored missing values are encoded in
-the 'Intensity' column. 'NA' (default) treats all NA intensities as
-left-censored (i.e., below the limit of detection). '0' treats zero
-intensities as left-censored; in this case NA intensities are assumed to be
-missing at random and are not censored. Skyline output should use '0'. NULL
-assumes that all missing values are missing at random — no values are treated
-as censored, and imputation is disabled (impute is ignored).}
+\item{censored_symbol}{how censored values are encoded: 'NA', '0', or NULL for none}
 
-\item{remove50missing}{only for summaryMethod = "TMP". TRUE removes the proteins
-where every run has at least 50\% missing values for each peptide. FALSE is default.}
+\item{remove50missing}{only for method = "TMP"; drops proteins missing >=50\% per peptide in every run}
 }
 \value{
 data.table
diff --git a/man/dot-setCensoredByThreshold.Rd b/man/dot-setCensoredByThreshold.Rd
index b3ca6e63..45e66237 100644
--- a/man/dot-setCensoredByThreshold.Rd
+++ b/man/dot-setCensoredByThreshold.Rd
@@ -13,7 +13,7 @@ possible values for left-censored data as the `time` input to the Surv function.
 
 \item{censored_symbol}{censoredInt parameter to `dataProcess`}
 
-\item{remove50missing}{if TRUE, features with at least 50% missing values
+\item{remove50missing}{if TRUE, features with at least 50\% missing values
 will be removed}
 }
 \description{
diff --git a/src/Makevars.win b/src/Makevars.win
index 6660c7f9..eac6e524 100644
--- a/src/Makevars.win
+++ b/src/Makevars.win
@@ -1,3 +1,3 @@
 CXX_STD = CXX14
 PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS)
-PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
+PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -lpsapi
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 5275017b..a27a4301 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -67,12 +67,23 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// peak_rss_mb
+double peak_rss_mb();
+RcppExport SEXP _MSstats_peak_rss_mb() {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    rcpp_result_gen = Rcpp::wrap(peak_rss_mb());
+    return rcpp_result_gen;
+END_RCPP
+}
 
 static const R_CallMethodDef CallEntries[] = {
     {"_MSstats_get_estimable_fixed_random", (DL_FUNC) &_MSstats_get_estimable_fixed_random, 2},
     {"_MSstats_make_contrast_run_quant", (DL_FUNC) &_MSstats_make_contrast_run_quant, 6},
     {"_MSstats_get_linear_summary", (DL_FUNC) &_MSstats_get_linear_summary, 5},
     {"_MSstats_median_polish_summary", (DL_FUNC) &_MSstats_median_polish_summary, 3},
+    {"_MSstats_peak_rss_mb", (DL_FUNC) &_MSstats_peak_rss_mb, 0},
     {NULL, NULL, 0}
 };
 
diff --git a/src/peak_rss.cpp b/src/peak_rss.cpp
new file mode 100644
index 00000000..35bcf7fc
--- /dev/null
+++ b/src/peak_rss.cpp
@@ -0,0 +1,41 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+#if defined(_WIN32)
+#include <windows.h>
+#include <psapi.h>
+#elif defined(__unix__) || defined(__APPLE__)
+#include <sys/resource.h>
+#endif
+
+// Goal:   Ask the operating system for the most memory (RAM) this process has
+//         ever used. It is a running maximum, not the amount in use right now,
+//         so freeing memory does not lower it. Takes no arguments.
+//
+// Output: One number: maximum memory used, in MB. NA if the platform is not
+//         supported, or if the Windows query fails.
+//
+// The three branches below do the same thing with different OS calls; the unit
+// each one returns differs (Windows and macOS give bytes, Linux gives kB),
+// which is why the divisors are not identical.
+//
+// [[Rcpp::export]]
+double peak_rss_mb() {
+#if defined(_WIN32)
+    PROCESS_MEMORY_COUNTERS pmc;
+    if (GetProcessMemoryInfo(GetCurrentProcess(), &pmc, sizeof(pmc))) {
+        return (double) pmc.PeakWorkingSetSize / (1024.0 * 1024.0);
+    }
+    return NA_REAL;
+#elif defined(__unix__) || defined(__APPLE__)
+    struct rusage ru;
+    getrusage(RUSAGE_SELF, &ru);
+#ifdef __APPLE__
+    return (double) ru.ru_maxrss / (1024.0 * 1024.0);
+#else
+    return (double) ru.ru_maxrss / 1024.0;
+#endif
+#else
+    return NA_REAL;
+#endif
+}