diff --git a/Addons/testLoadUnstructured.hpp b/Addons/testLoadUnstructured.hpp
index ab184ffa2..2f7a0430c 100644
--- a/Addons/testLoadUnstructured.hpp
+++ b/Addons/testLoadUnstructured.hpp
@@ -88,9 +88,10 @@ inline double coefficientDifference(TasmanianSparseGrid const &gA, TasmanianSpar
     size_t num_coeffs = Utils::size_mult(gA.getNumOutputs(), gA.getNumPoints());
     double const *c1 = gA.getHierarchicalCoefficients();
     double const *c2 = gB.getHierarchicalCoefficients();
-    return std::inner_product(c1, c1 + num_coeffs, c2, 0.0,
-                              [](double a, double b)->double{ return std::max(a, b); },
-                              [](double a, double b)->double{ return std::abs(a - b); });
+    double err = 0.0;
+    for(size_t i=0; i<num_coeffs; i++)
+        err = std::max(std::abs(c1[i] - c2[i]), err);
+    return err;
 }
 
 /*!
@@ -102,16 +103,18 @@ inline double evalDifference(std::vector<double> const&x, TasmanianSparseGrid co
     gA.evaluateBatch(x, yA);
     gB.evaluateBatch(x, yB);
 
-    return std::inner_product(yA.begin(), yA.end(), yB.begin(), 0.0,
-                              [](double a, double b)->double{ return std::max(a, b); },
-                              [](double a, double b)->double{ return std::abs(a - b); });
+    double err = 0.0;
+    for(size_t i=0; i<yA.size(); i++)
+        err = std::max(std::abs(yA[i] - yB[i]), err);
+    return err;
 }
 
 inline double vecDifference(std::vector<double> const &x, std::vector<double> const &y){
     if (x.size() != y.size()) return 1.E+20;
-    return std::inner_product(x.begin(), x.end(), y.begin(), 0.0,
-                              [](double a, double b)->double{ return std::max(a, b); },
-                              [](double a, double b)->double{ return std::abs(a - b); });
+    double err = 0.0;
+    for(size_t i=0; i<x.size(); i++)
+        err = std::max(std::abs(x[i] - y[i]), err);
+    return err;
 }
 
 /*!
diff --git a/SparseGrids/gridtestExternalTests.cpp b/SparseGrids/gridtestExternalTests.cpp
index 3f2de0a74..639e1e68c 100644
--- a/SparseGrids/gridtestExternalTests.cpp
+++ b/SparseGrids/gridtestExternalTests.cpp
@@ -1041,10 +1041,11 @@ bool ExternalTester::testDynamicRefinement(const BaseFunction *f, TasmanianSpars
     std::vector<double> xpnts = genRandom(10, grid.getNumDimensions());
     grid.evaluateBatch(xpnts, res1);
     grid2.evaluateBatch(xpnts, res2);
-    double err = std::inner_product(res1.begin(), res1.end(), res2.begin(), 0.0,
-                                   [](double a, double b)->double{ return std::max(a, b); },
-                                   [](double a, double b)->double{ return std::abs(a - b); });
-    if (err > Maths::num_tol){ cout << "ERROR: failed evaluate after loading batch points." << endl; return false; }
+    double err = 0.0;
+    for(size_t i=0; i<res1.size(); i++)
+        err = std::max(std::abs(res1[i] - res2[i]), err);
+
+    if (err > Maths::num_tol){ cout << "ERROR: failed evaluate after loading batch points. " << endl; return false; }
     return true;
 }
 
diff --git a/SparseGrids/tsgGridLocalPolynomial.cpp b/SparseGrids/tsgGridLocalPolynomial.cpp
index 3d62a71ba..62f75afe6 100644
--- a/SparseGrids/tsgGridLocalPolynomial.cpp
+++ b/SparseGrids/tsgGridLocalPolynomial.cpp
@@ -1266,83 +1266,117 @@ Data2D<int> GridLocalPolynomial::buildUpdateMap(double tolerance, TypeRefinement
         int max_1D_parents = rule->getMaxNumParents();
 
         HierarchyManipulations::SplitDirections split(points);
-        std::vector<int> global_to_pnts(num_points);
-        #pragma omp parallel for firstprivate(global_to_pnts)
-        for(int j=0; j<split.getNumJobs(); j++){ // split.getNumJobs() gives the number of 1D interpolants to construct
-            int d = split.getJobDirection(j);
-            int nump = split.getJobNumPoints(j);
-            const int *pnts = split.getJobPoints(j);
 
-            std::vector<int> levels(nump);
+        #pragma omp parallel
+        {
+            int max_nump = split.getMaxNumPoints();
+            std::vector<int> global_to_pnts(num_points);
+            std::vector<int> levels(max_nump);
+
+            std::vector<int> monkey_count;
+            std::vector<int> monkey_tail;
+            std::vector<bool> used;
+            if (max_1D_parents > 1) {
+                monkey_count.resize(top_level + 1);
+                monkey_tail.resize(top_level + 1);
+                used.resize(max_nump, false);
+            }
 
-            int max_level = 0;
+            #pragma omp for
+            for(int j=0; j<split.getNumJobs(); j++){ // split.getNumJobs() gives the number of 1D interpolants to construct
+                int d = split.getJobDirection(j);
+                int nump = split.getJobNumPoints(j);
+                const int *pnts = split.getJobPoints(j);
 
-            Data2D<double> vals(active_outputs, nump);
+                int max_level = 0;
 
-            for(int i=0; i<nump; i++){
-                const double* v = values.getValues(pnts[i]);
-                const int *p = points.getIndex(pnts[i]);
-                if (output == -1){
-                    std::copy(v, v + num_outputs, vals.getStrip(i));
-                }else{
-                    vals.getStrip(i)[0] = v[output];
+                Data2D<double> vals(active_outputs, nump);
+
+                if (output == -1) {
+                    for(int i=0; i<nump; i++){
+                        std::copy_n(values.getValues(pnts[i]), num_outputs, vals.getStrip(i));
+                        global_to_pnts[pnts[i]] = i;
+                        levels[i] = rule->getLevel(points.getIndex(pnts[i])[d]);
+                        if (max_level < levels[i]) max_level = levels[i];
+                    }
+                } else {
+                    for(int i=0; i<nump; i++){
+                        vals.getStrip(i)[0] = values.getValues(pnts[i])[output];
+                        global_to_pnts[pnts[i]] = i;
+                        levels[i] = rule->getLevel(points.getIndex(pnts[i])[d]);
+                        if (max_level < levels[i]) max_level = levels[i];
+                    }
                 }
-                global_to_pnts[pnts[i]] = i;
-                levels[i] = rule->getLevel(p[d]);
-                if (max_level < levels[i]) max_level = levels[i];
-            }
 
-            std::vector<int> monkey_count(max_level + 1);
-            std::vector<int> monkey_tail(max_level + 1);
+                if (max_1D_parents == 1) {
+                    for(int l=1; l<=max_level; l++){
+                        for(int i=0; i<nump; i++){
+                            if (levels[i] == l){
+                                double x = rule->getNode(points.getIndex(pnts[i])[d]);
+                                double *valsi = vals.getStrip(i);
 
-            for(int l=1; l<=max_level; l++){
-                for(int i=0; i<nump; i++){
-                    if (levels[i] == l){
-                        double x = rule->getNode(points.getIndex(pnts[i])[d]);
-                        double *valsi = vals.getStrip(i);
-
-                        int current = 0;
-                        monkey_count[0] = d * max_1D_parents;
-                        monkey_tail[0] = pnts[i]; // uses the global indexes
-                        std::vector<bool> used(nump, false);
-
-                        while(monkey_count[0] < (d+1) * max_1D_parents){
-                            if (monkey_count[current] < (d+1) * max_1D_parents){
-                                int branch = dagUp.getStrip(monkey_tail[current])[monkey_count[current]];
-                                if ((branch == -1) || (used[global_to_pnts[branch]])){
-                                    monkey_count[current]++;
-                                }else{
+                                int branch = dagUp.getStrip(pnts[i])[d];
+                                while(branch != -1) {
                                     const int *branch_point = points.getIndex(branch);
                                     double basis_value = rule->evalRaw(branch_point[d], x);
                                     const double *branch_vals = vals.getStrip(global_to_pnts[branch]);
                                     for(int k=0; k<active_outputs; k++) valsi[k] -= basis_value * branch_vals[k];
-
-                                    used[global_to_pnts[branch]] = true;
-                                    monkey_count[++current] = d * max_1D_parents;
-                                    monkey_tail[current] = branch;
+                                    branch = dagUp.getStrip(branch)[d];
+                                }
+                            }
+                        }
+                    }
+                } else {
+                    for(int l=1; l<=max_level; l++){
+                        for(int i=0; i<nump; i++){
+                            if (levels[i] == l){
+                                double x = rule->getNode(points.getIndex(pnts[i])[d]);
+                                double *valsi = vals.getStrip(i);
+
+                                int current = 0;
+                                monkey_count[0] = d * max_1D_parents;
+                                monkey_tail[0] = pnts[i]; // uses the global indexes
+                                std::fill_n(used.begin(), nump, false);
+
+                                while(monkey_count[0] < (d+1) * max_1D_parents){
+                                    if (monkey_count[current] < (d+1) * max_1D_parents){
+                                        int branch = dagUp.getStrip(monkey_tail[current])[monkey_count[current]];
+                                        if ((branch == -1) || (used[global_to_pnts[branch]])){
+                                            monkey_count[current]++;
+                                        }else{
+                                            const int *branch_point = points.getIndex(branch);
+                                            double basis_value = rule->evalRaw(branch_point[d], x);
+                                            const double *branch_vals = vals.getStrip(global_to_pnts[branch]);
+                                            for(int k=0; k<active_outputs; k++) valsi[k] -= basis_value * branch_vals[k];
+
+                                            used[global_to_pnts[branch]] = true;
+                                            monkey_count[++current] = d * max_1D_parents;
+                                            monkey_tail[current] = branch;
+                                        }
+                                    }else{
+                                        monkey_count[--current]++;
+                                    }
                                 }
-                            }else{
-                                monkey_count[--current]++;
                             }
                         }
                     }
                 }
-            }
 
-            // at this point, vals contains the one directional surpluses
-            for(int i=0; i<nump; i++){
-                const double *s = surpluses.getStrip(pnts[i]);
-                const double *c = scale.getStrip(pnts[i]);
-                const double *v = vals.getStrip(i);
-                bool small = true;
-                if (output == -1){
-                    for(int k=0; k<num_outputs; k++){
-                        small = small && (((c[k] * std::abs(s[k]) / norm[k]) <= tolerance) || ((c[k] * std::abs(v[k]) / norm[k]) <= tolerance));
+                // at this point, vals contains the one directional surpluses
+                for(int i=0; i<nump; i++){
+                    const double *s = surpluses.getStrip(pnts[i]);
+                    const double *c = scale.getStrip(pnts[i]);
+                    const double *v = vals.getStrip(i);
+                    bool small = true;
+                    if (output == -1){
+                        for(int k=0; k<num_outputs; k++){
+                            small = small && (((c[k] * std::abs(s[k]) / norm[k]) <= tolerance) || ((c[k] * std::abs(v[k]) / norm[k]) <= tolerance));
+                        }
+                    }else{
+                        small = ((c[0] * std::abs(s[output]) / norm[output]) <= tolerance) || ((c[0] * std::abs(v[0]) / norm[output]) <= tolerance);
                     }
-                }else{
-                    small = ((c[0] * std::abs(s[output]) / norm[output]) <= tolerance) || ((c[0] * std::abs(v[0]) / norm[output]) <= tolerance);
+                    pmap.getStrip(pnts[i])[d] = (small) ? 0 : 1;;
                 }
-                pmap.getStrip(pnts[i])[d] = (small) ? 0 : 1;;
             }
         }
     }
@@ -1357,28 +1391,40 @@ MultiIndexSet GridLocalPolynomial::getRefinementCanidates(double tolerance, Type
 
     int num_points = points.getNumIndexes();
 
-    if (level_limits.empty()){
-        for(int i=0; i<num_points; i++){
-            const int *map = pmap.getStrip(i);
-            for(int j=0; j<num_dimensions; j++){
-                if (map[j] == 1){ // if this dimension needs to be refined
-                    if (!(useParents && addParent(points.getIndex(i), j, points, refined))){
-                        addChild(points.getIndex(i), j, points, refined);
+    #pragma omp parallel
+    {
+        Data2D<int> lrefined(num_dimensions, 0);
+
+        if (level_limits.empty()){
+            #pragma omp for
+            for(int i=0; i<num_points; i++){
+                const int *map = pmap.getStrip(i);
+                for(int j=0; j<num_dimensions; j++){
+                    if (map[j] == 1){ // if this dimension needs to be refined
+                        if (!(useParents && addParent(points.getIndex(i), j, points, lrefined))){
+                            addChild(points.getIndex(i), j, points, lrefined);
+                        }
                     }
                 }
             }
-        }
-    }else{
-        for(int i=0; i<num_points; i++){
-            const int *map = pmap.getStrip(i);
-            for(int j=0; j<num_dimensions; j++){
-                if (map[j] == 1){ // if this dimension needs to be refined
-                    if (!(useParents && addParent(points.getIndex(i), j, points, refined))){
-                        addChildLimited(points.getIndex(i), j, points, level_limits, refined);
+        }else{
+            #pragma omp for
+            for(int i=0; i<num_points; i++){
+                const int *map = pmap.getStrip(i);
+                for(int j=0; j<num_dimensions; j++){
+                    if (map[j] == 1){ // if this dimension needs to be refined
+                        if (!(useParents && addParent(points.getIndex(i), j, points, lrefined))){
+                            addChildLimited(points.getIndex(i), j, points, level_limits, lrefined);
+                        }
                     }
                 }
             }
         }
+
+        #pragma omp critical
+        {
+            refined.append(lrefined);
+        }
     }
 
     MultiIndexSet result(refined);
diff --git a/SparseGrids/tsgHierarchyManipulator.hpp b/SparseGrids/tsgHierarchyManipulator.hpp
index 6815bc454..3c6f4eb72 100644
--- a/SparseGrids/tsgHierarchyManipulator.hpp
+++ b/SparseGrids/tsgHierarchyManipulator.hpp
@@ -309,6 +309,9 @@ class SplitDirections{
     int getJobNumPoints(int job) const{ return (int) job_pnts[job].size(); }
     //! \brief Return the indexes of the points associated with the job.
     const int* getJobPoints(int job) const{ return job_pnts[job].data(); }
+    //! \brief Return the max number of points for any job.
+    int getMaxNumPoints() const { return (job_pnts.size() > 0) ? (int) std::max_element(job_pnts.begin(), job_pnts.end(),
+        [&](std::vector<int> const &a, std::vector<int> const& b)->bool{ return (a.size() < b.size()); })->size() : 0; }
 
 private:
     std::vector<int> job_directions;
diff --git a/SparseGrids/tsgIndexSets.hpp b/SparseGrids/tsgIndexSets.hpp
index 8c6280ba4..3c238cb76 100644
--- a/SparseGrids/tsgIndexSets.hpp
+++ b/SparseGrids/tsgIndexSets.hpp
@@ -192,6 +192,12 @@ class Data2D{
         num_strips++;
     }
 
+    //! \brief Uses std::vector::insert to append all the data from the other to this
+    void append(Data2D<T> const &other) {
+        vec.insert(vec.end(), other.vec.begin(), other.vec.end());
+        num_strips += other.num_strips;
+    }
+
 private:
     size_t stride, num_strips;
     std::vector<T> vec;
diff --git a/SparseGrids/tsgSequenceOptimizer.cpp b/SparseGrids/tsgSequenceOptimizer.cpp
index f44889941..61b6727b5 100644
--- a/SparseGrids/tsgSequenceOptimizer.cpp
+++ b/SparseGrids/tsgSequenceOptimizer.cpp
@@ -346,8 +346,11 @@ template<> double getValue<rule_mindeltaodd>(CurrentNodes<rule_mindeltaodd> cons
     auto lag       = evalLagrange(current.nodes, current.coeff, x);
     auto lag_less1 = evalLagrange(current.nodes_less1, current.coeff_less1, x);
 
-    return std::abs(lag.back()) + std::inner_product(lag_less1.begin(), lag_less1.end(), lag.begin(), 0.0,
-                                     std::plus<double>(), [](double a, double b)->double{ return std::abs(a - b); });
+    double result = std::abs(lag.back());
+    for(size_t i=0; i<lag_less1.size(); i++) {
+        result += std::abs(lag_less1[i] - lag[i]);
+    }
+    return result;
 }
 
 /*!