diff --git a/bazel/ceres.bzl b/bazel/ceres.bzl
index 6ba0137..8963370 100644
--- a/bazel/ceres.bzl
+++ b/bazel/ceres.bzl
@@ -91,6 +91,7 @@
     "minimizer.cc",
     "normal_prior.cc",
     "parallel_for_cxx.cc",
+    "parallel_for_openmp.cc",
     "parallel_for_tbb.cc",
     "parallel_utils.cc",
     "parameter_block_ordering.cc",
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 8924173..1dad396 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -91,6 +91,7 @@
     low_rank_inverse_hessian.cc
     minimizer.cc
     normal_prior.cc
+    parallel_for_openmp.cc
     parallel_for_cxx.cc
     parallel_for_tbb.cc
     parallel_utils.cc
diff --git a/internal/ceres/coordinate_descent_minimizer.cc b/internal/ceres/coordinate_descent_minimizer.cc
index 6f20f20..558faed 100644
--- a/internal/ceres/coordinate_descent_minimizer.cc
+++ b/internal/ceres/coordinate_descent_minimizer.cc
@@ -30,10 +30,6 @@
 
 #include "ceres/coordinate_descent_minimizer.h"
 
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-#include "ceres/parallel_for.h"
-#endif
-
 #include <algorithm>
 #include <iterator>
 #include <memory>
@@ -43,14 +39,13 @@
 #include "ceres/evaluator.h"
 #include "ceres/linear_solver.h"
 #include "ceres/minimizer.h"
+#include "ceres/parallel_for.h"
 #include "ceres/parameter_block.h"
 #include "ceres/parameter_block_ordering.h"
 #include "ceres/problem_impl.h"
 #include "ceres/program.h"
 #include "ceres/residual_block.h"
-#include "ceres/scoped_thread_token.h"
 #include "ceres/solver.h"
-#include "ceres/thread_token_provider.h"
 #include "ceres/trust_region_minimizer.h"
 #include "ceres/trust_region_strategy.h"
 
@@ -164,62 +159,45 @@
     evaluator_options_.num_threads =
         max(1, options.num_threads / num_inner_iteration_threads);
 
-    ThreadTokenProvider thread_token_provider(num_inner_iteration_threads);
-
-#ifdef CERES_USE_OPENMP
     // The parameter blocks in each independent set can be optimized
     // in parallel, since they do not co-occur in any residual block.
-#pragma omp parallel for num_threads(num_inner_iteration_threads)
-#endif
+    ParallelFor(
+        context_,
+        independent_set_offsets_[i],
+        independent_set_offsets_[i + 1],
+        num_inner_iteration_threads,
+        [&](int thread_id, int j) {
+          ParameterBlock* parameter_block = parameter_blocks_[j];
+          const int old_index = parameter_block->index();
+          const int old_delta_offset = parameter_block->delta_offset();
+          parameter_block->SetVarying();
+          parameter_block->set_index(0);
+          parameter_block->set_delta_offset(0);
 
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-    for (int j = independent_set_offsets_[i];
-         j < independent_set_offsets_[i + 1];
-         ++j) {
-#else
-    ParallelFor(context_,
-                independent_set_offsets_[i],
-                independent_set_offsets_[i + 1],
-                num_inner_iteration_threads,
-                [&](int j) {
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
+          Program inner_program;
+          inner_program.mutable_parameter_blocks()->push_back(parameter_block);
+          *inner_program.mutable_residual_blocks() = residual_blocks_[j];
 
-      const ScopedThreadToken scoped_thread_token(&thread_token_provider);
-      const int thread_id = scoped_thread_token.token();
+          // TODO(sameeragarwal): Better error handling. Right now we
+          // assume that this is not going to lead to problems of any
+          // sort. Basically we should be checking for numerical failure
+          // of some sort.
+          //
+          // On the other hand, if the optimization is a failure, that in
+          // some ways is fine, since it won't change the parameters and
+          // we are fine.
+          Solver::Summary inner_summary;
+          Solve(&inner_program,
+                linear_solvers[thread_id],
+                parameters + parameter_block->state_offset(),
+                &inner_summary);
 
-      ParameterBlock* parameter_block = parameter_blocks_[j];
-      const int old_index = parameter_block->index();
-      const int old_delta_offset = parameter_block->delta_offset();
-      parameter_block->SetVarying();
-      parameter_block->set_index(0);
-      parameter_block->set_delta_offset(0);
-
-      Program inner_program;
-      inner_program.mutable_parameter_blocks()->push_back(parameter_block);
-      *inner_program.mutable_residual_blocks() = residual_blocks_[j];
-
-      // TODO(sameeragarwal): Better error handling. Right now we
-      // assume that this is not going to lead to problems of any
-      // sort. Basically we should be checking for numerical failure
-      // of some sort.
-      //
-      // On the other hand, if the optimization is a failure, that in
-      // some ways is fine, since it won't change the parameters and
-      // we are fine.
-      Solver::Summary inner_summary;
-      Solve(&inner_program,
-            linear_solvers[thread_id],
-            parameters + parameter_block->state_offset(),
-            &inner_summary);
-
-      parameter_block->set_index(old_index);
-      parameter_block->set_delta_offset(old_delta_offset);
-      parameter_block->SetState(parameters + parameter_block->state_offset());
-      parameter_block->SetConstant();
-    }
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-  );
-#endif
+          parameter_block->set_index(old_index);
+          parameter_block->set_delta_offset(old_delta_offset);
+          parameter_block->SetState(parameters +
+                                    parameter_block->state_offset());
+          parameter_block->SetConstant();
+        });
   }
 
   for (int i =  0; i < parameter_blocks_.size(); ++i) {
diff --git a/internal/ceres/covariance_impl.cc b/internal/ceres/covariance_impl.cc
index 2552416..7c63c2c 100644
--- a/internal/ceres/covariance_impl.cc
+++ b/internal/ceres/covariance_impl.cc
@@ -30,10 +30,6 @@
 
 #include "ceres/covariance_impl.h"
 
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-#include "ceres/parallel_for.h"
-#endif
-
 #include <algorithm>
 #include <cstdlib>
 #include <memory>
@@ -53,13 +49,12 @@
 #include "ceres/crs_matrix.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/map_util.h"
+#include "ceres/parallel_for.h"
 #include "ceres/parallel_utils.h"
 #include "ceres/parameter_block.h"
 #include "ceres/problem_impl.h"
 #include "ceres/residual_block.h"
-#include "ceres/scoped_thread_token.h"
 #include "ceres/suitesparse.h"
-#include "ceres/thread_token_provider.h"
 #include "ceres/wall_time.h"
 #include "glog/logging.h"
 
@@ -344,60 +339,41 @@
 
   bool success = true;
 
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-  ThreadTokenProvider thread_token_provider(num_threads);
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-
   // Technically the following code is a double nested loop where
   // i = 1:n, j = i:n.
   int iteration_count = (num_parameters * (num_parameters + 1)) / 2;
-#if defined(CERES_USE_OPENMP)
-#    pragma omp parallel for num_threads(num_threads) schedule(dynamic)
-#endif // CERES_USE_OPENMP
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-  for (int k = 0; k < iteration_count; ++k) {
-#else
   problem_->context()->EnsureMinimumThreads(num_threads);
-  ParallelFor(problem_->context(),
-              0,
-              iteration_count,
-              num_threads,
-              [&](int thread_id, int k) {
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-      int i, j;
-      LinearIndexToUpperTriangularIndex(k, num_parameters, &i, &j);
+  ParallelFor(
+      problem_->context(),
+      0,
+      iteration_count,
+      num_threads,
+      [&](int thread_id, int k) {
+        int i, j;
+        LinearIndexToUpperTriangularIndex(k, num_parameters, &i, &j);
 
-      int covariance_row_idx = cum_parameter_size[i];
-      int covariance_col_idx = cum_parameter_size[j];
-      int size_i = parameter_sizes[i];
-      int size_j = parameter_sizes[j];
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-      const ScopedThreadToken scoped_thread_token(&thread_token_provider);
-      const int thread_id = scoped_thread_token.token();
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-      double* covariance_block =
-          workspace.get() +
-          thread_id * max_covariance_block_size * max_covariance_block_size;
-      if (!GetCovarianceBlockInTangentOrAmbientSpace(
-              parameters[i], parameters[j], lift_covariance_to_ambient_space,
-              covariance_block)) {
-        success = false;
-      }
+        int covariance_row_idx = cum_parameter_size[i];
+        int covariance_col_idx = cum_parameter_size[j];
+        int size_i = parameter_sizes[i];
+        int size_j = parameter_sizes[j];
+        double* covariance_block =
+            workspace.get() + thread_id * max_covariance_block_size *
+            max_covariance_block_size;
+        if (!GetCovarianceBlockInTangentOrAmbientSpace(
+                parameters[i], parameters[j],
+                lift_covariance_to_ambient_space, covariance_block)) {
+          success = false;
+        }
 
-      covariance.block(covariance_row_idx, covariance_col_idx,
-                       size_i, size_j) =
-          MatrixRef(covariance_block, size_i, size_j);
+        covariance.block(covariance_row_idx, covariance_col_idx, size_i,
+                         size_j) = MatrixRef(covariance_block, size_i, size_j);
 
-      if (i != j) {
-        covariance.block(covariance_col_idx, covariance_row_idx,
-                         size_j, size_i) =
-            MatrixRef(covariance_block, size_i, size_j).transpose();
-
-      }
-    }
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-   );
-#endif // defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
+        if (i != j) {
+          covariance.block(covariance_col_idx, covariance_row_idx,
+                           size_j, size_i) =
+              MatrixRef(covariance_block, size_i, size_j).transpose();
+        }
+      });
   return success;
 }
 
@@ -709,50 +685,27 @@
   const int num_threads = options_.num_threads;
   std::unique_ptr<double[]> workspace(new double[num_threads * num_cols]);
 
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-  ThreadTokenProvider thread_token_provider(num_threads);
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-
-#ifdef CERES_USE_OPENMP
-#pragma omp parallel for num_threads(num_threads) schedule(dynamic)
-#endif // CERES_USE_OPENMP
-
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-  for (int r = 0; r < num_cols; ++r) {
-#else
   problem_->context()->EnsureMinimumThreads(num_threads);
-  ParallelFor(problem_->context(),
-              0,
-              num_cols,
-              num_threads,
-              [&](int thread_id, int r) {
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-
-    const int row_begin = rows[r];
-    const int row_end = rows[r + 1];
-    if (row_end != row_begin) {
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-      const ScopedThreadToken scoped_thread_token(&thread_token_provider);
-      const int thread_id = scoped_thread_token.token();
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-
-      double* solution = workspace.get() + thread_id * num_cols;
-      SolveRTRWithSparseRHS<SuiteSparse_long>(
-          num_cols,
-          static_cast<SuiteSparse_long*>(R->i),
-          static_cast<SuiteSparse_long*>(R->p),
-          static_cast<double*>(R->x),
-          inverse_permutation[r],
-          solution);
-      for (int idx = row_begin; idx < row_end; ++idx) {
-        const int c = cols[idx];
-        values[idx] = solution[inverse_permutation[c]];
-      }
-    }
-  }
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-  );
-#endif // defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
+  ParallelFor(
+      problem_->context(),
+      0,
+      num_cols,
+      num_threads,
+      [&](int thread_id, int r) {
+        const int row_begin = rows[r];
+        const int row_end = rows[r + 1];
+        if (row_end != row_begin) {
+          double* solution = workspace.get() + thread_id * num_cols;
+          SolveRTRWithSparseRHS<SuiteSparse_long>(
+              num_cols, static_cast<SuiteSparse_long*>(R->i),
+              static_cast<SuiteSparse_long*>(R->p), static_cast<double*>(R->x),
+              inverse_permutation[r], solution);
+          for (int idx = row_begin; idx < row_end; ++idx) {
+            const int c = cols[idx];
+            values[idx] = solution[inverse_permutation[c]];
+          }
+        }
+      });
 
   free(permutation);
   cholmod_l_free_sparse(&R, &cc);
@@ -918,54 +871,33 @@
   const int num_threads = options_.num_threads;
   std::unique_ptr<double[]> workspace(new double[num_threads * num_cols]);
 
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-  ThreadTokenProvider thread_token_provider(num_threads);
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-
-#ifdef CERES_USE_OPENMP
-#pragma omp parallel for num_threads(num_threads) schedule(dynamic)
-#endif // CERES_USE_OPENMP
-
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-  for (int r = 0; r < num_cols; ++r) {
-#else
   problem_->context()->EnsureMinimumThreads(num_threads);
-  ParallelFor(problem_->context(),
-              0,
+  ParallelFor(
+      problem_->context(),
+      0,
+      num_cols,
+      num_threads,
+      [&](int thread_id, int r) {
+        const int row_begin = rows[r];
+        const int row_end = rows[r + 1];
+        if (row_end != row_begin) {
+          double* solution = workspace.get() + thread_id * num_cols;
+          SolveRTRWithSparseRHS<int>(
               num_cols,
-              num_threads,
-              [&](int thread_id, int r) {
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
+              qr_solver.matrixR().innerIndexPtr(),
+              qr_solver.matrixR().outerIndexPtr(),
+              &qr_solver.matrixR().data().value(0),
+              inverse_permutation.indices().coeff(r),
+              solution);
 
-    const int row_begin = rows[r];
-    const int row_end = rows[r + 1];
-    if (row_end != row_begin) {
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-      const ScopedThreadToken scoped_thread_token(&thread_token_provider);
-      const int thread_id = scoped_thread_token.token();
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-
-      double* solution = workspace.get() + thread_id * num_cols;
-      SolveRTRWithSparseRHS<int>(
-          num_cols,
-          qr_solver.matrixR().innerIndexPtr(),
-          qr_solver.matrixR().outerIndexPtr(),
-          &qr_solver.matrixR().data().value(0),
-          inverse_permutation.indices().coeff(r),
-          solution);
-
-      // Assign the values of the computed covariance using the
-      // inverse permutation used in the QR factorization.
-      for (int idx = row_begin; idx < row_end; ++idx) {
-       const int c = cols[idx];
-       values[idx] = solution[inverse_permutation.indices().coeff(c)];
-      }
-    }
-  }
-
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-  );
-#endif // defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
+          // Assign the values of the computed covariance using the
+          // inverse permutation used in the QR factorization.
+          for (int idx = row_begin; idx < row_end; ++idx) {
+            const int c = cols[idx];
+            values[idx] = solution[inverse_permutation.indices().coeff(c)];
+          }
+        }
+      });
 
   event_logger.AddEvent("Inverse");
 
diff --git a/internal/ceres/parallel_for_openmp.cc b/internal/ceres/parallel_for_openmp.cc
new file mode 100644
index 0000000..b64f3d7
--- /dev/null
+++ b/internal/ceres/parallel_for_openmp.cc
@@ -0,0 +1,83 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2018 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: vitus@google.com (Michael Vitus)
+
+// This include must come before any #ifndef check on Ceres compile options.
+#include "ceres/internal/port.h"
+
+#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
+
+#include "ceres/parallel_for.h"
+
+#include "ceres/scoped_thread_token.h"
+#include "ceres/thread_token_provider.h"
+#include "glog/logging.h"
+
+namespace ceres {
+namespace internal {
+
+void ParallelFor(ContextImpl* context,
+                 int start,
+                 int end,
+                 int num_threads,
+                 const std::function<void(int)>& function) {
+  CHECK_GT(num_threads, 0);
+  CHECK(context != NULL);
+  if (end <= start) {
+    return;
+  }
+
+#ifdef CERES_USE_OPENMP
+#pragma omp parallel for num_threads(num_threads) \
+    schedule(dynamic) if (num_threads > 1)
+#endif  // CERES_USE_OPENMP
+  for (int i = start; i < end; ++i) {
+    function(i);
+  }
+}
+
+void ParallelFor(ContextImpl* context,
+                 int start,
+                 int end,
+                 int num_threads,
+                 const std::function<void(int thread_id, int i)>& function) {
+  CHECK(context != NULL);
+
+  ThreadTokenProvider thread_token_provider(num_threads);
+  ParallelFor(context, start, end, num_threads, [&](int i) {
+    const ScopedThreadToken scoped_thread_token(&thread_token_provider);
+    const int thread_id = scoped_thread_token.token();
+    function(thread_id, i);
+  });
+}
+
+}  // namespace internal
+}  // namespace ceres
+
+#endif  // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
diff --git a/internal/ceres/parallel_for_test.cc b/internal/ceres/parallel_for_test.cc
index 2de4865..04e5783 100644
--- a/internal/ceres/parallel_for_test.cc
+++ b/internal/ceres/parallel_for_test.cc
@@ -31,8 +31,6 @@
 // This include must come before any #ifndef check on Ceres compile options.
 #include "ceres/internal/port.h"
 
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-
 #include "ceres/parallel_for.h"
 
 #include <cmath>
@@ -129,6 +127,8 @@
   }
 }
 
+// This test is only valid when multithreading support is enabled.
+#ifndef CERES_NO_THREADS
 TEST(ParallelForWithThreadId, UniqueThreadIds) {
   // Ensure the hardware supports more than 1 thread to ensure the test will
   // pass.
@@ -157,8 +157,7 @@
 
   EXPECT_THAT(x, UnorderedElementsAreArray({0,1}));
 }
+#endif  // CERES_NO_THREADS
 
 }  // namespace internal
 }  // namespace ceres
-
-#endif // defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
diff --git a/internal/ceres/program_evaluator.h b/internal/ceres/program_evaluator.h
index 17175a6..13f2e83 100644
--- a/internal/ceres/program_evaluator.h
+++ b/internal/ceres/program_evaluator.h
@@ -82,25 +82,20 @@
 // This include must come before any #ifndef check on Ceres compile options.
 #include "ceres/internal/port.h"
 
+#include <atomic>
 #include <map>
 #include <memory>
 #include <string>
 #include <vector>
+
 #include "ceres/evaluation_callback.h"
 #include "ceres/execution_summary.h"
 #include "ceres/internal/eigen.h"
+#include "ceres/parallel_for.h"
 #include "ceres/parameter_block.h"
 #include "ceres/program.h"
 #include "ceres/residual_block.h"
-#include "ceres/scoped_thread_token.h"
 #include "ceres/small_blas.h"
-#include "ceres/thread_token_provider.h"
-
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-#include <atomic>
-
-#include "ceres/parallel_for.h"
-#endif
 
 namespace ceres {
 namespace internal {
@@ -183,128 +178,84 @@
     }
 
     const int num_residual_blocks = program_->NumResidualBlocks();
-
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-    ThreadTokenProvider thread_token_provider(options_.num_threads);
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-
-#ifdef CERES_USE_OPENMP
-    // This bool is used to disable the loop if an error is encountered
-    // without breaking out of it. The remaining loop iterations are still run,
-    // but with an empty body, and so will finish quickly.
-    bool abort = false;
-#pragma omp parallel for num_threads(options_.num_threads)
-    for (int i = 0; i < num_residual_blocks; ++i) {
-// Disable the loop instead of breaking, as required by OpenMP.
-#pragma omp flush(abort)
-#endif // CERES_USE_OPENMP
-
-#ifdef CERES_NO_THREADS
-    bool abort = false;
-    for (int i = 0; i < num_residual_blocks; ++i) {
-#endif // CERES_NO_THREADS
-
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
+    // This bool is used to disable the loop if an error is encountered without
+    // breaking out of it. The remaining loop iterations are still run, but with
+    // an empty body, and so will finish quickly.
     std::atomic_bool abort(false);
-
-    ParallelFor(options_.context,
-                0,
-                num_residual_blocks,
-                options_.num_threads,
-                [&](int thread_id, int i) {
-#endif // defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-
-      if (abort) {
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-        return;
-#else
-        continue;
-#endif // defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-      }
-
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-      const ScopedThreadToken scoped_thread_token(&thread_token_provider);
-      const int thread_id = scoped_thread_token.token();
-#endif  // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-
-      EvaluatePreparer* preparer = &evaluate_preparers_[thread_id];
-      EvaluateScratch* scratch = &evaluate_scratch_[thread_id];
-
-      // Prepare block residuals if requested.
-      const ResidualBlock* residual_block = program_->residual_blocks()[i];
-      double* block_residuals = NULL;
-      if (residuals != NULL) {
-        block_residuals = residuals + residual_layout_[i];
-      } else if (gradient != NULL) {
-        block_residuals = scratch->residual_block_residuals.get();
-      }
-
-      // Prepare block jacobians if requested.
-      double** block_jacobians = NULL;
-      if (jacobian != NULL || gradient != NULL) {
-        preparer->Prepare(residual_block,
-                          i,
-                          jacobian,
-                          scratch->jacobian_block_ptrs.get());
-        block_jacobians = scratch->jacobian_block_ptrs.get();
-      }
-
-      // Evaluate the cost, residuals, and jacobians.
-      double block_cost;
-      if (!residual_block->Evaluate(
-              evaluate_options.apply_loss_function,
-              &block_cost,
-              block_residuals,
-              block_jacobians,
-              scratch->residual_block_evaluate_scratch.get())) {
-        abort = true;
-#ifdef CERES_USE_OPENMP
-// This ensures that the OpenMP threads have a consistent view of 'abort'. Do
-// the flush inside the failure case so that there is usually only one
-// synchronization point per loop iteration instead of two.
-#pragma omp flush(abort)
-#endif // CERES_USE_OPENMP
-
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-        return;
-#else
-        continue;
-#endif // defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-      }
-
-      scratch->cost += block_cost;
-
-      // Store the jacobians, if they were requested.
-      if (jacobian != NULL) {
-        jacobian_writer_.Write(i,
-                               residual_layout_[i],
-                               block_jacobians,
-                               jacobian);
-      }
-
-      // Compute and store the gradient, if it was requested.
-      if (gradient != NULL) {
-        int num_residuals = residual_block->NumResiduals();
-        int num_parameter_blocks = residual_block->NumParameterBlocks();
-        for (int j = 0; j < num_parameter_blocks; ++j) {
-          const ParameterBlock* parameter_block =
-              residual_block->parameter_blocks()[j];
-          if (parameter_block->IsConstant()) {
-            continue;
+    ParallelFor(
+        options_.context,
+        0,
+        num_residual_blocks,
+        options_.num_threads,
+        [&](int thread_id, int i) {
+          if (abort) {
+            return;
           }
 
-          MatrixTransposeVectorMultiply<Eigen::Dynamic, Eigen::Dynamic, 1>(
-              block_jacobians[j],
-              num_residuals,
-              parameter_block->LocalSize(),
-              block_residuals,
-              scratch->gradient.get() + parameter_block->delta_offset());
-        }
-      }
-    }
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-    );
-#endif // defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
+          EvaluatePreparer* preparer = &evaluate_preparers_[thread_id];
+          EvaluateScratch* scratch = &evaluate_scratch_[thread_id];
+
+          // Prepare block residuals if requested.
+          const ResidualBlock* residual_block = program_->residual_blocks()[i];
+          double* block_residuals = NULL;
+          if (residuals != NULL) {
+            block_residuals = residuals + residual_layout_[i];
+          } else if (gradient != NULL) {
+            block_residuals = scratch->residual_block_residuals.get();
+          }
+
+          // Prepare block jacobians if requested.
+          double** block_jacobians = NULL;
+          if (jacobian != NULL || gradient != NULL) {
+            preparer->Prepare(residual_block,
+                              i,
+                              jacobian,
+                              scratch->jacobian_block_ptrs.get());
+            block_jacobians = scratch->jacobian_block_ptrs.get();
+          }
+
+          // Evaluate the cost, residuals, and jacobians.
+          double block_cost;
+          if (!residual_block->Evaluate(
+                  evaluate_options.apply_loss_function,
+                  &block_cost,
+                  block_residuals,
+                  block_jacobians,
+                  scratch->residual_block_evaluate_scratch.get())) {
+            abort = true;
+            return;
+          }
+
+          scratch->cost += block_cost;
+
+          // Store the jacobians, if they were requested.
+          if (jacobian != NULL) {
+            jacobian_writer_.Write(i,
+                                   residual_layout_[i],
+                                   block_jacobians,
+                                   jacobian);
+          }
+
+          // Compute and store the gradient, if it was requested.
+          if (gradient != NULL) {
+            int num_residuals = residual_block->NumResiduals();
+            int num_parameter_blocks = residual_block->NumParameterBlocks();
+            for (int j = 0; j < num_parameter_blocks; ++j) {
+              const ParameterBlock* parameter_block =
+                  residual_block->parameter_blocks()[j];
+              if (parameter_block->IsConstant()) {
+                continue;
+              }
+
+              MatrixTransposeVectorMultiply<Eigen::Dynamic, Eigen::Dynamic, 1>(
+                  block_jacobians[j],
+                  num_residuals,
+                  parameter_block->LocalSize(),
+                  block_residuals,
+                  scratch->gradient.get() + parameter_block->delta_offset());
+            }
+          }
+        });
 
     if (!abort) {
       const int num_parameters = program_->NumEffectiveParameters();
diff --git a/internal/ceres/schur_eliminator_impl.h b/internal/ceres/schur_eliminator_impl.h
index da5d922..203dcc9 100644
--- a/internal/ceres/schur_eliminator_impl.h
+++ b/internal/ceres/schur_eliminator_impl.h
@@ -51,6 +51,7 @@
 #include <algorithm>
 #include <map>
 
+#include "Eigen/Dense"
 #include "ceres/block_random_access_matrix.h"
 #include "ceres/block_sparse_matrix.h"
 #include "ceres/block_structure.h"
@@ -58,18 +59,14 @@
 #include "ceres/internal/fixed_array.h"
 #include "ceres/invert_psd_matrix.h"
 #include "ceres/map_util.h"
+#include "ceres/parallel_for.h"
 #include "ceres/schur_eliminator.h"
 #include "ceres/scoped_thread_token.h"
 #include "ceres/small_blas.h"
 #include "ceres/stl_util.h"
 #include "ceres/thread_token_provider.h"
-#include "Eigen/Dense"
 #include "glog/logging.h"
 
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-#include "ceres/parallel_for.h"
-#endif
-
 namespace ceres {
 namespace internal {
 
@@ -190,43 +187,29 @@
 
   // Add the diagonal to the schur complement.
   if (D != NULL) {
-#ifdef CERES_USE_OPENMP
-#pragma omp parallel for num_threads(num_threads_) schedule(dynamic)
-#endif // CERES_USE_OPENMP
+    ParallelFor(
+        context_,
+        num_eliminate_blocks_,
+        num_col_blocks,
+        num_threads_,
+        [&](int i) {
+          const int block_id = i - num_eliminate_blocks_;
+          int r, c, row_stride, col_stride;
+          CellInfo* cell_info = lhs->GetCell(block_id, block_id, &r, &c,
+                                             &row_stride, &col_stride);
+          if (cell_info != NULL) {
+            const int block_size = bs->cols[i].size;
+            typename EigenTypes<Eigen::Dynamic>::ConstVectorRef diag(
+                D + bs->cols[i].position, block_size);
 
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-    for (int i = num_eliminate_blocks_; i < num_col_blocks; ++i) {
-#else
-    ParallelFor(context_, num_eliminate_blocks_, num_col_blocks, num_threads_,
-                [&](int i) {
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-
-      const int block_id = i - num_eliminate_blocks_;
-      int r, c, row_stride, col_stride;
-      CellInfo* cell_info = lhs->GetCell(block_id, block_id,
-                                         &r, &c,
-                                         &row_stride, &col_stride);
-      if (cell_info != NULL) {
-        const int block_size = bs->cols[i].size;
-        typename EigenTypes<Eigen::Dynamic>::ConstVectorRef
-            diag(D + bs->cols[i].position, block_size);
-
-        std::lock_guard<std::mutex> l(cell_info->m);
-        MatrixRef m(cell_info->values, row_stride, col_stride);
-        m.block(r, c, block_size, block_size).diagonal()
-            += diag.array().square().matrix();
-      }
-    }
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-    );
-#endif // defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
+            std::lock_guard<std::mutex> l(cell_info->m);
+            MatrixRef m(cell_info->values, row_stride, col_stride);
+            m.block(r, c, block_size, block_size).diagonal() +=
+                diag.array().square().matrix();
+          }
+        });
   }
 
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-  ThreadTokenProvider thread_token_provider(num_threads_);
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-
-#ifdef CERES_USE_OPENMP
   // Eliminate y blocks one chunk at a time.  For each chunk, compute
   // the entries of the normal equations and the gradient vector block
   // corresponding to the y block and then apply Gaussian elimination
@@ -240,88 +223,76 @@
   // z blocks that share a row block/residual term with the y
   // block. EliminateRowOuterProduct does the corresponding operation
   // for the lhs of the reduced linear system.
-#pragma omp parallel for num_threads(num_threads_) schedule(dynamic)
-#endif // CERES_USE_OPENMP
+  ParallelFor(
+      context_,
+      0,
+      int(chunks_.size()),
+      num_threads_,
+      [&](int thread_id, int i) {
+        double* buffer = buffer_.get() + thread_id * buffer_size_;
+        const Chunk& chunk = chunks_[i];
+        const int e_block_id = bs->rows[chunk.start].cells.front().block_id;
+        const int e_block_size = bs->cols[e_block_id].size;
 
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-  for (int i = 0; i < chunks_.size(); ++i) {
-    const ScopedThreadToken scoped_thread_token(&thread_token_provider);
-    const int thread_id = scoped_thread_token.token();
-#else
-    ParallelFor(context_,
-                0,
-                int(chunks_.size()),
-                num_threads_,
-                [&](int thread_id, int i) {
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
+        VectorRef(buffer, buffer_size_).setZero();
 
-    double* buffer = buffer_.get() + thread_id * buffer_size_;
-    const Chunk& chunk = chunks_[i];
-    const int e_block_id = bs->rows[chunk.start].cells.front().block_id;
-    const int e_block_size = bs->cols[e_block_id].size;
+        typename EigenTypes<kEBlockSize, kEBlockSize>::Matrix
+            ete(e_block_size, e_block_size);
 
-    VectorRef(buffer, buffer_size_).setZero();
+        if (D != NULL) {
+          const typename EigenTypes<kEBlockSize>::ConstVectorRef
+              diag(D + bs->cols[e_block_id].position, e_block_size);
+          ete = diag.array().square().matrix().asDiagonal();
+        } else {
+          ete.setZero();
+        }
 
-    typename EigenTypes<kEBlockSize, kEBlockSize>::Matrix
-        ete(e_block_size, e_block_size);
+        FixedArray<double, 8> g(e_block_size);
+        typename EigenTypes<kEBlockSize>::VectorRef gref(g.get(), e_block_size);
+        gref.setZero();
 
-    if (D != NULL) {
-      const typename EigenTypes<kEBlockSize>::ConstVectorRef
-          diag(D + bs->cols[e_block_id].position, e_block_size);
-      ete = diag.array().square().matrix().asDiagonal();
-    } else {
-      ete.setZero();
-    }
+        // We are going to be computing
+        //
+        //   S += F'F - F'E(E'E)^{-1}E'F
+        //
+        // for each Chunk. The computation is broken down into a number of
+        // function calls as below.
 
-    FixedArray<double, 8> g(e_block_size);
-    typename EigenTypes<kEBlockSize>::VectorRef gref(g.get(), e_block_size);
-    gref.setZero();
+        // Compute the outer product of the e_blocks with themselves (ete
+        // = E'E). Compute the product of the e_blocks with the
+        // corresonding f_blocks (buffer = E'F), the gradient of the terms
+        // in this chunk (g) and add the outer product of the f_blocks to
+        // Schur complement (S += F'F).
+        ChunkDiagonalBlockAndGradient(
+            chunk, A, b, chunk.start, &ete, g.get(), buffer, lhs);
 
-    // We are going to be computing
-    //
-    //   S += F'F - F'E(E'E)^{-1}E'F
-    //
-    // for each Chunk. The computation is broken down into a number of
-    // function calls as below.
+        // Normally one wouldn't compute the inverse explicitly, but
+        // e_block_size will typically be a small number like 3, in
+        // which case its much faster to compute the inverse once and
+        // use it to multiply other matrices/vectors instead of doing a
+        // Solve call over and over again.
+        typename EigenTypes<kEBlockSize, kEBlockSize>::Matrix inverse_ete =
+            InvertPSDMatrix<kEBlockSize>(assume_full_rank_ete_, ete);
 
-    // Compute the outer product of the e_blocks with themselves (ete
-    // = E'E). Compute the product of the e_blocks with the
-    // corresonding f_blocks (buffer = E'F), the gradient of the terms
-    // in this chunk (g) and add the outer product of the f_blocks to
-    // Schur complement (S += F'F).
-    ChunkDiagonalBlockAndGradient(
-        chunk, A, b, chunk.start, &ete, g.get(), buffer, lhs);
+        // For the current chunk compute and update the rhs of the reduced
+        // linear system.
+        //
+        //   rhs = F'b - F'E(E'E)^(-1) E'b
 
-    // Normally one wouldn't compute the inverse explicitly, but
-    // e_block_size will typically be a small number like 3, in
-    // which case its much faster to compute the inverse once and
-    // use it to multiply other matrices/vectors instead of doing a
-    // Solve call over and over again.
-    typename EigenTypes<kEBlockSize, kEBlockSize>::Matrix inverse_ete =
-        InvertPSDMatrix<kEBlockSize>(assume_full_rank_ete_, ete);
+        FixedArray<double, 8> inverse_ete_g(e_block_size);
+        MatrixVectorMultiply<kEBlockSize, kEBlockSize, 0>(
+            inverse_ete.data(),
+            e_block_size,
+            e_block_size,
+            g.get(),
+            inverse_ete_g.get());
 
-    // For the current chunk compute and update the rhs of the reduced
-    // linear system.
-    //
-    //   rhs = F'b - F'E(E'E)^(-1) E'b
+        UpdateRhs(chunk, A, b, chunk.start, inverse_ete_g.get(), rhs);
 
-    FixedArray<double, 8> inverse_ete_g(e_block_size);
-    MatrixVectorMultiply<kEBlockSize, kEBlockSize, 0>(
-        inverse_ete.data(),
-        e_block_size,
-        e_block_size,
-        g.get(),
-        inverse_ete_g.get());
-
-    UpdateRhs(chunk, A, b, chunk.start, inverse_ete_g.get(), rhs);
-
-    // S -= F'E(E'E)^{-1}E'F
-    ChunkOuterProduct(
+        // S -= F'E(E'E)^{-1}E'F
+        ChunkOuterProduct(
         thread_id, bs, inverse_ete, buffer, chunk.buffer_layout, lhs);
-  }
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-  );
-#endif // defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
+      });
 
   // For rows with no e_blocks, the schur complement update reduces to
   // S += F'F.
@@ -338,21 +309,17 @@
                double* y) {
   const CompressedRowBlockStructure* bs = A->block_structure();
 
-#ifdef CERES_USE_OPENMP
-#pragma omp parallel for num_threads(num_threads_) schedule(dynamic)
-#endif // CERES_USE_OPENMP
-
-#if !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-  for (int i = 0; i < chunks_.size(); ++i) {
-#else
-  ParallelFor(context_, 0, int(chunks_.size()), num_threads_, [&](int i) {
-#endif // !(defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS))
-
+  ParallelFor(
+      context_,
+      0,
+      int(chunks_.size()),
+      num_threads_,
+      [&](int i) {
     const Chunk& chunk = chunks_[i];
     const int e_block_id = bs->rows[chunk.start].cells.front().block_id;
     const int e_block_size = bs->cols[e_block_id].size;
 
-    double* y_ptr = y +  bs->cols[e_block_id].position;
+    double* y_ptr = y + bs->cols[e_block_id].position;
     typename EigenTypes<kEBlockSize>::VectorRef y_block(y_ptr, e_block_size);
 
     typename EigenTypes<kEBlockSize, kEBlockSize>::Matrix
@@ -395,17 +362,14 @@
 
       MatrixTransposeMatrixMultiply
           <kRowBlockSize, kEBlockSize, kRowBlockSize, kEBlockSize, 1>(
-              values + e_cell.position, row.block.size, e_block_size,
-              values + e_cell.position, row.block.size, e_block_size,
-              ete.data(), 0, 0, e_block_size, e_block_size);
+          values + e_cell.position, row.block.size, e_block_size,
+          values + e_cell.position, row.block.size, e_block_size,
+          ete.data(), 0, 0, e_block_size, e_block_size);
     }
 
-    y_block = InvertPSDMatrix<kEBlockSize>(assume_full_rank_ete_, ete)
-        * y_block;
-  }
-#if defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
-  );
-#endif // defined(CERES_USE_TBB) || defined(CERES_USE_CXX11_THREADS)
+    y_block =
+        InvertPSDMatrix<kEBlockSize>(assume_full_rank_ete_, ete) * y_block;
+  });
 }
 
 // Update the rhs of the reduced linear system. Compute
diff --git a/jni/Android.mk b/jni/Android.mk
index 2c58147..80861dd 100644
--- a/jni/Android.mk
+++ b/jni/Android.mk
@@ -178,6 +178,7 @@
                    $(CERES_SRC_PATH)/low_rank_inverse_hessian.cc \
                    $(CERES_SRC_PATH)/minimizer.cc \
                    $(CERES_SRC_PATH)/normal_prior.cc \
+                   $(CERES_SRC_PATH)/parallel_for_openmp.cc \
                    $(CERES_SRC_PATH)/parameter_block_ordering.cc \
                    $(CERES_SRC_PATH)/partitioned_matrix_view.cc \
                    $(CERES_SRC_PATH)/polynomial.cc \