diff --git a/UserManual/.DS_Store b/UserManual/.DS_Store
deleted file mode 100644
index 4bc01a0f..00000000
Binary files a/UserManual/.DS_Store and /dev/null differ
diff --git a/nCompiler/.DS_Store b/nCompiler/.DS_Store
deleted file mode 100644
index 5008ddfc..00000000
Binary files a/nCompiler/.DS_Store and /dev/null differ
diff --git a/nCompiler/R/Rexecution.R b/nCompiler/R/Rexecution.R
index fabae62a..c6ea54ba 100644
--- a/nCompiler/R/Rexecution.R
+++ b/nCompiler/R/Rexecution.R
@@ -14,6 +14,18 @@ parallel_for <- function(index, range, body, ...) {
 
 #' @export
 parallel_reduce <- function(f, x, init, ...) {
+  if(is.character(f)) {  # Not clear how to convert to char ...
+    operatorDef <- operatorDefEnv[[f]]
+    if(!is.null(operatorDef) && is.null(operatorDef$reduction))
+        stop("`", f, "` is not a valid reduction function/operator")
+  }
+  if(missing(init)) {
+    if(!is.character(f) || is.null(operatorDef) || is.null(operatorDef$reduction))
+      stop("`init` argument is missing and no default value provided for reduction function/operator")
+    init <- operatorDef$reduction
+  }
+  if(identical(f, "pairmin")) f <- "pmin"
+  if(identical(f, "pairmax")) f <- "pmax"
   Reduce(f, x, init)
 }
 
@@ -406,18 +418,18 @@ nRep <- function(x, ...) {
 }
 
 #' Sequence Generation
-#'
+#' 
 #' In a \code{nFunction}, \code{nSeq} is (mostly) equivalent to \code{base::seq}
 #'
-#' @details This function is similar to R's \code{\link{seq}} function, but
-#'   can be used in a nFunction and compiled using \code{nCompile}.
-#'
+#' @details This function is similar to R's \code{\link{seq}} function, but 
+#'   can be used in a nFunction and compiled using \code{nCompile}.  
+#' 
 #' @param from the starting value of the sequence.
-#'
+#' 
 #' @param to the ending value of the sequence.
-#'
+#' 
 #' @param by increment of the sequence
-#'
+#' 
 #' @param length.out desired length of the sequence
 #'
 #' @details
@@ -428,7 +440,6 @@ nRep <- function(x, ...) {
 #' There are no nCompiler versions of \code{seq.int}, \code{seq_along} or \code{seq_len}.
 
 #' @export
-#'
 nSeq <- function(...) {
   base::seq(...)
 }
@@ -513,6 +524,6 @@ nVar <- function(x) {
 #' Wrapper for sd
 #'
 #' @export
-nSd <- function(x) {
-  sd(x)
+nSd <- function(x) { 
+  sd(x) 
 }
diff --git a/nCompiler/R/compile_aaa_operatorLists.R b/nCompiler/R/compile_aaa_operatorLists.R
index 088ea49e..451705ef 100644
--- a/nCompiler/R/compile_aaa_operatorLists.R
+++ b/nCompiler/R/compile_aaa_operatorLists.R
@@ -12,6 +12,8 @@ returnTypeCodes <- list(
   promoteToDoubleOrAD = 6L,
   promoteNoLogical = 7L)
 
+liftedBlockOperatorsArg <- c("parallel_for" = 3, "parallel_reduce" = 1) # These are used for flagging when methods in lifted code block will need to have their object reference them explicitly via `obj__.<methodname>`. The values are the argID where the method needs to occur.
+
 returnTypeString2Code <- function(returnTypeString) {
   if(is.character(returnTypeString))
     do.call('switch', c(list("double"), returnTypeCodes))
@@ -162,6 +164,16 @@ assignOperatorDef(
   )
 )
 
+assignOperatorDef(
+  'nClass_method_in_lifted', # This is used for local method calls in the body of lifted code blocks (currently `parallel_{for,reduce}`).
+  list(
+    labelAbstractTypes = list(
+      handler = 'nFunction_or_method_call'),
+    cppOutput = list(
+      handler = 'nClass_method_in_lifted')
+  )
+)
+
 assignOperatorDef(
   'custom_default',
   list(
@@ -238,6 +250,10 @@ updateOperatorDef(
   c('nMatrix', 'nArray'),
   'labelAbstractTypes', 'returnTypeCode', returnTypeCodes$promote
 )
+updateOperatorDef(
+  'nMatrix',
+  'matchDef', val = function(value = 0, nrow = NA, ncol = NA, init = TRUE, fillZeros = TRUE, recycle = TRUE, type = 'double') {}
+)
 
 assignOperatorDef(
   'type_is',
@@ -301,20 +317,26 @@ assignOperatorDef(
 assignOperatorDef(
   c('parallel_for'),
   list(
+    matchDef = function(index, range, body, copyVars, shareVars, nThreads=0) {},
     labelAbstractTypes = list(
       handler = 'ParallelFor'),
     finalTransformations = list(
-      handler = "ParallelFor") ## Creates GeneralFor in the parallel_loop_body class
+      handler = "ParallelFor"), ## Creates GeneralFor in the parallel_loop_body class
+    cppOutput = list(
+      handler = 'ParallelExpr')  
   )
 )
 
 assignOperatorDef(
   c('parallel_reduce'),
-  list(
+  list(  
+    matchDef = function(operator, object, init, nThreads=0) {},
     labelAbstractTypes = list(
       handler = 'ParallelReduce'),
     finalTransformations = list(
-      handler = 'ParallelReduce')
+      handler = 'ParallelReduce'),
+    cppOutput = list(
+      handler = 'ParallelExpr')      
   )
 )
 
@@ -554,9 +576,11 @@ assignOperatorDef(
           )
         ),
         cppOutput = list(
-          handler = 'BinaryOrUnary')
+          handler = 'BinaryOrUnary'),
+        reduction = 0
     )
 )
+updateOperatorDef('-', 'reduction', val = NULL)
 
 assignOperatorDef(
   c('inprod'),
@@ -603,11 +627,13 @@ assignOperatorDef(
     labelAbstractTypes = list(
         handler = 'BinaryCwise',
         returnTypeCode = returnTypeCodes$promoteNoLogical),
-    cppOutput = list()
+    cppOutput = list(),
+    reduction = Inf
   )
 )
 updateOperatorDef('pairmax', 'cppOutput', 'cppString', 'std::max')
 updateOperatorDef('pairmin', 'cppOutput', 'cppString', 'std::min')
+updateOperatorDef('pairmax', 'reduction', val = -Inf)
 
 assignOperatorDef(
   c('pmin', 'pmax'),
@@ -984,7 +1010,8 @@ assignOperatorDef(
       )
     ),
     cppOutput = list(
-      handler = 'MidOperator')
+      handler = 'MidOperator'),
+    reduction = 1
   )
 )
 
diff --git a/nCompiler/R/compile_eigenization.R b/nCompiler/R/compile_eigenization.R
index 85853a36..d5606115 100644
--- a/nCompiler/R/compile_eigenization.R
+++ b/nCompiler/R/compile_eigenization.R
@@ -5,7 +5,8 @@
 eigenizeUseArgs <- c(
     list(
         setWhich = c(FALSE, TRUE),
-        setRepVectorTimes = c(FALSE, TRUE, TRUE)
+        setRepVectorTimes = c(FALSE, TRUE, TRUE),
+        parallel_reduce = c(FALSE, TRUE, TRUE)
         ))
 
 eigenizeEnv <- new.env()
diff --git a/nCompiler/R/compile_exprClass.R b/nCompiler/R/compile_exprClass.R
index 49565624..240daebd 100644
--- a/nCompiler/R/compile_exprClass.R
+++ b/nCompiler/R/compile_exprClass.R
@@ -218,6 +218,7 @@ wrapExprClassOperator <- function(code, funName, isName = FALSE, isCall = TRUE,
   newExpr
 }
 
+          
 insertIndexingBracket <- function(code, argID, index) {
   insertExprClassLayer(code, argID, 'index[')
   setArg(code$args[[argID]], 2, index)
diff --git a/nCompiler/R/compile_finalTransformations.R b/nCompiler/R/compile_finalTransformations.R
index 5ecf424f..9ec23d25 100644
--- a/nCompiler/R/compile_finalTransformations.R
+++ b/nCompiler/R/compile_finalTransformations.R
@@ -27,20 +27,7 @@ inFinalTransformationsEnv(
 
 inFinalTransformationsEnv(
   ParallelExpr <- function(parallel_expr_name, loop_body_name, auxEnv_field,
-                           code, symTab, auxEnv, info) {
-    copyVars <- eval(nDeparse(code$args[[4]], toR = TRUE), 
-                     envir = auxEnv$where)
-    shareVars <- eval(nDeparse(code$args[[5]], toR = TRUE),
-                      envir = auxEnv$where)
-    ## Look for a mangled argument name in nameSubList.
-    ## It is unfortunate to have to do this here instead of earlier
-    ## when other names are replaced, but here the names are given
-    ## as character objects (potentially from R evaluation).
-    copyVars <- replace_nameSubList(copyVars, auxEnv$nameSubList)
-    shareVars <- replace_nameSubList(shareVars, auxEnv$nameSubList)
-
-    code$args[[4]] <- copyVars ## This is no longer an exprClass
-    code$args[[5]] <- shareVars ## Ditto
+                           code, symTab, auxEnv, allVars, info) {
     auxEnv[[auxEnv_field]] <- c(auxEnv[[auxEnv_field]], code)
     ##  parallel_for(blocked_range<size_t>(0, n), parallel_loop_body(x));
     ## blocked_range_expr will be blocked_range<int>(start, end + 1)
@@ -62,11 +49,10 @@ inFinalTransformationsEnv(
                                    isName = FALSE, isLiteral = FALSE,
                                    isAssign = FALSE)
     setArg(parallel_expr, 1, blocked_range_expr)
-    ## loop_body_expr will be parallel_loop_body(var1, var2, etc.)
+    ## loop_body_expr will be parallel_loop_body_<id>(var1, var2, etc.)
     loop_body_expr <- exprClass$new(name = loop_body_name,
                                     isCall = TRUE,
                                     isName = FALSE, isLiteral = FALSE, isAssign = FALSE)
-    allVars <- c(copyVars, shareVars)
     for(iv in seq_along(allVars)) {
       ## Look for a mangled argument name in nameSubList.
       ## It is unfortunate to have to do this here instead of earlier
@@ -77,16 +63,103 @@ inFinalTransformationsEnv(
              exprClass$new(name = thisVar, 
                            isCall = FALSE, isName = TRUE, isLiteral = FALSE, isAssign = FALSE))
     }
+    if(length(code$aux$localMethods) || isTRUE(code$aux$liftedSelf))
+      setArg(loop_body_expr, iv+1, nParse('cppLiteral("*this")'))
     setArg(parallel_expr, 2, loop_body_expr)
     setArg(code$caller, code$callerArgID, parallel_expr)
+
+    nThreads_arg <- removeArg(code, 'nThreads')
+    setArg(parallel_expr, 3, nThreads_arg)
     NULL
   }
 )
 
 inFinalTransformationsEnv(
   ParallelFor <- function(code, symTab, auxEnv, info) {
-    ParallelExpr('parallel_for', 'parallel_loop_body', 'parallelContent', code,
-                 symTab, auxEnv, info)
+    nThreads_arg <- removeArg(code, 'nThreads')
+    ## TODO: not sure if we will do more work on arg matching such that
+    ## code$args[[4]] and code$args[[5]] will always exist and correspond to `copyVars` and `shareVars`
+    ## respectively for `parallel_for`. But if so, we might rework/simplify this.
+    ## Check for "" is because it seems valid to do: `parallel_for(i,1:5,{},}` or `parallel_for(i,1:5,{},,}`.
+    if(!'copyVars' %in% names(code$args) || nDeparse(code$args[['copyVars']]) == "") {
+      copyVars <- NULL
+    } else copyVars <- eval(nDeparse(code$args[['copyVars']], toR = TRUE), 
+                     envir = auxEnv$where)
+    if(!'shareVars' %in% names(code$args) || nDeparse(code$args[['shareVars']]) == "") {
+      shareVars <- NULL
+    } else shareVars <- eval(nDeparse(code$args[['shareVars']], toR = TRUE),
+                     envir = auxEnv$where)
+    if(any(shareVars %in% copyVars))
+              stop(exprClassProcessingErrorMsg(
+      code,
+      paste('In finalTransformations handler ParallelExpr:',
+            'arguments `shareVars` and `copyVars` to `parallel_for`',
+            'both contain the same variable')), call. = FALSE)
+    ## Look for a mangled argument name in nameSubList.
+    ## It is unfortunate to have to do this here instead of earlier
+    ## when other names are replaced, but here the names are given
+    ## as character objects (potentially from R evaluation).
+    copyVars <- replace_nameSubList(copyVars, auxEnv$nameSubList)
+    shareVars <- replace_nameSubList(shareVars, auxEnv$nameSubList)
+
+    ## Add default vars:
+    ## Any argument, class member variable, nFunction local variable by default is shared.
+    ## Any local variable in the loop body by default is copied.
+    vars <- all.vars(code$args[[3]]$Rexpr)
+    nms <- all.names(code$args[[3]]$Rexpr)
+
+    vars2 <- vars[!vars %in% c("self", nDeparse(code$args[[1]]))]  # Last item is index variable.
+    # Omit vars only referenced by `self`.
+    # This reliance on ordering of result of `all.names` feels fragile.
+    if("self" %in% nms) {
+      nonSelfNames <- nms[-(which(nms == "self") + 1)]
+    } else nonSelfNames <- nms
+    vars2 <- vars2[vars2 %in% nonSelfNames]
+    
+    inST <- vars2 %in% c(symTab$getSymbolNames(), symTab$parentST$getSymbolNames())
+    defaultCopyVars <- code$aux$localVars  # Local vars in for loop body.
+    defaultCopyVars <- defaultCopyVars[!defaultCopyVars %in% shareVars]
+    defaultShareVars <- vars2[inST]     # All other vars.
+    defaultShareVars <- defaultShareVars[!defaultShareVars %in% code$aux$localVars]
+    defaultShareVars <- defaultShareVars[!defaultShareVars %in% copyVars]
+
+    ## Find nClass objects (if methods are used; members would have been found above).
+    nms <- nms[!nms %in% c(vars, "self")]
+    objects <- nms[nms %in% c(symTab$getSymbolNames(), symTab$parentST$getSymbolNames())]
+    ## Make sure the items are actually nClass objects.
+    if(length(objects))  
+      objects <- objects[sapply(objects, 
+         function(x) !is.null(symTab$getSymbol(x)$NCgenerator) || !is.null(symTab$parentST$getSymbol(x)$NCgenerator))]
+    
+    shareVars <- unique(c(shareVars, defaultShareVars, objects))
+    copyVars  <- unique(c(copyVars, defaultCopyVars))
+      
+    ## NULL cannot hold a position in `code$args`.
+    if(is.null(copyVars)) copyVars <- character(0)
+    if(is.null(shareVars)) shareVars <- character(0)
+        
+    code$args[[4]] <- copyVars ## This is no longer an exprClass
+    code$args[[5]] <- shareVars ## Ditto
+    setArg(code, 6, nThreads_arg)
+    names(code$args)[4:6] <- c('copyVars','shareVars','nThreads')
+
+      
+    ## We have already found the local method calls and set the `opInfo$case` to be 'nClass_method_in_lifted',
+    ## such that C++ calls to the method will be handled by cppOutput handler.
+    ## The following checks for such methods in a different way (so perhaps worry an inconsistency could arise).
+    ## Perhaps there is a better way to get this information.
+    ## This information is used to ensure that the self object is passed into the lifted TBB code.
+    ## Currently we don't use the actual identified `localMethods` values, just whether there are any.
+    nms <- all.names(code$args[[3]]$Rexpr)
+    code$aux$localMethods <- nms[nms %in% c(names(auxEnv$where$public_methods), names(auxEnv$where$private_methods))]
+    if("self" %in% nms)  # This is needed to catch case of `self$<field>` but is redundant with the above line for `self$<method>`.
+      code$aux$liftedSelf <- TRUE
+    code$aux$class <- auxEnv$where$classname
+
+    code$aux$bodyName <- parallelForBodyLabelMaker()
+      
+    ParallelExpr('parallel_for', code$aux$bodyName, 'parallelContent', code,
+                 symTab, auxEnv, allVars = c(copyVars, shareVars), info)
   }
 )
 
@@ -104,11 +177,25 @@ inFinalTransformationsEnv(
     ## shareVars <- replace_nameSubList(shareVars, auxEnv$nameSubList)
     ## code$args[[4]] <- copyVars ## This is no longer an exprClass
     ## code$args[[5]] <- shareVars ## Ditto
-            
-    ## remove the vector and initial value arg and save for later
-    vector_arg <- removeArg(code, 2)
-    ## TODO: don't remove the init arg unless isTRUE(code$caller$isAssign)
-    init_arg <- removeArg(code, 2)
+
+    ## We have already found the local method calls and set the `opInfo$case` to be 'nClass_method_in_lifted',
+    ## such that C++ calls to the method will be handled by cppOutput handler.
+    ## The following checks for such methods in a different way (so perhaps worry an inconsistency could arise).
+    ## Perhaps there is a better way to get this information.
+    ## This information is used to ensure that the self object is passed into the lifted TBB code.
+    ## Currently we don't use the actual identified `localMethods` values, just whether there are any.
+    if ("self" %in% all.names(code$Rexpr))
+      code$aux$liftedSelf <- TRUE
+    nm <- code$args[[1]]$Rexpr
+    if(is.character(nm) && nm %in% c(names(auxEnv$where$public_methods), names(auxEnv$where$private_methods)))
+      code$aux$localMethods <- nm else code$aux$localMethods <- character(0)
+      
+    code$aux$class <- auxEnv$where$classname
+    code$aux$bodyName <- parallelReduceBodyLabelMaker()
+    ## remove the vector, initial value, and nThreads args and save for later
+    vector_arg <- removeArg(code, 'object')
+    init_arg <- removeArg(code, 'init')
+    nThreads_arg <- removeArg(code, 'nThreads')  
     ## add an index var
     index_arg <- exprClass$new(name = 'i__', isName = TRUE, isCall = FALSE,
                                isLiteral = FALSE, isAssign = FALSE)
@@ -128,13 +215,17 @@ inFinalTransformationsEnv(
     setArg(colon, 1, exprClass$new(name = 1, isLiteral = TRUE, isCall = FALSE,
                                    isName = FALSE, isAssign = FALSE))
     size_expr <- setArg(
-      colon, 2, nParse(paste0('cppLiteral("', vector_arg$name, '.size();")')))
+      colon, 2, nParse(paste0('cppLiteral("', vector_arg$name, '.size()")')))
     ## make the vector an argument of the reduce op and index it
     reduce_op <- code$args[[3]]
-    setArg(reduce_op, 1, copyExprClass(vector_arg))
-    insertIndexingBracket(reduce_op, 1, copyExprClass(index_arg))
+    inc <- 0  
+    if(reduce_op$name == 'chainedCall')
+      inc <- 1
+          
+    setArg(reduce_op, 1+inc, copyExprClass(vector_arg))
+    insertIndexingBracket(reduce_op, 1+inc, copyExprClass(index_arg))
     ## the other arg to the reduce op is a local aggregation var called 'val__'
-    val <- setArg(reduce_op, 2, exprClass$new(name = 'val__', isName = TRUE,
+    val <- setArg(reduce_op, 2+inc, exprClass$new(name = 'val__', isName = TRUE,
                                               isCall = FALSE,
                                               isLiteral = FALSE,
                                               isAssign = FALSE))
@@ -162,40 +253,102 @@ inFinalTransformationsEnv(
       symTab$addSymbol(value_type)
     }
 
-    ## The class name is hard-wired expecting only a single case of parallel
-    ## reduce content.
-    ## TO-DO: generalize the name with unique identifier.
+    instName <- sub("_body", "_inst__", code$aux$bodyName)
+
+    inputVar <-  eval(nDeparse(code$args[[4]], toR = TRUE), 
+                   envir = auxEnv$where)
+    outputVar <- eval(nDeparse(code$args[[5]], toR = TRUE),
+                    envir = auxEnv$where)
+    ## Look for a mangled argument name in nameSubList.
+    ## It is unfortunate to have to do this here instead of earlier
+    ## when other names are replaced, but here the names are given
+    ## as character objects (potentially from R evaluation).
+    inputVar <- replace_nameSubList(inputVar, auxEnv$nameSubList)
+
+    nms <- all.vars(code$Rexpr)
+    nClass_object <- nms[nms %in% c(symTab$getSymbolNames(), symTab$parentST$getSymbolNames()) &
+                  !nms %in% c(inputVar, "self")]
+    if(length(nClass_object) > 1)
+            stop(exprClassProcessingErrorMsg(
+                code$Rexpr,
+                paste('In finalTransformations handler ParallelReduce:',
+                      'Unexpectedly found multiple objects in parallel_reduce reduction function')),
+                call. = FALSE)
+    ## Make sure the items are actually nClass objects.
+    if(length(nClass_object) && is.null(symTab$getSymbol(nClass_object)$NCgenerator) &&
+       is.null(symTab$parentST$getSymbol(nClass_object)$NCgenerator))
+      nClass_object <- character(0)
+    
+    ## TODO: consider reworking how we handle these items as it doesn't map cleanly onto the
+    ## `args`, which was really set up for `parallel_for`.
+    code$args[[4]] <- inputVar  ## This is no longer an exprClass
+    code$args[[5]] <- outputVar ## Ditto
+    code$args[[6]] <- nClass_object    ## Ditto
+    setArg(code, 7, nThreads_arg)
+    names(code$args)[4:7] <- c('input','output','nClass_object','nThreads')    
+
     ParallelExpr('parallel_reduce',
-                 'parallel_reduce_body parallel_reduce_inst__',
-                 'parallelReduceContent', code, symTab, auxEnv, info)
-
-    if (isTRUE(code$caller$isAssign)) {
-      assign_argID <- code$caller$callerArgID
-      parallel_reduce_expr <- removeArg(code$caller, 2)
-      ## the instantiation of the parallel_reduce_body object will happen
-      ## before the call to parallel_reduce
-      instance_expr <- removeArg(parallel_reduce_expr, 2)
-      ## the second argument should be the initial value provided by the user
-      setArg(instance_expr, 2, init_arg)
-      ## TODO: this doesn't have the effect I hoped for... is there a way to
-      ## add type annotation to a call (such as object instantiation)?
-      instance_expr$type <- symbolBase$new(name = 'parallel_reduce_body',
-                                           type = 'parallel_reduce_body')
-      ## the parallel_reduce_body instance name is the second arg to the
-      ## parallel_reduce call (note that this isn't an exprClass)
-      setArg(parallel_reduce_expr, 2,
-             exprClass$new(name = 'parallel_reduce_inst__', isName = TRUE,
-                           isCall = FALSE, isLiteral = FALSE, isAssign = FALSE))
-      ## move the parallel_reduce_body instantiation to before the assignment
-      insertArg(code$caller$caller, assign_argID, instance_expr)
-      ## put the parallel_reduce call between the parallel_reduce_body
-      ## instantiation and the assign
-      insertArg(code$caller$caller, assign_argID + 1, parallel_reduce_expr)
-      ## now the RHS of the assign is the aggregation value after the
-      ## parallel_reduce
-      setArg(code$caller, 2,
-             nParse(paste0('cppLiteral("parallel_reduce_inst__.value__;")')))
+                 paste(code$aux$bodyName, instName, collapse = ' '),
+                 'parallelReduceContent', code, symTab, auxEnv, allVars = c(inputVar, outputVar, nClass_object), info)
+  
+    outerCall <- code$caller
+    level <- 1  
+    while(!isTRUE(outerCall$isAssign) && !outerCall$name == "return" && !outerCall$name == "{") {
+      outerCall <- outerCall$caller  # Find correct level to insert the reduction code.
+      level <- level + 1
+      if(level > 100)  ## Not sure what situation could lead to this.
+        stop(exprClassProcessingErrorMsg(
+           code$Rexpr,
+           paste('In finalTransformations handler ParallelReduce:',
+                 'Unexpected levels of nesting in use of parallel_reduce')),
+           call. = FALSE)
     }
+    if(outerCall$name == "{")  # No assignment or return. Handle these gracefully but no known use cases.
+      if(outerCall$args[[1]]$name == 'parallel_reduce') {  # A lone `parallel_reduce()` 
+        ## Add layer so that parallel_reduce call is within a call so can be handled
+        ## as other cases are handled.
+        code$caller <- wrapInExprClass(code$caller, "{")
+        setCaller(code, code$caller$args[[1]], 1)
+      } else outerCall <- outerCall$args[[1]]  # A case like `3 + parallel_reduce()`
+      
+    code$aux$init <- init_arg
+    assign_argID <- outerCall$callerArgID  # Always 1, presumably.
+
+    ## Check for `tbb::blocked_range<int>` handles cases such as `parallel_reduce() + parallel_reduce()`,
+    ## distinguishing which one is currently being processed.
+    reduce_argID <- which(sapply(code$caller$args, function(x) 
+        x$name == "parallel_reduce" && x$args[[1]]$name == "tbb::blocked_range<int>"))
+      
+    if(length(reduce_argID) != 1)  
+      stop(exprClassProcessingErrorMsg(
+         code$Rexpr,
+         paste('In finalTransformations handler ParallelReduce:',
+               'Unable to process code - missing or too many uses of parallel_reduce')),
+         call. = FALSE)
+    parallel_reduce_expr <- removeArg(code$caller, reduce_argID)
+    ## the instantiation of the parallel_reduce_body object will happen
+    ## before the call to parallel_reduce
+    instance_expr <- removeArg(parallel_reduce_expr, 2)
+    ## the second argument should be the initial value provided by the user
+    setArg(instance_expr, 2, init_arg)
+    ## TODO: this doesn't have the effect I hoped for... is there a way to
+    ## add type annotation to a call (such as object instantiation)?
+    instance_expr$type <- symbolBase$new(name = code$aux$bodyName,
+                                         type = 'parallel_reduce_body')
+    ## the parallel_reduce_body instance name is the second arg to the
+    ## parallel_reduce call (note that this isn't an exprClass)
+    insertArg(parallel_reduce_expr, 2,
+           exprClass$new(name = instName, isName = TRUE,
+                         isCall = FALSE, isLiteral = FALSE, isAssign = FALSE))
+    ## move the parallel_reduce_body instantiation to before the assignment    
+    insertArg(outerCall$caller, assign_argID, instance_expr)
+    ## put the parallel_reduce call between the parallel_reduce_body
+    ## instantiation and the assign
+    insertArg(outerCall$caller, assign_argID + 1, parallel_reduce_expr)
+    ## now the RHS of the assign is the aggregation value after the
+    ## parallel_reduce
+    insertArg(code$caller, reduce_argID,
+              nParse(paste0('cppLiteral("', instName, '.value__")')))
     NULL
   }
 )
diff --git a/nCompiler/R/compile_generateCpp.R b/nCompiler/R/compile_generateCpp.R
index 9bdd5b54..bc3d5fbd 100644
--- a/nCompiler/R/compile_generateCpp.R
+++ b/nCompiler/R/compile_generateCpp.R
@@ -258,6 +258,20 @@ inGenCppEnv(
   }
 )
 
+inGenCppEnv(
+  nClass_method_in_lifted <- function(code, symTab) {
+    cpp_code_name <- code$aux$cachedOpInfo$obj_internals$cpp_code_name
+    paste0(selfNameInLiftedBlock, ".", cpp_code_name,
+           '(', paste0(unlist(lapply(code$args,
+                                     compile_generateCpp,
+                                     symTab,
+                                     asArg = TRUE) ),
+                       collapse = ', '),
+           ')' )
+  }
+)
+
+
 inGenCppEnv(
   nClass_constructor <- function(code, symTab) {
     paste0("nClass_builder<" , code$type$name ,">()")
@@ -277,7 +291,7 @@ inGenCppEnv(
 )
 
 inGenCppEnv(
-  MidOperator <- function(code, symTab) {
+    MidOperator <- function(code, symTab) {
     if(length(code$args) != 2) stop('Error: expecting 2 arguments for operator ',code$name)
     if(is.null(code$caller)) useParens <- FALSE
     else {
@@ -411,15 +425,20 @@ inGenCppEnv(
   ## Member(A, x) -> A.x
   Member <- function(code, symTab, connector = '.') {
     isSelf <- code$args[[1]]$name == 'self' && inherits(code$args[[1]]$type, "symbolSelf")
-    objOutput <- if(isSelf) 'this' else compile_generateCpp(code$args[[1]], symTab)
-    paste0( '(',
-           objOutput,
-           ')', connector, code$args[[2]]$name)
+    isLiftedSelf <- code$args[[1]]$name == 'liftedSelf' && inherits(code$args[[1]]$type, "symbolSelf")
+    if(!isLiftedSelf) {    
+      objOutput <- if(isSelf) 'this' else compile_generateCpp(code$args[[1]], symTab)
+      return(paste0( '(',
+                    objOutput,
+                    ')', connector, code$args[[2]]$name)
+             )
+    }
+    return(paste0(selfNameInLiftedBlock, ".", code$args[[2]]$name))
   }
 )
 
 inGenCppEnv(
-  ## Member(A, x) -> A.x
+  ## PtrMember(A, x) -> A->x
   PtrMember <- function(code, symTab) {
     Member(code, symTab, connector = '->')
   }
@@ -711,3 +730,29 @@ inGenCppEnv(
     res
   }
 )
+
+inGenCppEnv(
+  ParallelExpr <- function(code, symTab) {
+    nThreads_arg <- removeArg(code, 3)
+    paste0('{',
+      paste0('tbb::global_control gc(tbb::global_control::max_allowed_parallelism, getNumThreads(',
+             compile_generateCpp(nThreads_arg, symTab),
+             '));'),
+      paste0(eval(call("AsIs", code, symTab), envir = genCppEnv), ';'),
+      '}', collapse = '\n')
+  }
+)
+
+inGenCppEnv(
+  ParallelExpr <- function(code, symTab) {
+    nThreads_arg <- removeArg(code, 3)
+    paste0('{',
+           "depth_counter mydc;",   ## This increments TBB_DEPTH and destructor will decrement.
+           paste0('tbb::global_control gc(tbb::global_control::max_allowed_parallelism, getNumThreads(',
+                  compile_generateCpp(nThreads_arg, symTab),
+                  '));'),
+           AsIs(code, symTab),
+           ';}',
+           collapse = '\n')
+  }
+)
diff --git a/nCompiler/R/compile_labelAbstractTypes.R b/nCompiler/R/compile_labelAbstractTypes.R
index 9e62a7df..60645cc7 100644
--- a/nCompiler/R/compile_labelAbstractTypes.R
+++ b/nCompiler/R/compile_labelAbstractTypes.R
@@ -320,9 +320,16 @@ inLabelAbstractTypesEnv(
     function(code, symTab, auxEnv, handlingInfo) {
       ## TO-DO: Check for exactly 2 arguments
       # Special handling for "new"
+      
       inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv,
                                             handlingInfo,
                                             useArgs = c(TRUE, FALSE))
+
+      # We need to handle `self$` differently in parallel body code
+      # as the code is lifted out of the class method in a new class.
+      if(isTRUE(auxEnv$inParallelBody) && code$args[[1]]$name == "self")
+        code$args[[1]]$name <- "liftedSelf"
+
       ## TO-DO: Check that LHS type is symbolNC or symbolNCgenerator
       ## TO-DO: Improve these error messages
       if(inherits(code$args[[1]]$type, "symbolNCgenerator")) {
@@ -574,7 +581,7 @@ inLabelAbstractTypesEnv(
 )
 
 inLabelAbstractTypesEnv(
-  InitData <- function(code, symTab, auxEnv, handlingInfo) {
+    InitData <- function(code, symTab, auxEnv, handlingInfo) {
     ## TODO: handle 'init' arg
     ## defaults:
     ## n{Numeric|Integer|Logical}(length = 0, value = 0, init = TRUE)
@@ -770,11 +777,12 @@ inLabelAbstractTypesEnv(
 
 inLabelAbstractTypesEnv(
   ParallelFor <- function(code, symTab, auxEnv, handlingInfo) {
-    if(length(code$args) != 5) 
+    if(length(code$args) < 3 || !identical(names(code$args)[1:3],
+                                               c('index','range','body')))
       stop(exprClassProcessingErrorMsg(
         code,
         paste('In labelAbstractTypes handler ParallelFor:',
-              'expected 5 arguments to a parallel_for-loop')), call. = FALSE)
+              'expected arguments `index`, `range`, and `body` to a parallel_for-loop')), call. = FALSE)
     ## first handle type of the indexing variable
     if(!inherits(code$args[[2]], 'exprClass'))
       stop(
@@ -794,71 +802,148 @@ inLabelAbstractTypesEnv(
     if (!symTab$symbolExists(code$args[[1]]$name, inherits = TRUE))
       if (TRUE) 
         symTab$addSymbol(code$args[[1]]$type)
-
-    ## Now the 3rd arg, the body of the loop, can be processed
+        
+    ## Now the 3rd arg, the body of the loop, can be processed.
+    ## For now, we will handle local vars in body as `copyVars` that become vars
+    ## in the encompassing method, but consider setting up local symbol table for
+    ## the loop body with the loop body C++ function declaring its own variables.
+    auxEnv$inParallelBody <- TRUE  # Needed for replacement of `self`.
+    symbolsNoBody <- symTab$getSymbolNames()
     inserts <- c(inserts, compile_labelAbstractTypes(code$args[[3]], symTab, auxEnv))
+    auxEnv$inParallelBody <- NULL
     
     auxEnv$uses_TBB <- TRUE
     nCompiler_pluginEnv$uses_TBB <- TRUE
+
+    symbols <- symTab$getSymbolNames()
+    code$aux$localVars <- symbols[!symbols %in% symbolsNoBody]
+
+    symbols <- symTab$getSymbolNames()
+    code$aux$localVars <- symbols[!symbols %in% symbolsNoBody]
+
     ## I think there shouldn't be any inserts returned since the body should be a bracket expression.
+    inserts <- c(inserts, compile_labelAbstractTypes(code$args[['nThreads']], symTab, auxEnv))
+    
     return(if (length(inserts) == 0) invisible(NULL) else inserts)
   }
 )
 
 inLabelAbstractTypesEnv(
   ParallelReduce <- function(code, symTab, auxEnv, handlingInfo) {
-    if (length(code$args) != 3)
+    inserts <- NULL
+    if(is.null(symTab$parentST))   #  TODO: this seems kludgey and perhaps should be done at a different processing stage.
       stop(exprClassProcessingErrorMsg(
         code,
-        paste('In labelAbstractTypes handler ParallelReduce:',
-              'expected 3 arguments but got', length(code$args))),
-        call. = FALSE)
-    ## process the initial value
-    inserts <- compile_labelAbstractTypes(code$args[[3]], symTab, auxEnv)
-    if (code$args[[3]]$type$nDim != 0)
+        paste0('In labelAbstractTypes handler ParallelReduce: ',
+               'parallel_reduce must be used in a method of an nClass, not in a stand-alone nFunction.')),
+        call. = FALSE)          
+    operatorDef <- operatorDefEnv[[code$args[['operator']]$name]]
+    if (code$args[['operator']]$name != '$' && !is.null(operatorDef) && is.null(operatorDef$reduction))   # Check for validity only for our operators.
+    # TODO: perhaps this should just be a warning.        
       stop(exprClassProcessingErrorMsg(
         code,
-        paste('In labelAbstractTypes handler ParallelReduce:',
-              'initial value for parallel_reduce should be scalar but got',
-              ' nDim = ', code$args[[3]]$type$nDim)),
+        paste0('In labelAbstractTypes handler ParallelReduce: ',
+               'function/operator `', code$args[['operator']]$name, '` is not a valid reduction function/operator.')),
         call. = FALSE)
-    if (isFALSE(code$args[[3]]$isLiteral))
+    if (length(code$args) < 3 || !'operator' %in% names(code$args) || !'object' %in% names(code$args))
       stop(exprClassProcessingErrorMsg(
         code,
         paste('In labelAbstractTypes handler ParallelReduce:',
-              'initial value for parallel_reduce must be a literal')),
+              'unexpected arguments -- at least two arguments (`operator` and `object`) are required.')),
         call. = FALSE)
-    ## process the reduce operator
-    if (isTRUE(code$args[[1]]$isLiteral)) {
-      if (!is.character(code$args[[1]]$name))
+    if(code$args[['operator']]$isName) {  ## Handle reduction function as function not char.        
+      code$args[['operator']]$isName <- FALSE
+      code$args[['operator']]$isLiteral <- TRUE
+      code$args[['operator']]$Rexpr <- deparse(code$args[['operator']]$Rexpr)
+    }
+    ## Process the reduce operator.
+    if (isTRUE(code$args[['operator']]$isLiteral)) {
+      if (!is.character(code$args[['operator']]$name))
         stop(exprClassProcessingErrorMsg(
           code,
           paste('In labelAbstractTypes handler ParallelReduce:',
                 'do not know how to use a reduce operator of type',
-                typeof(code$args[[1]]$name))),
+                typeof(code$args[['operator']]$name))),
           call. = FALSE)
-      code$args[[1]]$isLiteral <- FALSE
-      code$args[[1]]$isCall <- TRUE
+      code$args[['operator']]$isLiteral <- FALSE
+      code$args[['operator']]$isCall <- TRUE
+    }
+    if(code$args[['operator']]$name == "$") {
+      if(code$args[['operator']]$args[[1]]$name == "$")
+        stop(exprClassProcessingErrorMsg(
+          code,
+          paste('In labelAbstractTypes handler ParallelReduce:',
+                'too many levels of class hierarchy in reduction operator',
+                deparse(code$args[['operator']]$Rexpr))),
+          call. = FALSE)  
+      code$args[['operator']] <- wrapInExprClass(code$args[['operator']], 'chainedCall')
+      auxEnv$inParallelBody <- TRUE  # Needed for replacement of `self`.
+      inserts <- c(inserts, compile_labelAbstractTypes(code$args[['operator']], symTab, auxEnv))
+      auxEnv$inParallelBody <- NULL
     }
-    ## give reduce operator the same return type as the initial value
+    
+    ## Give reduce operator the same return type as the input vector.
     ## TODO: Maybe symbolNF is the right type for the reduction op.
-    code$args[[1]]$type <-
-      symbolBasic$new(name = code$args[[1]]$name,
-                      nDim = 0, type = code$args[[3]]$type$type)
-    ## finish by processing the vector arg
-    inserts <- c(inserts, compile_labelAbstractTypes(code$args[[2]], symTab, auxEnv))
-    if (code$args[[2]]$type$nDim != 1)
+    code$args[['operator']]$type <-
+      symbolBasic$new(name = code$args[['operator']]$name,
+                      nDim = 0, type = code$args[['object']]$type$type)
+    ## Process the vector arg.
+    ## TODO: we want to handle if vector is an expression (including obj$x),
+    ## presumably by lifting.
+    if(!code$args[['object']]$isName)
+      stop(exprClassProcessingErrorMsg(
+        code,
+        paste('In labelAbstractTypes handler ParallelReduce:',
+              'vector argument for parallel_reduce must be a variable, but found an expression `',
+              deparse(code$args[['object']]$Rexpr),
+              '`. Please create a temporary variable to use as the second argument.')),
+        call. = FALSE)
+            
+    inserts <- c(inserts, compile_labelAbstractTypes(code$args[['object']], symTab, auxEnv))
+    if (code$args[['object']]$type$nDim != 1)
       stop(exprClassProcessingErrorMsg(
         code,
         paste('In labelAbstractTypes handler ParallelReduce:',
-              'expected the second argument to be a vector but got nDim = ',
-              code$args[[2]]$type$nDim)),
+              'expected the `object` argument to be a vector but got nDim = ',
+              code$args[['object']]$type$nDim)),
         call. = FALSE)
     code$type <- symbolBasic$new(name = code$name, nDim = 0,
-                                 type = code$args[[3]]$type$type)
-    
+                                 type = code$args[['object']]$type$type)
+
+
+    ## Process the initial value.
+    if(!'init' %in% names(code$args) && !is.null(operatorDef$reduction))
+      setArg(code, 'init', nParse(operatorDef$reduction), add = TRUE)
+    if(!'init' %in% names(code$args))
+      stop(exprClassProcessingErrorMsg(
+        code,
+        paste('In labelAbstractTypes handler ParallelReduce:',
+              '`init` argument not provided and no default `init` is available for the operator')),
+        call. = FALSE)
+            
+    inserts <- c(inserts, compile_labelAbstractTypes(code$args[['init']], symTab, auxEnv))
+    if (code$args[['init']]$type$nDim != 0)
+      stop(exprClassProcessingErrorMsg(
+        code,
+        paste('In labelAbstractTypes handler ParallelReduce:',
+              'initial value for parallel_reduce should be scalar but got',
+              ' nDim = ', code$args[[3]]$type$nDim)),
+        call. = FALSE)
+    if (isFALSE(code$args[['init']]$isLiteral)) {
+      if(!(code$args[['init']]$name == "-" && isTRUE(code$args[['init']]$args[[1]]$isLiteral)))  # Handle negative init.
+        stop(exprClassProcessingErrorMsg(
+          code,
+          paste('In labelAbstractTypes handler ParallelReduce:',
+                'initial value for parallel_reduce must be a literal value, not a variable or expression')),
+          call. = FALSE)
+    }
+       
+
+    inserts <- c(inserts, compile_labelAbstractTypes(code$args[['nThreads']], symTab, auxEnv))
+
     auxEnv$uses_TBB <- TRUE
     nCompiler_pluginEnv$uses_TBB <- TRUE
+
     return(if (length(inserts) == 0) invisible(NULL) else inserts)
   }
 )
diff --git a/nCompiler/R/compile_normalizeCalls.R b/nCompiler/R/compile_normalizeCalls.R
index a4c77f63..11e27249 100644
--- a/nCompiler/R/compile_normalizeCalls.R
+++ b/nCompiler/R/compile_normalizeCalls.R
@@ -1,3 +1,9 @@
+## Special cases placed here by analogy with `eigenizeUseArgs`,
+## but perhaps should be in handler list.
+normalizeCallsFunctionArgs <- list(
+    parallel_reduce = 1
+)
+
 normalizeCallsEnv <- new.env()
 normalizeCallsEnv$.debug <- FALSE
 
@@ -24,9 +30,19 @@ compile_normalizeCalls <- function(code,
   logging <- get_nOption('compilerOptions')[['logging']]
   if (logging) appendToLog(paste('###', nErrorEnv$stateInfo, '###'))
 
-  if(code$isLiteral) return(NULL)
-  if(code$isName) return(NULL)
-  if(code$isCall) {
+  ## Handle arguments that are functions (`parallel_reduce`).
+  ## Do this by looking for parallel_reduce as the caller rather than as the call
+  ## so that other args are handled as usual.
+  fxnArg <- NULL
+  if(!is.null(code$caller)) {
+    fxnArg <- normalizeCallsFunctionArgs[[code$caller$name]]
+    if(!is.null(fxnArg) && fxnArg != code$callerArgID)
+      fxnArg <- NULL
+  }
+  
+  if(code$isLiteral && is.null(fxnArg)) return(NULL)
+  if(code$isName && is.null(fxnArg)) return(NULL)
+  if(code$isCall || !is.null(fxnArg)) {
     if(code$name == '{') {
       ## recurse over lines
       for(i in seq_along(code$args)) {
@@ -88,13 +104,15 @@ compile_normalizeCalls <- function(code,
       return(NULL)
     }
 
-    opDef <- cachedOpInfo$opDef
-    matchDef <- opDef[["matchDef"]]
-    if(is.null(matchDef))
-      matchDef <- cachedOpInfo$obj_internals$default_matchDef
-    if(!is.null(matchDef)) {
-      exprClass_put_args_in_order(matchDef, code, opDef$compileArgs)
-      # code <- replaceArgInCaller(code, matched_code)
+    if(is.null(fxnArg)) {
+      opDef <- cachedOpInfo$opDef
+      matchDef <- opDef[["matchDef"]]
+      if(is.null(matchDef))
+        matchDef <- cachedOpInfo$obj_internals$default_matchDef
+      if(!is.null(matchDef)) {
+        exprClass_put_args_in_order(matchDef, code, opDef$compileArgs)
+        # code <- replaceArgInCaller(code, matched_code)
+      }
     }
     normalizeCallsEnv$recurse_normalizeCalls(code, symTab, auxEnv, handlingInfo)
   }
@@ -152,7 +170,9 @@ update_cachedOpInfo <- function(code, where, allowFail=FALSE) {
         obj <- NC_find_method(where, code$name, inherits=TRUE)
         if(!is.null(obj)) {
           if(isNF(obj)) {
-            cachedOpInfo$case <- "nClass method" # possibly disambiguate method from keyword
+            if(!checkForLiftedBody(code)) {  # Current lifted body cases are `parallel_{for,reduce}`.
+              cachedOpInfo$case <- "nClass method" # possibly disambiguate method from keyword
+            } else cachedOpInfo$case <- "nClass method in lifted"  # a method call in a code block that will be lifted out of the class def and wil
             cachedOpInfo$obj_internals <- NFinternals(obj)
             opDef <- cachedOpInfo$obj_internals$compileInfo$opDef # might be NULL
           } else {
@@ -213,6 +233,10 @@ update_cachedOpInfo <- function(code, where, allowFail=FALSE) {
       if(cachedOpInfo$case == "nFunction" || cachedOpInfo$case == "nClass method") {
         opDef <- getOperatorDef("nFunction_default")
       }
+      if(cachedOpInfo$case == "nClass method in lifted") {
+        opDef <- getOperatorDef("nClass_method_in_lifted")
+      }
+      
     }
   }
   if(is.null(opDef)) {
@@ -234,3 +258,14 @@ update_cachedOpInfo <- function(code, where, allowFail=FALSE) {
   code$aux$cachedOpInfo <- cachedOpInfo
   cachedOpInfo
 }
+
+checkForLiftedBody <- function(code) {
+  while(!is.null(code$caller)) {  # Caller needs to be in specific set of operators and method in particular argument.
+    if(code$caller$isCall && code$caller$name %in% names(liftedBlockOperatorsArg))
+      if(code$callerArgID %in% liftedBlockOperatorsArg[[code$caller$name]])
+       return(TRUE) else return(FALSE)
+    return(checkForLiftedBody(code$caller))
+  }
+  return(FALSE)
+}
+    
diff --git a/nCompiler/R/cppDefs_R_interface_calls.R b/nCompiler/R/cppDefs_R_interface_calls.R
index 015123ab..b1792d4a 100644
--- a/nCompiler/R/cppDefs_R_interface_calls.R
+++ b/nCompiler/R/cppDefs_R_interface_calls.R
@@ -46,7 +46,7 @@ global_R_interface_cppDef <-
       "    get_genericInterfaceBaseC(Xptr);\n",
       "  //  std::cout << name << std::endl;\n",
       "  return(obj->call_method( name, Sargs ));\n",
-      "}\n"),
+      "};\n"),
     name = "R_interfaces"
   )
 
diff --git a/nCompiler/R/cppDefs_TBB.R b/nCompiler/R/cppDefs_TBB.R
index 7ec9d27c..7f04c715 100644
--- a/nCompiler/R/cppDefs_TBB.R
+++ b/nCompiler/R/cppDefs_TBB.R
@@ -1,6 +1,10 @@
-# not working
 ## cppDefs for parallel loop bodies for TBB
 
+selfNameInLiftedBlock <- "obj__"
+
+parallelForBodyLabelMaker <- labelFunctionCreator('parallel_loop_body')
+parallelReduceBodyLabelMaker <- labelFunctionCreator('parallel_reduce_body')   
+
 cppParallelBodyClass <- R6::R6Class(
   'cppParallelBodyClass',
   inherit = cppClassClass,
@@ -10,13 +14,15 @@ cppParallelBodyClass <- R6::R6Class(
                           loop_var,
                           symbolTable,
                           copyVars = character(),
-                          noncopyVars = character()) {
+                          noncopyVars = character(),
+                          aux = list()) {
       cppParallelBodyClass_init_impl(self,
                                      loop_body = loop_body,
                                      loop_var = loop_var,
                                      symbolTable = symbolTable,
                                      copyVars = copyVars,
-                                     noncopyVars = noncopyVars)
+                                     noncopyVars = noncopyVars,
+                                     aux = aux)
     },
     generate = function(declaration = FALSE, ...) {
       ## This version of generate creates a fully inlined version
@@ -47,13 +53,12 @@ cppParallelBodyClass <- R6::R6Class(
 )
 
 cppParallelBodyClass_init_impl <- function(cppDef,
-                                           name = "parallel_loop_body",
-                                           orig_loop_code = orig_loop_code,
-                                           loop_body = orig_loop_code$args[[3]],
-                                           loop_var = orig_loop_code$args[[1]],
+                                           loop_body,
+                                           loop_var,
                                            symbolTable,
                                            copyVars,
-                                           noncopyVars) {
+                                           noncopyVars,
+                                           aux) {
   ## 1. Create symbolTable for copyVars + noncopyVars
   ## 2. Create operator()
   ## 3. Create constructor
@@ -99,6 +104,11 @@ cppParallelBodyClass_init_impl <- function(cppDef,
     sym$ref <- TRUE
     newSymTab$addSymbol(sym)  
   }
+  if(length(aux$localMethods) || isTRUE(aux$liftedSelf)) 
+    newSymTab$addSymbol(cppVarFullClass$new(name = selfNameInLiftedBlock,
+                                            baseType = aux$class,
+                                            ref = TRUE))
+    
   ## Create operator()
   generalForExpr <- exprClass$new(name = 'GeneralFor', isCall = TRUE,
                                   isName = FALSE, isAssign = FALSE, isLiteral = FALSE)
@@ -134,7 +144,7 @@ cppParallelBodyClass_init_impl <- function(cppDef,
                                                  list(X = as.name(thisSymName),
                                                       X_ = as.name(thisArgName))))
   }
-  constructor <- cppFunctionClass$new(name = name,
+  constructor <- cppFunctionClass$new(name = aux$bodyName,
                                       args = ctorArgSymTab,
                                       code = cppCodeBlockClass$new(
                                         code = nParse(quote({})),
@@ -142,7 +152,7 @@ cppParallelBodyClass_init_impl <- function(cppDef,
                                       ),
                                       initializerList = initializerList,
                                       returnType = cppBlank())
-  cppDef$name <- name
+  cppDef$name <- aux$bodyName
   cppDef$memberCppDefs <- list(`operator()` = `operator()`,
                                  constructor = constructor)
   cppDef$symbolTable <- newSymTab
@@ -158,13 +168,17 @@ cppParallelReduceBodyClass <- R6::R6Class(
                           loop_var,
                           symbolTable,
                           copyVars = character(),
-                          noncopyVars = character()) {
+                          noncopyVars = character(),
+                          aux = list(),
+                          lastUse = FALSE) {
       cppParallelReduceBodyClass_init_impl(self,
                                            loop_body = loop_body,
                                            loop_var = loop_var,
                                            symbolTable = symbolTable,
                                            copyVars = copyVars,
-                                           noncopyVars = noncopyVars)
+                                           noncopyVars = noncopyVars,
+                                           aux = aux,
+                                           lastUse = lastUse)
     },
     generate = function(declaration = FALSE, ...) {
       ## This version of generate creates a fully inlined version
@@ -190,20 +204,18 @@ cppParallelReduceBodyClass <- R6::R6Class(
         unlist(output)
       } else
         ""
-      ## TODO: C++ generation of the original nFunction where parallel_reduce appears gets
-      ## 'int i__;' and 'double value__;', seemingly because of symbol table sharing.
     }
   )
 )
 
 cppParallelReduceBodyClass_init_impl <- function(cppDef,
-                                                 name = "parallel_reduce_body",
-                                                 orig_loop_code = orig_loop_code,
-                                                 loop_body = orig_loop_code$args[[3]],
-                                                 loop_var = orig_loop_code$args[[1]],
+                                                 loop_body,
+                                                 loop_var,
                                                  symbolTable,
                                                  copyVars,
-                                                 noncopyVars) {
+                                                 noncopyVars,
+                                                 aux,
+                                                 lastUse) {
   ## 1. call cppParallelBodyClass_init_impl which creates GeneralFor
   ## 2. make some minor alterations to the body of `operator()`
   ## 3. Create split constructor
@@ -211,10 +223,10 @@ cppParallelReduceBodyClass_init_impl <- function(cppDef,
 
   ## need to save this here because cppParallelBodyClass_init_impl will change
   ## loop_body's caller
-  orig_caller <- loop_body$caller
+  orig_caller <- copyExprClass(loop_body$caller)
 
-  cppParallelBodyClass_init_impl(cppDef, name, orig_loop_code, loop_body,
-                                 loop_var, symbolTable, copyVars, noncopyVars)
+  cppParallelBodyClass_init_impl(cppDef, loop_body, loop_var,
+                                 symbolTable, copyVars, noncopyVars, aux)
 
   ## get the local aggregation var copy variable
   val_expr <- copyExprClass(loop_body$args[[1]])
@@ -242,13 +254,9 @@ cppParallelReduceBodyClass_init_impl <- function(cppDef,
     cppVarClass$new(name = val_expr$name, baseType = val_expr$type$type))
   ## remove 'const' from the `operator()` declaration
   cppDef$memberCppDefs[['operator()']]$const <- FALSE
-
-  ## get the reduce op's identity element which is guaranteed to be a literal
-  ## by the labelAbstractTypes ParallelReduce handler
-  init_arg <- copyExprClass(orig_caller$caller$caller$args[[1]]$args[[2]])
-
+    
   split_ctor_symTab <- symbolTableClass$new()
-  split_ctor_symTab$addSymbol(cppVarClass$new(name = 'parent', baseType = name,
+  split_ctor_symTab$addSymbol(cppVarClass$new(name = 'parent', baseType = aux$bodyName,
                                               ref = TRUE))
   split_ctor_symTab$addSymbol(cppVarClass$new(name = 'tbb::split'))
   ## Get the name of the vector we're working with, which together with the
@@ -256,13 +264,24 @@ cppParallelReduceBodyClass_init_impl <- function(cppDef,
   vector_name <- orig_caller$args[[4]] ## should be a string
   initializerList <- list()
   initializerList[[1]] <- nParse(
-    substitute(X(X_), list(X = as.name(value_name),
-                           X_ = as.name(init_arg$name))))
+    substitute(X(X_), list(X = as.name(value_name))))
+  ## Need to directly parse the init value to handle various numeric cases, e.g., `Inf`.
+  setArg(initializerList[[1]], 1, orig_caller$aux$init) 
   initializerList[[2]] <- nParse(
     substitute(X(X_), list(X = as.name(vector_name),
                            X_ = as.name(paste0('parent.', vector_name)))))
+  if(length(aux$localMethods) || isTRUE(aux$liftedSelf)) 
+    initializerList[[3]] <- nParse(
+      substitute(X(X_), list(X = selfNameInLiftedBlock, 
+                             X_ = as.name(paste0('parent.', selfNameInLiftedBlock)))))
+  if(length(orig_caller$args) == 6 && length(orig_caller$args[[6]]))  # This is the object if using an object method as the operator.
+    initializerList[[3]] <- nParse(
+      substitute(X(X_), list(X = orig_caller$args[[6]], 
+                             X_ = as.name(paste0('parent.', orig_caller$args[[6]])))))
+        
 
-  split_constructor <- cppFunctionClass$new(name = name,
+    
+  split_constructor <- cppFunctionClass$new(name = aux$bodyName,
                                             args = split_ctor_symTab,
                                             code = cppCodeBlockClass$new(
                                               code = nParse(quote({})),
@@ -274,18 +293,28 @@ cppParallelReduceBodyClass_init_impl <- function(cppDef,
   ## join_symTab is the symbolTable for the arguments to join
   join_symTab <- symbolTableClass$new()
   join_symTab$addSymbol(cppVarFullClass$new(name = 'target',
-                                            baseType = name,
+                                            baseType = aux$bodyName,
                                             ref = TRUE,
                                             const = TRUE))
+
   ## make the reduce code
-  reduce_op <- exprClass$new(name = loop_body$args[[2]]$name, isCall = TRUE,
-                             isName = FALSE, isAssign = FALSE,
+  ## `aux` needed so that user-defined reduction functions will be replaced with `cpp_code_name`.
+  if(loop_body$args[[2]]$name == 'chainedCall') {
+    reduce_op <- copyExprClass(loop_body$args[[2]]) 
+    inc <- 1
+  } else {
+    reduce_op <- exprClass$new(name = loop_body$args[[2]]$name, aux = loop_body$args[[2]]$aux,
+                             isCall = TRUE, isName = FALSE, isAssign = FALSE,
                              isLiteral = FALSE)
-  setArg(reduce_op, 1, copyExprClass(value_expr))
-  setArg(reduce_op, 2, nParse(paste0('cppLiteral("target.', value_name, ';")')))
+    inc <- 0
+  }
+  setArg(reduce_op, 1+inc, copyExprClass(value_expr))
+  setArg(reduce_op, 2+inc, nParse(paste0('cppLiteral("target.', value_name, '")')))
   join_code <- newAssignmentExpression()
   setArg(join_code, 1, copyExprClass(value_expr))
   setArg(join_code, 2, reduce_op)
+  ## Put code in {} so handled by full processing later, in particular adding ending `;`.  
+  join_code <- newBracketExpr(list(join_code))  
   ## create the join cppFunctionClass definition
   join_body <- cppCodeBlockClass$new(code = join_code,
                                      ## TODO: any symbols ever needed?
@@ -299,5 +328,13 @@ cppParallelReduceBodyClass_init_impl <- function(cppDef,
   cppDef$memberCppDefs <- c(cppDef$memberCppDefs,
                               list(split_constructor = split_constructor,
                                    join = join))
+
+  ## Remove index and value variables so not defined in calling method.
+  if(lastUse) {
+    symbolTable$removeSymbol(loop_var$name)
+    symbolTable$removeSymbol(value_name)
+  }
   invisible(NULL)
 }
+
+
diff --git a/nCompiler/R/cppDefs_nClass.R b/nCompiler/R/cppDefs_nClass.R
index 6e0d6087..39b94cbc 100644
--- a/nCompiler/R/cppDefs_nClass.R
+++ b/nCompiler/R/cppDefs_nClass.R
@@ -258,23 +258,18 @@ cpp_nClassClass <- R6::R6Class(
     buildParallelClassDefs = function() {
       for(i in seq_along(Compiler$NFcompilers)) {
         parallelContent <- Compiler$NFcompilers[[i]]$auxEnv$parallelContent
-        if(!is.null(parallelContent)) {
-          for(j in seq_along(parallelContent)) {
-            cppDef_TBB <- cppParallelBodyClass$new(loop_body = parallelContent[[j]]$args[[3]],
+        if(length(parallelContent)) {
+            for(j in seq_along(parallelContent)) {
+             cppDef_TBB <- cppParallelBodyClass$new(loop_body = parallelContent[[j]]$args[[3]],
                                                    loop_var = parallelContent[[j]]$args[[1]],
                                                    symbolTable = memberCppDefs[[i]]$code$symbolTable,
-                                                   copyVars = parallelContent[[j]]$args[[4]],
-                                                   noncopyVars = parallelContent[[j]]$args[[5]])
-            ## The name is hard-wired expecting only a single case of parallel content.
-            ## TO-DO: generalize the name with unique identifier.
-            self$memberCppDefs[["parallel_loop_body"]] <<- cppDef_TBB
+                                                   copyVars = parallelContent[[j]]$args[['copyVars']],
+                                                   noncopyVars = parallelContent[[j]]$args[['shareVars']],
+                                                   aux = parallelContent[[j]]$aux)
+            self$memberCppDefs[[parallelContent[[j]]$aux$bodyName]] <<- cppDef_TBB
           }
         }
         parallelReduceContent <- Compiler$NFcompilers[[i]]$auxEnv$parallelReduceContent
-        ## TODO: if there are multiple nFunctions in a class, they share the
-        ## auxEnv and as a result we'll find the parallelReduceContent even for
-        ## the wrong NFcompiler (and thus pass the wrong symbolTable to
-        ## cppParallelReduceBodyClass$new
         if(!is.null(parallelReduceContent)) {
           for(j in seq_along(parallelReduceContent)) {
             cppDef_TBB <- cppParallelReduceBodyClass$new(
@@ -282,12 +277,12 @@ cpp_nClassClass <- R6::R6Class(
               loop_var = parallelReduceContent[[j]]$args[[1]],
               symbolTable = memberCppDefs[[i]]$code$symbolTable,
               copyVars = list(),
-              noncopyVars = list(parallelReduceContent[[j]]$args[[4]],
-                                 parallelReduceContent[[j]]$args[[5]])
-            )
-            ## The name is hard-wired expecting only a single case of parallel content.
-            ## TO-DO: generalize the name with unique identifier.
-            self$memberCppDefs[["parallel_reduce_body"]] <<- cppDef_TBB
+              noncopyVars = as.list(c(parallelReduceContent[[j]]$args[['input']],
+                                      parallelReduceContent[[j]]$args[['output']],
+                                      parallelReduceContent[[j]]$args[['nClass_object']])),              
+              aux = parallelReduceContent[[j]]$aux,
+              lastUse = (j == length(parallelReduceContent)))
+            self$memberCppDefs[[parallelReduceContent[[j]]$aux$bodyName]] <<- cppDef_TBB
           }
         }
       }
diff --git a/nCompiler/R/options.R b/nCompiler/R/options.R
index 6c9c8885..0a575d18 100644
--- a/nCompiler/R/options.R
+++ b/nCompiler/R/options.R
@@ -52,7 +52,8 @@ updateDefaults <- function(defaults, control) {
     sourceCpp_verbose = FALSE,
     nimble = FALSE, ## ensure all backward compatibility
     dropSingleSizes = FALSE, ## backward compatibility
-    useSafeDeparse = TRUE
+    useSafeDeparse = TRUE,
+    nThreads = 0  ## 0 corresponds to TBB using all available threads
   )
 )
 
diff --git a/nCompiler/R/symbolTable.R b/nCompiler/R/symbolTable.R
index a8bf1b86..b43bf1ee 100644
--- a/nCompiler/R/symbolTable.R
+++ b/nCompiler/R/symbolTable.R
@@ -383,7 +383,7 @@ symbolSelf <- R6::R6Class(
   inherit = symbolNC,
   portable = TRUE,
   public = list(
-    initialize = function(name, type, NCgenerator, isArg = FALSE) {
+      initialize = function(name, type, NCgenerator, isArg = FALSE) {
       super$initialize(name = name,
                        type = type,
                        NCgenerator = NCgenerator,
@@ -395,7 +395,7 @@ symbolSelf <- R6::R6Class(
     # The cppVar for self is only used if self is used as a value, such as in an argument or return value.
     genCppVar = function() cppVarSelfClass$new()
   )
-)
+  )
 
 ## type is the unique ID of the NCgenerator.
 ## same value as for a symbolNC for an object of the class.
diff --git a/nCompiler/inst/include/.DS_Store b/nCompiler/inst/include/.DS_Store
deleted file mode 100644
index d96ea5ae..00000000
Binary files a/nCompiler/inst/include/.DS_Store and /dev/null differ
diff --git a/nCompiler/inst/include/nCompiler/.DS_Store b/nCompiler/inst/include/nCompiler/.DS_Store
deleted file mode 100644
index 3e8e0d9e..00000000
Binary files a/nCompiler/inst/include/nCompiler/.DS_Store and /dev/null differ
diff --git a/nCompiler/inst/include/nCompiler/ET_Rcpp_ext/.DS_Store b/nCompiler/inst/include/nCompiler/ET_Rcpp_ext/.DS_Store
deleted file mode 100644
index 5008ddfc..00000000
Binary files a/nCompiler/inst/include/nCompiler/ET_Rcpp_ext/.DS_Store and /dev/null differ
diff --git a/nCompiler/inst/include/nCompiler/Rcpp_extensions/.DS_Store b/nCompiler/inst/include/nCompiler/Rcpp_extensions/.DS_Store
deleted file mode 100644
index 5008ddfc..00000000
Binary files a/nCompiler/inst/include/nCompiler/Rcpp_extensions/.DS_Store and /dev/null differ
diff --git a/nCompiler/inst/include/nCompiler/nC_inter/.DS_Store b/nCompiler/inst/include/nCompiler/nC_inter/.DS_Store
deleted file mode 100644
index 5008ddfc..00000000
Binary files a/nCompiler/inst/include/nCompiler/nC_inter/.DS_Store and /dev/null differ
diff --git a/nCompiler/inst/include/nCompiler/nCompiler_omnibus.h b/nCompiler/inst/include/nCompiler/nCompiler_omnibus.h
index a565860c..5bf16a92 100644
--- a/nCompiler/inst/include/nCompiler/nCompiler_omnibus.h
+++ b/nCompiler/inst/include/nCompiler/nCompiler_omnibus.h
@@ -23,6 +23,23 @@
 //#ifndef NCOMPILER_OMNIBUS_FIRST_CPP_
 // #define NCOMPILER_OMNIBUS_FIRST_CPP_
 
+// Despite use of `inline` and C++17, TBB_DEPTH is being defined multiple times.
+// Also, I think I need the include guard for the class definition anyway.
+// The header guard is needed for the class definition and to
+// avoid multiple definitions of `TBB_DEPTH` (despite the `inline`; I
+// don't fully understand that).
+#ifndef TBB_DEPTH_H
+#define TBB_DEPTH_H
+
+inline int TBB_DEPTH = 0;
+
+class depth_counter {
+  public:
+    depth_counter() { TBB_DEPTH++; }
+    ~depth_counter() { TBB_DEPTH--; }
+};
+#endif
+   
 #include "nCompiler_omnibus_pre_Rcpp.h" // should always be redundant, but it is here to be clear
 
 // We shall see if these two "first" files (_h and _cpp) are really needed or can be consolidated.
diff --git a/nCompiler/inst/include/nCompiler/predef/.DS_Store b/nCompiler/inst/include/nCompiler/predef/.DS_Store
deleted file mode 100644
index fc6fdc43..00000000
Binary files a/nCompiler/inst/include/nCompiler/predef/.DS_Store and /dev/null differ
diff --git a/nCompiler/inst/include/nCompiler/utils.h b/nCompiler/inst/include/nCompiler/utils.h
index a3fc4c09..e2b6e543 100644
--- a/nCompiler/inst/include/nCompiler/utils.h
+++ b/nCompiler/inst/include/nCompiler/utils.h
@@ -1,3 +1,18 @@
+inline int getNumThreads(double value_) {
+  static int nThreads_nOption=0;   // This variable is preserved across calls.
+  if(TBB_DEPTH==1) {   // Avoid calls to R API within threads.
+     Rcpp::Environment nc = Rcpp::Environment::namespace_env("nCompiler");
+     Rcpp::Function get_nOption = nc["get_nOption"];
+     nThreads_nOption=Rcpp::as<int>(get_nOption("nThreads"));
+  }
+  int value = (int) value_;
+  if (nThreads_nOption > 0)
+    value = nThreads_nOption;
+  if(value == 0)
+    value = 100000;
+  return value;
+}
+
 #ifndef _NC_UTILS_
 #define _NC_UTILS_
 
diff --git a/nCompiler/tests/.DS_Store b/nCompiler/tests/.DS_Store
deleted file mode 100644
index a1f14fc9..00000000
Binary files a/nCompiler/tests/.DS_Store and /dev/null differ
diff --git a/nCompiler/tests/testthat/.DS_Store b/nCompiler/tests/testthat/.DS_Store
deleted file mode 100644
index 5008ddfc..00000000
Binary files a/nCompiler/tests/testthat/.DS_Store and /dev/null differ
diff --git a/nCompiler/tests/testthat/nCompile_tests/.DS_Store b/nCompiler/tests/testthat/nCompile_tests/.DS_Store
deleted file mode 100644
index 5008ddfc..00000000
Binary files a/nCompiler/tests/testthat/nCompile_tests/.DS_Store and /dev/null differ
diff --git a/nCompiler/tests/testthat/tbb_tests/test-parallel_for.R b/nCompiler/tests/testthat/tbb_tests/test-parallel_for.R
new file mode 100644
index 00000000..55d7ff6e
--- /dev/null
+++ b/nCompiler/tests/testthat/tbb_tests/test-parallel_for.R
@@ -0,0 +1,495 @@
+test_that("basic usage of parallel_for", {
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- x
+                    mult <- 2
+                    parallel_for(i, 1:length(x), {y[i] <- mult*x[i]})
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+    obj <- nc$new()
+    Cobj <- Cnc$new()
+    expect_identical(obj$go(2:6), as.numeric(2*(2:6)))
+    expect_identical(Cobj$go(2:6), as.numeric(2*(2:6)))
+    
+    nc <- nClass(
+        Cpublic = list(
+            myconst = 'numericScalar',
+            twice = nFunction(
+                fun=function(x = 'numericScalar') {
+                    return(2*x)
+                }, returnType = 'numericScalar'
+            ),
+            thrice = nFunction(
+                fun=function(x = 'numericScalar') {
+                    return(3*x)
+                }, returnType = 'numericScalar'
+            ),
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- x
+                    parallel_for(i, 1:length(x), {y[i] <- myconst + twice(thrice(x[i]))})
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+    obj <- nc$new()
+    obj$myconst <- 7
+    Cobj <- Cnc$new()
+    Cobj$myconst <- 7
+    expect_identical(obj$go(2:6), as.numeric(6*(2:6)+7))
+    expect_identical(Cobj$go(2:6), as.numeric(6*(2:6)+7))
+
+    mult = nFunction(
+        fun=function(x = 'numericScalar', c = 'numericScalar') {
+            return(c*x)
+        }, returnType = 'numericScalar'
+    )
+    
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector'){
+                    y <- x
+                    parallel_for(i, 1:length(x), {y[i] <- mult(x[i], 3)})
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc,mult)[[1]]
+    obj <- nc$new()
+    Cobj <- Cnc$new()
+    expect_identical(obj$go(2:6), as.numeric(3*(2:6)))
+    expect_identical(Cobj$go(2:6), as.numeric(3*(2:6)))
+    
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- x
+                    mult <- 2
+                    parallel_for(i, 2:length(x), {y[i] <- mult*x[i]})
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+    obj <- nc$new()
+    Cobj <- Cnc$new()
+    expect_identical(obj$go(2:6), as.numeric(c(2, 2*(3:6))))
+    expect_identical(Cobj$go(2:6), as.numeric(c(2, 2*(3:6))))
+ })
+
+test_that("use of `self`", {
+  nc <- nClass(
+        Cpublic = list(
+            myconst = 'numericScalar',
+            twice = nFunction(
+                fun=function(x = 'numericScalar') {
+                    return(2*x)
+                }, returnType = 'numericScalar'
+            ),
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- x
+                    parallel_for(i, 1:length(x), {y[i] <- self$myconst + twice(x[i])})
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+  )
+  Cnc <- nCompile(nc)
+  obj <- nc$new()
+  obj$myconst <- 1
+  Cobj <- Cnc$new()
+  Cobj$myconst <- 1
+  expect_identical(obj$go(2:6), as.numeric(1 + 2*2:6)) 
+  expect_identical(Cobj$go(2:6), as.numeric(1 + 2*2:6)) 
+
+  nc <- nClass(
+        Cpublic = list(
+            myconst = 'numericScalar',
+            twice = nFunction(
+                fun=function(x = 'numericScalar') {
+                    return(2*x)
+                }, returnType = 'numericScalar'
+            ),
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- x
+                    parallel_for(i, 1:length(x), {y[i] <- myconst + self$twice(x[i])})
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+  )
+  Cnc <- nCompile(nc)
+  obj <- nc$new()
+  obj$myconst <- 1
+  Cobj <- Cnc$new()
+  Cobj$myconst <- 1
+  expect_identical(obj$go(2:6), as.numeric(1 + 2*2:6))
+  expect_identical(Cobj$go(2:6), as.numeric(1 + 2*2:6))
+
+  nc <- nClass(
+        Cpublic = list(
+            myconst = 'numericScalar',
+            twice = nFunction(
+                fun=function(x = 'numericScalar') {
+                    return(2*x)
+                }, returnType = 'numericScalar'
+            ),
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- x
+                    parallel_for(i, 1:length(x), {y[i] <- self$twice(self$myconst + x[i])})
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+  )
+  Cnc <- nCompile(nc)
+  Cnc <- nCompile(nc)
+  obj <- nc$new()
+  obj$myconst <- 1
+  Cobj <- Cnc$new()
+  Cobj$myconst <- 1
+  expect_identical(obj$go(2:6), as.numeric(2*3:7))
+  expect_identical(Cobj$go(2:6), as.numeric(2*3:7))
+
+  nc <- nClass(
+        Cpublic = list(
+            myconst = 'numericScalar',
+            myotherconst = 'numericScalar',
+            twice = nFunction(
+                fun=function(x = 'numericScalar') {
+                    return(2*x)
+                }, returnType = 'numericScalar'
+            ),
+            thrice = nFunction(
+                fun=function(x = 'numericScalar') {
+                    return(3*x)
+                }, returnType = 'numericScalar'
+            ),
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                  y <- x
+                  tmp <- self$myotherconst + 3   # Use `self` outside of parallel body.
+                  parallel_for(i, 1:length(x), {y[i] <- tmp + self$myconst + myotherconst + self$twice(thrice(x[i]))})
+                  return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    
+    Cnc <- nCompile(nc)
+    obj <- nc$new()
+    obj$myconst <- 7
+    obj$myotherconst <- 1
+    Cobj <- Cnc$new()
+    Cobj$myconst <- 7
+    Cobj$myotherconst <- 1
+    expect_identical(obj$go(2:6), as.numeric(6*(2:6)+12))
+    expect_identical(Cobj$go(2:6), as.numeric(6*(2:6)+12))
+})
+
+test_that("use of object from another class", {
+    nc0 <- nClass(
+        Cpublic = list(
+            foo = 'numericScalar',
+            twice = nFunction(
+                fun=function(x = 'numericScalar') {
+                    return(2*x)
+                }, returnType = 'numericScalar'
+            )
+        )
+    )
+    
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector', o = 'nc0') {
+                    y <- x
+                    parallel_for(i, 1:length(x),  
+                    {y[i] <- o$twice( x[i]) + o$foo },
+                    )
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc,nc0)
+
+    tmp1 <- nc0$new()
+    tmp1$foo <- 5
+    nc1 <- nc$new()
+    expect_identical(nc1$go(1:5, tmp1), as.numeric(2*(1:5)+5))
+    
+    Ctmp1 <- Cnc[[2]]$new()
+    Ctmp1$foo <- 5
+    Cnc1 <- Cnc[[1]]$new()
+    expect_identical(Cnc1$go(1:5, Ctmp1), as.numeric(2*(1:5)+5))
+
+    nc <- nClass(
+        Cpublic = list(
+            o = 'nc0',
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- x
+                    parallel_for(i, 1:length(x),  
+                    {y[i] <- o$twice( x[i]) + o$foo },
+                    )
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc,nc0)
+
+    nc1 <- nc$new()
+    o <- nc0$new()
+    o$foo <- 5
+    nc1$o <- o
+    expect_identical(nc1$go(1:5), as.numeric(2*(1:5)+5))
+    Cnc1 <- Cnc[[1]]$new()
+    Co <- Cnc[[2]]$new()
+    Co$foo <- 5
+    Cnc1$o <- Co
+    expect_identical(Cnc1$go(1:5), as.numeric(2*(1:5)+5))
+
+})
+
+test_that("specifying {copy,share}Vars", {
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- x
+                    parallel_for(i, 1:length(x),  
+                    {y[i] <- 2*x[i]; x[i] <- 0},
+                    copyVars = 'x',
+                    )
+                    return(x)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+    Cnc1 <- Cnc$new()
+    expect_identical(Cnc1$go(1:5), as.numeric(1:5))
+    
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- x
+                    parallel_for(i, 1:length(x),  
+                    {y[i] <- 2*x[i]; x[i] <- 0},
+                    shareVars = 'x'  # Also the default, so this is not really needed.
+                    )
+                    return(x)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+    Cnc1 <- Cnc$new()
+    expect_identical(Cnc1$go(1:5), as.numeric(rep(0,5)))
+    
+    
+    nc <- nClass(
+        Cpublic = list(
+            x = 'numericVector',
+            go = nFunction(
+                fun = function() {
+                    y <- x
+                    parallel_for(i, 1:length(x),  
+                    {y[i] <- 2*x[i]; x[i] <- 0}  # By default 'x' is shared (and therefore modified).
+                    )
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+    Cnc1 <- Cnc$new()
+    Cnc1$x <- 1:5
+    Cnc1$go()
+    expect_identical(Cnc1$x, as.numeric(rep(0,5)))
+
+    nc <- nClass(
+        Cpublic = list(
+            x = 'numericVector',
+            go = nFunction(
+                fun = function() {
+                    y <- x
+                    parallel_for(i, 1:length(x),  
+                    {y[i] <- 2*x[i]; x[i] <- 0},
+                    copyVars = 'x'
+                    )
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+    Cnc1 <- Cnc$new()
+    Cnc1$x <- 1:5
+    Cnc1$go()
+    expect_identical(Cnc1$x, as.numeric(1:5))
+    
+})
+
+
+test_that("lookup precedence", {
+    twice = nFunction(
+        fun=function(x = 'numericScalar') {
+            return(3*x)
+        }, returnType = 'numericScalar'
+    )
+    
+    nc <- nClass(
+        Cpublic = list(
+            twice = nFunction(
+                fun=function(x = 'numericScalar') {
+                    return(2*x)
+                }, returnType = 'numericScalar'
+            ),
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- x
+                    parallel_for(i, 1:length(x), {y[i] <- twice(x[i])})
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc,twice)[[1]]
+    obj <- nc$new()
+    Cobj <- Cnc$new()
+    expect_identical(obj$go(2:6), as.numeric(2*(2:6)))
+    expect_identical(Cobj$go(2:6), as.numeric(2*(2:6)))
+    
+})
+    
+
+
+test_that("multiple non-nested loops", {
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- x
+                    parallel_for(j, 1:length(x), {y[j] <- 2*x[j]})
+                    parallel_for(j, 1:length(x), {y[j] <- 3*y[j]})
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+    obj <- nc$new()
+    Cobj <- Cnc$new()
+    expect_identical(obj$go(2:6), as.numeric(6*(2:6)))
+    expect_identical(Cobj$go(2:6), as.numeric(6*(2:6)))
+
+    twice = nFunction(
+        fun=function(x = 'numericScalar') {
+            return(2*x)
+        }, returnType = 'numericScalar'
+    )
+    
+    nc <- nClass(
+        Cpublic = list(
+            thrice = nFunction(
+                fun=function(x = 'numericScalar') {
+                    return(3*x)
+                }, returnType = 'numericScalar'
+            ),
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- x
+                    parallel_for(i, 1:length(x),
+                    {y[i] <- twice(x[i])})
+                    parallel_for(i, 1:length(y),
+                    {y[i] <- thrice(y[i])})
+                    return(y)
+                },
+                returnType = 'numericVector'
+            ),
+            go2 = nFunction(
+                fun = function(x = 'numericVector',z='numericVector') {
+                    y <- x
+                    parallel_for(i, 1:length(x),
+                    {y[i] <- twice(x[i])+z[i]+3})
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    
+    Cnc <- nCompile(nc, twice)[[1]]
+    obj <- nc$new()
+    Cobj <- Cnc$new()
+    expect_identical(obj$go(2:6), as.numeric(6*(2:6)))
+    expect_identical(Cobj$go(2:6), as.numeric(6*(2:6)))
+    expect_identical(obj$go2(2:6, 100:104), as.numeric(2*(2:6)+100:104+3))
+    expect_identical(Cobj$go2(2:6,100:104), as.numeric(2*(2:6)+100:104+3))
+})
+ 
+test_that("nested loops", {   # See issue 152.
+
+    nc <- nClass(
+        Cpublic = list(
+            thrice = nFunction(
+                fun=function(x = 'numericVector') {
+                    y <- x
+                    parallel_for(i, 1:length(x), {y[i] <- 3*x[i]})
+                    return(y)
+                }, returnType = 'numericVector'
+            ),
+            go = nFunction(
+                fun = function(x = 'numericMatrix') {
+                    y <- x
+                    nc <- dim(x)[2]
+                    nr <- dim(x)[1]
+                    parallel_for(i, 1:nr, {y[i,1:nc] <- thrice(x[i,1:nc])})
+                    return(y)
+                },
+                returnType = 'numericMatrix'
+            )
+        ))
+
+    Cnc <- nCompile(nc)
+    obj <- nc$new()
+    Cobj <- Cnc$new()
+    input <- matrix(as.numeric(1:6), nrow=2)
+    expect_identical(obj$go(input), 3*input)
+    expect_identical(Cobj$go(input), 3*input)
+})
diff --git a/nCompiler/tests/testthat/tbb_tests/test-parallel_reduce.R b/nCompiler/tests/testthat/tbb_tests/test-parallel_reduce.R
new file mode 100644
index 00000000..84d90561
--- /dev/null
+++ b/nCompiler/tests/testthat/tbb_tests/test-parallel_reduce.R
@@ -0,0 +1,538 @@
+test_that("basic usage of parallel_reduce", {
+  nc <- nClass(
+    Cpublic = list(
+      go = nFunction(
+        fun = function(x = 'numericVector') {
+          y <- parallel_reduce('+', x)
+          return(y)
+        },
+        returnType = 'numericScalar'
+      )
+    )
+  )
+  Cnc <- nCompile(nc)
+  obj <- nc$new()
+  Cobj <- Cnc$new()
+  expect_identical(obj$go(101:110), as.numeric(sum(101:110)))
+  expect_identical(Cobj$go(101:110), as.numeric(sum(101:110)))
+
+  nc <- nClass(
+    Cpublic = list(
+      go = nFunction(
+        fun = function(x = 'numericVector') {
+          y <- parallel_reduce('+', x, 5)
+          return(y)
+        },
+        returnType = 'numericScalar'
+      )
+    )
+  )
+  Cnc <- nCompile(nc)
+  obj <- nc$new()
+  Cobj <- Cnc$new()
+  expect_identical(obj$go(101:110), as.numeric(5+sum(101:110)))
+  expect_identical(Cobj$go(101:110), as.numeric(5+sum(101:110)))
+
+  ## Negative values required some additional processing, so test that case explicitly.
+  nc <- nClass(
+    Cpublic = list(
+      go = nFunction(
+        fun = function(x = 'numericVector') {
+          y <- parallel_reduce('+', x, -5)
+          return(y)
+        },
+        returnType = 'numericScalar'
+      )
+    )
+  )
+  Cnc <- nCompile(nc)
+  obj <- nc$new()
+  Cobj <- Cnc$new()
+  expect_identical(obj$go(101:110), as.numeric(sum(101:110)-5))
+  expect_identical(Cobj$go(101:110), as.numeric(sum(101:110)-5))
+
+  nc <- nClass(
+    Cpublic = list(
+      go = nFunction(
+        fun = function(x = 'numericVector') {
+          y <- parallel_reduce('pairmin', x)
+          return(y)
+        },
+        returnType = 'numericScalar'
+      )
+    )
+  )
+  Cnc <- nCompile(nc)
+  obj <- nc$new()
+  Cobj <- Cnc$new()
+  x <- c(3.7, 2.5, 4.9, 3.1)
+  expect_identical(obj$go(x), 2.5)
+  expect_identical(Cobj$go(x), 2.5)
+
+
+  ## With additional code.
+  nc <- nClass(
+      Cpublic = list(
+          go = nFunction(
+              fun = function(x = 'numericVector') {
+                  tmp <- 7    
+                  y <- 3 + exp(parallel_reduce('+', x, 0))
+                  return(y)
+              },
+              returnType = 'numericScalar'
+          )
+      )
+  )
+  Cnc <- nCompile(nc)
+  obj <- nc$new()
+  Cobj <- Cnc$new()
+  expect_identical(obj$go(1:3), exp(6)+3)
+  expect_identical(Cobj$go(1:3), exp(6)+3)
+
+  ## Use in `return()`.
+  nc <- nClass(
+    Cpublic = list(
+        go = nFunction(
+        fun = function(x = 'numericVector') {
+          return(3+parallel_reduce('+', x, 0))
+        },
+        returnType = 'numericScalar'
+      )
+    )
+  )
+  Cnc <- nCompile(nc)
+  obj <- nc$new()
+  Cobj <- Cnc$new()
+  expect_identical(Cobj$go(1:3), 9)
+  expect_identical(obj$go(1:3), 9)
+})
+
+test_that("error trapping for parallel_reduce", {
+  nc <- nClass(
+    Cpublic = list(
+      go = nFunction(
+        fun = function(x = 'numericVector') {
+          y <- parallel_reduce('-', x)
+          return(y)
+        },
+        returnType = 'numericScalar'
+      )
+    )
+  )
+  ## The error message is not silent. 
+  expect_error(Cnc <- nCompile(nc), "is not a valid reduction")  ## Compile-time error.
+  obj <- nc$new()
+  expect_error(obj$go(1:5), "not a valid reduction")  ## Run-time error.
+
+  ## No init for user-defined reduction function.
+  mypairmin <- nFunction(
+      fun = function(x = 'numericScalar', y = 'numericScalar') {
+          return(pmin(x,y))
+      }, returnType = 'numericScalar'
+  )
+  
+  nc <- nClass(
+    Cpublic = list(
+      go = nFunction(
+        fun = function(x = 'numericVector') {
+          y <- parallel_reduce(mypairmin, x)
+          return(y)
+        },
+        returnType = 'numericScalar'
+      )
+    )
+  )
+  expect_error(Cnc <- nCompile(nc, mypairmin), "no default `init`")
+  obj <- nc$new()
+  expect_error(obj$go(1:5), "no default value provided")
+
+  nc <- nClass(
+      Cpublic = list(
+          go = nFunction(
+              fun = function(x = 'numericVector', start = 'numericScalar') {
+                  y <- parallel_reduce('pairmin', x, start)
+                  return(y)
+              },
+              returnType = 'numericScalar'
+          )
+      )
+  )
+  expect_error(Cnc <- nCompile(nc), "must be a literal")
+
+  go = nFunction(
+      fun = function(x = 'numericVector') {
+          y <- parallel_reduce('+', x)
+          return(y)
+      },
+      returnType = 'numericScalar'
+  )
+  expect_error(Cgo <- nCompile(go), "must be used in a method of an nClass")
+
+})
+
+## Could add check for user-defined reduction function with defined default init via operatorDef.
+
+test_that("user-defined reduction functions", {
+    ## User-defined nFunction
+    reduction_fun <- nFunction(
+        fun = function(x = 'numericScalar', y = 'numericScalar') {
+            ans <- x + y
+            return(ans)
+        },
+        returnType = 'numericScalar'
+    )
+    assign("reduction_fun", reduction_fun, envir = .GlobalEnv)
+    
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- parallel_reduce('reduction_fun', x, 0) 
+                    return(y)
+                },
+                returnType = 'numericScalar'
+            )
+        )
+    )
+    Cnc <- nCompile(nc, reduction_fun)[[1]]
+    obj <- nc$new()
+    Cobj <- Cnc$new()
+    expect_identical(obj$go(101:110), as.numeric(sum(101:110)))
+    expect_identical(Cobj$go(101:110), as.numeric(sum(101:110)))
+
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- parallel_reduce(reduction_fun, x, 0) 
+                    return(y)
+                },
+                returnType = 'numericScalar'
+            )
+        )
+    )
+    Cnc <- nCompile(nc, reduction_fun)[[1]]
+    obj <- nc$new()
+    Cobj <- Cnc$new()
+    expect_identical(obj$go(101:110), as.numeric(sum(101:110)))
+    expect_identical(Cobj$go(101:110), as.numeric(sum(101:110)))
+
+    ## See issue 133. This should work now.
+    nc <- nClass(
+        Cpublic = list(
+            reduction_fun = nFunction(
+                fun = function(x = 'numericScalar', y = 'numericScalar') {
+                    ans <- x + y
+                    return(ans)
+                },
+                returnType = 'numericScalar'
+            ),
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- parallel_reduce(reduction_fun, x, 0) 
+                    return(y)
+                },
+                returnType = 'numericScalar'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+    obj <- nc$new()
+    Cobj = Cnc$new()
+    expect_identical(obj$go(1:5), 15)
+    expect_identical(Cobj$go(1:5), 15)
+    
+    nc0 <- nClass(
+        Cpublic = list(
+            reduction_fun = nFunction(
+                fun = function(x = 'numericScalar', y = 'numericScalar') {
+                    ans <- x + y
+                    return(ans)
+                },
+                returnType = 'numericScalar'
+            )
+        ))
+    
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector', obj = 'nc0') {
+                    y <- parallel_reduce(obj$reduction_fun, x, 0) 
+                    return(y)
+                },
+                returnType = 'numericScalar'
+            )
+        )
+    )
+
+    obj0 <- nc0$new()
+    obj = nc$new()
+    expect_identical(obj$go(1:5, obj0), 15)
+    
+    Cnc <- nCompile(nc, nc0)
+    Cobj0 <- Cnc[[2]]$new()
+    Cobj = Cnc[[1]]$new()
+    expect_identical(Cobj$go(1:5, Cobj0), 15)
+
+    nc <- nClass(
+        Cpublic = list(
+            obj = 'nc0',
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- parallel_reduce(obj$reduction_fun, x, 0) 
+                    return(y)
+                },
+                returnType = 'numericScalar'
+            )
+        )
+    )
+    obj <- nc$new()
+    obj$obj <- nc0$new()
+    expect_identical(obj$go(1:5), 15)
+    
+    Cnc <- nCompile(nc, nc0)
+    Cobj = Cnc[[1]]$new()
+    expect_identical(Cobj$go(1:5), 15)
+    
+    nc <- nClass(
+        Cpublic = list(
+            obj = 'nc0',
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    obj <<- nc0$new()
+                    y <- parallel_reduce(obj$reduction_fun, x, 0) 
+                    return(y)
+                },
+                returnType = 'numericScalar'
+            )
+        )
+    )
+    obj <- nc$new()
+    expect_identical(obj$go(1:5), 15)
+    
+    Cnc <- nCompile(nc, nc0)
+    Cobj = Cnc[[1]]$new()
+    expect_identical(Cobj$go(1:5), 15)
+
+    # Using `self`. Currently we can't handle `self$x` as the vector just as we can't handle `obj$x`.
+    # We also can't handle the `init` not being a literal.
+    nc <- nClass(
+      Cpublic = list(
+        reduction_fun = nFunction(
+          fun = function(x = 'numericScalar', y = 'numericScalar') {
+            ans <- x + y
+            return(ans)
+          },
+          returnType = 'numericScalar'
+        ),
+        go = nFunction(
+          fun = function(x = 'numericVector') {
+            y <- parallel_reduce(self$reduction_fun, x, 0)
+            return(y)
+          },
+          returnType = 'numericScalar'
+        )
+      )
+    )
+    Cnc <- nCompile(nc)
+    obj <- nc$new()
+    Cobj = Cnc$new()
+    expect_identical(obj$go(1:5), 15)
+    expect_identical(Cobj$go(1:5), 15)
+
+})
+
+
+test_that("reduction cases that don't work", {
+  ## This doesn't work at present but we should make it work, presumably by lifting the `object` expression.  
+  nc <- nClass(
+    Cpublic = list(
+      go = nFunction(
+        fun = function(x = 'numericVector', y = 'numericVector') {
+          z <- parallel_reduce('+', x+y)
+          return(z)
+        },
+        returnType = 'numericScalar'
+      )
+    )
+  )
+  expect_error(Cnc <- nCompile(nc), 'found an expression')
+
+  ## Similar issue for this use case.
+  nc1 <- nClass(
+      Cpublic = list(
+          x = 'numericVector'
+      )
+  )
+  
+  nc <- nClass(
+      Cpublic = list(
+          go = nFunction(
+              fun = function(o = 'nc1') {
+                  y <- parallel_reduce('+', o$x, 0)
+                  return(y)
+              },
+              returnType = 'numericScalar'
+          )
+      )
+  )
+  expect_error(Cnc <- nCompile(nc,nc1), 'found an expression')
+  
+  ## Issue 136. This should work, but some type issue.
+  nc <- nClass(
+    Cpublic = list(
+      go = nFunction(
+        fun = function(x = 'integerVector') {
+          z <- parallel_reduce('pairmin', x)
+          return(z)
+        },
+        returnType = 'integerScalar'
+      )
+    )
+  )
+  expect_error(out <- capture.output(nCompile(nc)))  ## Lots of C++ compiler output.
+
+  nc1 <- nClass(
+      Cpublic = list(
+          y = 'nc2'
+    ))
+
+  nc2 <- nClass(
+      Cpublic = list(
+          plus = nFunction(
+              fun = function(x = 'numericScalar', y = 'numericScalar') {
+                  ans <- x + y
+                  return(ans)
+              },
+              returnType = 'numericScalar'
+          )
+      )
+  )
+  
+  nc <- nClass(
+      Cpublic = list(
+          go = nFunction(
+              fun = function(x = 'numericVector', o = 'nc1') {
+                  y <- parallel_reduce(o$y$plus, x)
+                  return(y)
+              },
+              returnType = 'numericScalar'
+          )
+      )
+  )
+  expect_error(Cnc <- nCompile(nc,nc1, nc2),  'too many levels of class hierarchy')
+
+})
+
+## Not clear why one would do this, but we do allow it.
+test_that("reduction without assignment", {
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    tmp <- 7
+                    parallel_reduce('+', x)
+                    return(0)
+                },
+                returnType = 'numericScalar'
+            )
+        )
+    )
+    obj <- nc$new()
+    expect_identical(obj$go(1:3), 0)
+    
+    Cnc <- nCompile(nc)
+    Cobj <- Cnc$new()
+    expect_identical(Cobj$go(1:3), 0)
+    
+    
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    3 + exp(parallel_reduce('+', x))
+                    return(0)
+                },
+                returnType = 'numericScalar'
+            )
+        )
+    )
+    obj <- nc$new()
+    expect_identical(obj$go(1:3), 0)
+    Cnc <- nCompile(nc)
+    Cobj <- Cnc$new()
+    expect_identical(Cobj$go(1:3), 0)
+})
+
+test_that("multiple reduction functions", {
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector', z = 'numericVector') {
+                    y <- parallel_reduce('+', x) + 3*parallel_reduce('+',z)
+                    return(y)
+                },
+                returnType = 'numericScalar'
+            )
+        )
+    )
+    obj <- nc$new()
+    Cnc <- nCompile(nc)
+    Cobj <- Cnc$new()
+    expect_identical(obj$go(1:3, 4:7), sum(1:3)+3*sum(4:7))
+    expect_identical(Cobj$go(1:3, 4:7), sum(1:3)+3*sum(4:7))
+
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(x = 'numericVector', z = 'numericVector') {
+                    y <- parallel_reduce('+', x) 
+                    return(y + 3*parallel_reduce('+',z,0))
+                },
+                returnType = 'numericScalar'
+            )
+        )
+    )
+    obj <- nc$new()
+    Cnc <- nCompile(nc)
+    Cobj <- Cnc$new()
+    expect_identical(obj$go(1:3, 4:7), sum(1:3)+3*sum(4:7))
+    expect_identical(Cobj$go(1:3, 4:7), sum(1:3)+3*sum(4:7))
+
+    # Avoid issue perhaps caused by persistence of TBB arena/pool/global_control object,
+    # resulting in `longjmp causes uninitialized stack frame ***: terminated`
+    rm(Cnc, Cobj); gc()
+})
+
+test_that("nested reduction functions", {
+    nc <- nClass(
+        Cpublic = list(
+            adder = nFunction(
+                fun = function(x = 'numericScalar', y = 'numericScalar') {
+                    z <- 1:5
+                    tmp <- parallel_reduce('+',z)
+                    ans <- x + y + tmp
+                    return(ans)
+                },
+                returnType = 'numericScalar'
+            ),
+            go = nFunction(
+                fun = function(x = 'numericVector') {
+                    y <- parallel_reduce(adder, x, 0)
+                    return(y)
+                },
+                returnType = 'numericScalar'
+            )
+        )
+    )
+    obj <- nc$new() 
+    Cnc <- nCompile(nc)
+    Cobj = Cnc$new()
+    expect_identical(obj$go(1:5), 90) # Note that this would not work if use `parallel_reduce('adder', x, 0)`.
+    expect_identical(Cobj$go(1:5), 90)
+    
+})
+
+
+
+
diff --git a/nCompiler/tests/testthat/tbb_tests/test-threads.R b/nCompiler/tests/testthat/tbb_tests/test-threads.R
new file mode 100644
index 00000000..aa4d3e9c
--- /dev/null
+++ b/nCompiler/tests/testthat/tbb_tests/test-threads.R
@@ -0,0 +1,145 @@
+## Note that it is hard to test that the number of threads is actually being
+## modified as on Linux TBB is not scaling well so testing timing won't be
+## useful. Therefore these tests mostly check that code compiles and runs.
+## But the design of the parallelized code is such that we could check
+## run time in the future.
+
+test_that("basic usage of setting number of threads", {
+
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(n = 'numericScalar', m = 'numericScalar') {
+                    y <- numeric(length=n)
+                    parallel_for(i, 1:n, {y[i] <- mean(rnorm(m))})
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+
+    Cnc1 <- Cnc$new()
+    expect_silent(out <- Cnc1$go(100, 100))
+    
+    val <- set_nOption('nThreads', 2)
+    expect_silent(out <- Cnc1$go(100, 100))
+    
+    val <- set_nOption('nThreads', 100000)
+    expect_silent(out <- Cnc1$go(100, 100))
+
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(n = 'numericScalar', m = 'numericScalar') {
+                    y <- numeric(length=n)
+                    parallel_for(i, 1:n, {y[i] <- mean(rnorm(m))}, nThreads = 2)
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+
+    Cnc1 <- Cnc$new()
+    expect_silent(out <- Cnc1$go(100, 100))
+
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(n = 'numericScalar', m = 'numericScalar', nThreads = 'numericScalar') {
+                    y <- numeric(length=n)
+                    parallel_for(i, 1:n, {y[i] <- mean(rnorm(m))}, nThreads = nThreads)
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+
+    Cnc1 <- Cnc$new()
+    expect_silent(out <- Cnc1$go(100, 100, 2))
+    expect_silent(out <- Cnc1$go(100, 100, 0))
+    val <- set_nOption('nThreads', 4)
+    expect_silent(out <- Cnc1$go(100, 100, 2))
+
+    ## Should execute in R but will ignore threads argument.
+    nc1 <- nc$new()
+    expect_silent(out <- nc1$go(100, 100, 2))
+
+    ## Multiple loops with different numbers of threads.
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(n = 'numericScalar', m = 'numericScalar') { 
+                    y <- numeric(length=n)
+                    parallel_for(i, 1:n,
+                    {y[i] <- mean(rnorm(m))}, nThreads=2)
+                    
+                    parallel_for(i, 1:n,
+                    {y[i] <- mean(rnorm(m))}, nThreads=8)
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    
+    Cnc <- nCompile(nc)
+    
+    Cnc1 <- Cnc$new()
+    expect_silent(out <- Cnc1$go(100, 100))
+
+    nc <- nClass(
+        Cpublic = list(
+            go = nFunction(
+                fun = function(n = 'numericScalar', m = 'numericScalar', nThreads = 'numericScalar') {
+                    y <- numeric(length=n)
+                    parallel_for(i, 1:n, {y[i] <- mean(rnorm(m))}, nThreads = nThreads+2)
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    Cnc <- nCompile(nc)
+
+    Cnc1 <- Cnc$new()
+    expect_silent(out <- Cnc1$go(100, 100, 2))
+
+    ## Nested loops
+    nc <- nClass(
+        Cpublic = list(
+            myfun = nFunction(
+                fun = function(m = 'numericScalar', k = 'numericScalar') {
+                    y <- numeric(length=m)
+                    parallel_for(i, 1:m, {y[i] = mean(rnorm(k))}, nThreads = 4)
+                    return(y)
+                }, returnType = 'numericVector'
+            ),                                 
+            go = nFunction(
+                fun = function(n = 'numericScalar', m = 'numericScalar', k = 'numericScalar') { 
+                    y <- numeric(length=n)
+                    parallel_for(i, 1:n, {
+                        y[i] <- mean(myfun(m, k))
+                    }, nThreads=2)
+                    return(y)
+                },
+                returnType = 'numericVector'
+            )
+        )
+    )
+    
+    Cnc <- nCompile(nc)
+    
+    Cnc1 <- Cnc$new()
+    ## Issue 152.
+    ## expect_silent(out <- Cnc1$go(100, 100, 10))
+    ## Error: C stack usage  14983114845556 is too close to the limit
+    ## Error: no more error handlers available (recursive errors?); invoking 'abort' restart
+    ## *** longjmp causes uninitialized stack frame ***: terminated
+
+})