diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 4a85e8b..06c0914 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -603,9 +603,6 @@
   add_executable(spmv_benchmark spmv_benchmark.cc)
   add_dependencies_to_benchmark(spmv_benchmark)
 
-  add_executable(partitioned_matrix_view_benchmark partitioned_matrix_view_benchmark.cc)
-  add_dependencies_to_benchmark(partitioned_matrix_view_benchmark)
-
   add_executable(block_jacobi_preconditioner_benchmark
     block_jacobi_preconditioner_benchmark.cc)
   add_dependencies_to_benchmark(block_jacobi_preconditioner_benchmark)
diff --git a/internal/ceres/evaluation_benchmark.cc b/internal/ceres/evaluation_benchmark.cc
index 2f9fff0..9849b5a 100644
--- a/internal/ceres/evaluation_benchmark.cc
+++ b/internal/ceres/evaluation_benchmark.cc
@@ -145,23 +145,20 @@
     return crs_jacobian.get();
   }
 
-  const PartitionedView* PartitionedMatrixViewJacobian(ContextImpl* context) {
-    if (!partitioned_view_jacobian) {
-      auto block_sparse = BlockSparseJacobian(context);
-      partitioned_view_jacobian = std::make_unique<PartitionedView>(
-          *block_sparse, bal_problem->num_points());
-    }
+  const PartitionedView* PartitionedMatrixViewJacobian(
+      const LinearSolver::Options& options) {
+    auto block_sparse = BlockSparseJacobian(options.context);
+    partitioned_view_jacobian =
+        std::make_unique<PartitionedView>(options, *block_sparse);
     return partitioned_view_jacobian.get();
   }
 
   const PartitionedView* PartitionedMatrixViewJacobianWithTranspose(
-      ContextImpl* context) {
-    if (!partitioned_view_jacobian_with_transpose) {
-      auto block_sparse_transpose = BlockSparseJacobianWithTranspose(context);
-      partitioned_view_jacobian_with_transpose =
-          std::make_unique<PartitionedView>(*block_sparse_transpose,
-                                            bal_problem->num_points());
-    }
+      const LinearSolver::Options& options) {
+    auto block_sparse_transpose =
+        BlockSparseJacobianWithTranspose(options.context);
+    partitioned_view_jacobian_with_transpose =
+        std::make_unique<PartitionedView>(options, *block_sparse_transpose);
     return partitioned_view_jacobian_with_transpose.get();
   }
 
@@ -243,16 +240,17 @@
 static void PMVRightMultiplyAndAccumulateF(benchmark::State& state,
                                            BALData* data,
                                            ContextImpl* context) {
-  const int num_threads = state.range(0);
-
-  auto jacobian = data->PartitionedMatrixViewJacobian(context);
+  LinearSolver::Options options;
+  options.num_threads = state.range(0);
+  options.elimination_groups.push_back(data->bal_problem->num_points());
+  options.context = context;
+  auto jacobian = data->PartitionedMatrixViewJacobian(options);
 
   Vector y = Vector::Zero(jacobian->num_rows());
   Vector x = Vector::Random(jacobian->num_cols_f());
 
   for (auto _ : state) {
-    jacobian->RightMultiplyAndAccumulateF(
-        x.data(), y.data(), context, num_threads);
+    jacobian->RightMultiplyAndAccumulateF(x.data(), y.data());
   }
   CHECK_GT(y.squaredNorm(), 0.);
 }
@@ -260,49 +258,35 @@
 static void PMVLeftMultiplyAndAccumulateF(benchmark::State& state,
                                           BALData* data,
                                           ContextImpl* context) {
-  const int num_threads = state.range(0);
-
-  auto jacobian = data->PartitionedMatrixViewJacobian(context);
+  LinearSolver::Options options;
+  options.num_threads = state.range(0);
+  options.elimination_groups.push_back(data->bal_problem->num_points());
+  options.context = context;
+  auto jacobian = data->PartitionedMatrixViewJacobianWithTranspose(options);
 
   Vector y = Vector::Zero(jacobian->num_cols_f());
   Vector x = Vector::Random(jacobian->num_rows());
 
   for (auto _ : state) {
-    jacobian->LeftMultiplyAndAccumulateF(
-        x.data(), y.data(), context, num_threads);
+    jacobian->LeftMultiplyAndAccumulateF(x.data(), y.data());
   }
   CHECK_GT(y.squaredNorm(), 0.);
 }
 
-static void PMVLeftMultiplyAndAccumulateWithTransposeF(benchmark::State& state,
-                                                       BALData* data,
-                                                       ContextImpl* context) {
-  const int num_threads = state.range(0);
-
-  auto jacobian = data->PartitionedMatrixViewJacobianWithTranspose(context);
-
-  Vector y = Vector::Zero(jacobian->num_cols_f());
-  Vector x = Vector::Random(jacobian->num_rows());
-
-  for (auto _ : state) {
-    jacobian->LeftMultiplyAndAccumulateF(
-        x.data(), y.data(), context, num_threads);
-  }
-  CHECK_GT(y.squaredNorm(), 0.);
-}
 static void PMVRightMultiplyAndAccumulateE(benchmark::State& state,
                                            BALData* data,
                                            ContextImpl* context) {
-  const int num_threads = state.range(0);
-
-  auto jacobian = data->PartitionedMatrixViewJacobian(context);
+  LinearSolver::Options options;
+  options.num_threads = state.range(0);
+  options.elimination_groups.push_back(data->bal_problem->num_points());
+  options.context = context;
+  auto jacobian = data->PartitionedMatrixViewJacobian(options);
 
   Vector y = Vector::Zero(jacobian->num_rows());
   Vector x = Vector::Random(jacobian->num_cols_e());
 
   for (auto _ : state) {
-    jacobian->RightMultiplyAndAccumulateE(
-        x.data(), y.data(), context, num_threads);
+    jacobian->RightMultiplyAndAccumulateE(x.data(), y.data());
   }
   CHECK_GT(y.squaredNorm(), 0.);
 }
@@ -310,33 +294,17 @@
 static void PMVLeftMultiplyAndAccumulateE(benchmark::State& state,
                                           BALData* data,
                                           ContextImpl* context) {
-  const int num_threads = state.range(0);
-
-  auto jacobian = data->PartitionedMatrixViewJacobian(context);
+  LinearSolver::Options options;
+  options.num_threads = state.range(0);
+  options.elimination_groups.push_back(data->bal_problem->num_points());
+  options.context = context;
+  auto jacobian = data->PartitionedMatrixViewJacobianWithTranspose(options);
 
   Vector y = Vector::Zero(jacobian->num_cols_e());
   Vector x = Vector::Random(jacobian->num_rows());
 
   for (auto _ : state) {
-    jacobian->LeftMultiplyAndAccumulateE(
-        x.data(), y.data(), context, num_threads);
-  }
-  CHECK_GT(y.squaredNorm(), 0.);
-}
-
-static void PMVLeftMultiplyAndAccumulateWithTransposeE(benchmark::State& state,
-                                                       BALData* data,
-                                                       ContextImpl* context) {
-  const int num_threads = state.range(0);
-
-  auto jacobian = data->PartitionedMatrixViewJacobianWithTranspose(context);
-
-  Vector y = Vector::Zero(jacobian->num_cols_e());
-  Vector x = Vector::Random(jacobian->num_rows());
-
-  for (auto _ : state) {
-    jacobian->LeftMultiplyAndAccumulateE(
-        x.data(), y.data(), context, num_threads);
+    jacobian->LeftMultiplyAndAccumulateE(x.data(), y.data());
   }
   CHECK_GT(y.squaredNorm(), 0.);
 }
@@ -363,22 +331,6 @@
                                               ContextImpl* context) {
   const int num_threads = state.range(0);
 
-  auto jacobian = data->BlockSparseJacobian(context);
-
-  Vector y = Vector::Zero(jacobian->num_cols());
-  Vector x = Vector::Random(jacobian->num_rows());
-
-  for (auto _ : state) {
-    jacobian->LeftMultiplyAndAccumulate(
-        x.data(), y.data(), context, num_threads);
-  }
-  CHECK_GT(y.squaredNorm(), 0.);
-}
-
-static void JacobianLeftMultiplyAndAccumulateWithTranspose(
-    benchmark::State& state, BALData* data, ContextImpl* context) {
-  const int num_threads = state.range(0);
-
   auto jacobian = data->BlockSparseJacobianWithTranspose(context);
 
   Vector y = Vector::Zero(jacobian->num_cols());
@@ -551,15 +503,6 @@
         ceres::internal::JacobianLeftMultiplyAndAccumulate,
         data,
         &context)
-        ->Arg(1);
-
-    const std::string name_left_product_transpose =
-        "JacobianLeftMultiplyAndAccumulateWithTranspose<" + path + ">";
-    ::benchmark::RegisterBenchmark(
-        name_left_product_transpose.c_str(),
-        ceres::internal::JacobianLeftMultiplyAndAccumulateWithTranspose,
-        data,
-        &context)
         ->Arg(1)
         ->Arg(2)
         ->Arg(4)
@@ -573,14 +516,6 @@
         ceres::internal::PMVLeftMultiplyAndAccumulateF,
         data,
         &context)
-        ->Arg(1);
-    const std::string name_left_product_partitioned_transpose_f =
-        "PMVLeftMultiplyAndAccumulateWithTransposeF<" + path + ">";
-    ::benchmark::RegisterBenchmark(
-        name_left_product_partitioned_transpose_f.c_str(),
-        ceres::internal::PMVLeftMultiplyAndAccumulateWithTransposeF,
-        data,
-        &context)
         ->Arg(1)
         ->Arg(2)
         ->Arg(4)
@@ -594,15 +529,6 @@
         ceres::internal::PMVLeftMultiplyAndAccumulateE,
         data,
         &context)
-        ->Arg(1);
-
-    const std::string name_left_product_partitioned_transpose_e =
-        "PMVLeftMultiplyAndAccumulateWithTransposeE<" + path + ">";
-    ::benchmark::RegisterBenchmark(
-        name_left_product_partitioned_transpose_e.c_str(),
-        ceres::internal::PMVLeftMultiplyAndAccumulateWithTransposeE,
-        data,
-        &context)
         ->Arg(1)
         ->Arg(2)
         ->Arg(4)
diff --git a/internal/ceres/fake_bundle_adjustment_jacobian.cc b/internal/ceres/fake_bundle_adjustment_jacobian.cc
index af1d03c..d8f3c01 100644
--- a/internal/ceres/fake_bundle_adjustment_jacobian.cc
+++ b/internal/ceres/fake_bundle_adjustment_jacobian.cc
@@ -109,8 +109,10 @@
       PartitionedMatrixView<2, Eigen::Dynamic, Eigen::Dynamic>;
   auto block_sparse_matrix = CreateFakeBundleAdjustmentJacobian(
       num_cameras, num_points, camera_size, landmark_size, visibility, rng);
+  LinearSolver::Options options;
+  options.elimination_groups.push_back(num_points);
   auto partitioned_view =
-      std::make_unique<PartitionedView>(*block_sparse_matrix, num_points);
+      std::make_unique<PartitionedView>(options, *block_sparse_matrix);
   return std::make_pair(std::move(partitioned_view),
                         std::move(block_sparse_matrix));
 }
diff --git a/internal/ceres/implicit_schur_complement.cc b/internal/ceres/implicit_schur_complement.cc
index 6249c8c..4e36048 100644
--- a/internal/ceres/implicit_schur_complement.cc
+++ b/internal/ceres/implicit_schur_complement.cc
@@ -105,15 +105,11 @@
                                                          double* y) const {
   // y1 = F x
   tmp_rows_.setZero();
-  A_->RightMultiplyAndAccumulateF(
-      x, tmp_rows_.data(), options_.context, options_.num_threads);
+  A_->RightMultiplyAndAccumulateF(x, tmp_rows_.data());
 
   // y2 = E' y1
   tmp_e_cols_.setZero();
-  A_->LeftMultiplyAndAccumulateE(tmp_rows_.data(),
-                                 tmp_e_cols_.data(),
-                                 options_.context,
-                                 options_.num_threads);
+  A_->LeftMultiplyAndAccumulateE(tmp_rows_.data(), tmp_e_cols_.data());
 
   // y3 = -(E'E)^-1 y2
   tmp_e_cols_2_.setZero();
@@ -121,13 +117,11 @@
                                                           tmp_e_cols_2_.data(),
                                                           options_.context,
                                                           options_.num_threads);
+
   tmp_e_cols_2_ *= -1.0;
 
   // y1 = y1 + E y3
-  A_->RightMultiplyAndAccumulateE(tmp_e_cols_2_.data(),
-                                  tmp_rows_.data(),
-                                  options_.context,
-                                  options_.num_threads);
+  A_->RightMultiplyAndAccumulateE(tmp_e_cols_2_.data(), tmp_rows_.data());
 
   // y5 = D * x
   if (D_ != nullptr) {
@@ -140,8 +134,7 @@
   }
 
   // y = y5 + F' y1
-  A_->LeftMultiplyAndAccumulateF(
-      tmp_rows_.data(), y, options_.context, options_.num_threads);
+  A_->LeftMultiplyAndAccumulateF(tmp_rows_.data(), y);
 }
 
 void ImplicitSchurComplement::InversePowerSeriesOperatorRightMultiplyAccumulate(
@@ -149,15 +142,11 @@
   CHECK(compute_ftf_inverse_);
   // y1 = F x
   tmp_rows_.setZero();
-  A_->RightMultiplyAndAccumulateF(
-      x, tmp_rows_.data(), options_.context, options_.num_threads);
+  A_->RightMultiplyAndAccumulateF(x, tmp_rows_.data());
 
   // y2 = E' y1
   tmp_e_cols_.setZero();
-  A_->LeftMultiplyAndAccumulateE(tmp_rows_.data(),
-                                 tmp_e_cols_.data(),
-                                 options_.context,
-                                 options_.num_threads);
+  A_->LeftMultiplyAndAccumulateE(tmp_rows_.data(), tmp_e_cols_.data());
 
   // y3 = (E'E)^-1 y2
   tmp_e_cols_2_.setZero();
@@ -167,17 +156,11 @@
                                                           options_.num_threads);
   // y1 = E y3
   tmp_rows_.setZero();
-  A_->RightMultiplyAndAccumulateE(tmp_e_cols_2_.data(),
-                                  tmp_rows_.data(),
-                                  options_.context,
-                                  options_.num_threads);
+  A_->RightMultiplyAndAccumulateE(tmp_e_cols_2_.data(), tmp_rows_.data());
 
   // y4 = F' y1
   tmp_f_cols_.setZero();
-  A_->LeftMultiplyAndAccumulateF(tmp_rows_.data(),
-                                 tmp_f_cols_.data(),
-                                 options_.context,
-                                 options_.num_threads);
+  A_->LeftMultiplyAndAccumulateF(tmp_rows_.data(), tmp_f_cols_.data());
 
   // y += (F'F)^-1 y4
   block_diagonal_FtF_inverse_->RightMultiplyAndAccumulate(
@@ -219,18 +202,14 @@
 
   // y1 = F x
   tmp_rows_.setZero();
-  A_->RightMultiplyAndAccumulateF(
-      x, tmp_rows_.data(), options_.context, options_.num_threads);
+  A_->RightMultiplyAndAccumulateF(x, tmp_rows_.data());
 
   // y2 = b - y1
   tmp_rows_ = ConstVectorRef(b_, num_rows) - tmp_rows_;
 
   // y3 = E' y2
   tmp_e_cols_.setZero();
-  A_->LeftMultiplyAndAccumulateE(tmp_rows_.data(),
-                                 tmp_e_cols_.data(),
-                                 options_.context,
-                                 options_.num_threads);
+  A_->LeftMultiplyAndAccumulateE(tmp_rows_.data(), tmp_e_cols_.data());
 
   // y = (E'E)^-1 y3
   VectorRef(y, num_cols).setZero();
@@ -254,8 +233,7 @@
 void ImplicitSchurComplement::UpdateRhs() {
   // y1 = E'b
   tmp_e_cols_.setZero();
-  A_->LeftMultiplyAndAccumulateE(
-      b_, tmp_e_cols_.data(), options_.context, options_.num_threads);
+  A_->LeftMultiplyAndAccumulateE(b_, tmp_e_cols_.data());
 
   // y2 = (E'E)^-1 y1
   Vector y2 = Vector::Zero(A_->num_cols_e());
@@ -264,16 +242,14 @@
 
   // y3 = E y2
   tmp_rows_.setZero();
-  A_->RightMultiplyAndAccumulateE(
-      y2.data(), tmp_rows_.data(), options_.context, options_.num_threads);
+  A_->RightMultiplyAndAccumulateE(y2.data(), tmp_rows_.data());
 
   // y3 = b - y3
   tmp_rows_ = ConstVectorRef(b_, A_->num_rows()) - tmp_rows_;
 
   // rhs = F' y3
   rhs_.setZero();
-  A_->LeftMultiplyAndAccumulateF(
-      tmp_rows_.data(), rhs_.data(), options_.context, options_.num_threads);
+  A_->LeftMultiplyAndAccumulateF(tmp_rows_.data(), rhs_.data());
 }
 
 }  // namespace ceres::internal
diff --git a/internal/ceres/parallel_for.h b/internal/ceres/parallel_for.h
index 107d19c..fe5d242 100644
--- a/internal/ceres/parallel_for.h
+++ b/internal/ceres/parallel_for.h
@@ -299,6 +299,42 @@
               });
 }
 
+// Execute function for every element in the range [start, end) with at most
+// num_threads, using the user provided partitioning. taking into account
+// user-provided integer cumulative costs of iterations.
+template <typename F>
+void ParallelFor(ContextImpl* context,
+                 int start,
+                 int end,
+                 int num_threads,
+                 const F& function,
+                 const std::vector<int>& partitions) {
+  using namespace parallel_for_details;
+  CHECK_GT(num_threads, 0);
+  if (start >= end) {
+    return;
+  }
+  CHECK_EQ(partitions.front(), start);
+  CHECK_EQ(partitions.back(), end);
+  if (num_threads == 1 || end - start <= num_threads) {
+    ParallelFor(context, start, end, num_threads, function);
+    return;
+  }
+  CHECK_GT(partitions.size(), 1);
+  const int num_partitions = partitions.size() - 1;
+  ParallelFor(context,
+              0,
+              num_partitions,
+              num_threads,
+              [&function, &partitions](int thread_id, int partition_id) {
+                const int partition_start = partitions[partition_id];
+                const int partition_end = partitions[partition_id + 1];
+
+                for (int i = partition_start; i < partition_end; ++i) {
+                  Invoke<F>(thread_id, i, function);
+                }
+              });
+}
 }  // namespace ceres::internal
 
 // Backend-specific implementations of ParallelInvoke
diff --git a/internal/ceres/partitioned_matrix_view.cc b/internal/ceres/partitioned_matrix_view.cc
index 2710366..d952a2a 100644
--- a/internal/ceres/partitioned_matrix_view.cc
+++ b/internal/ceres/partitioned_matrix_view.cc
@@ -55,121 +55,121 @@
      (options.e_block_size == 2) &&
      (options.f_block_size == 2)) {
     return std::make_unique<PartitionedMatrixView<2,2, 2>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 2) &&
      (options.f_block_size == 3)) {
     return std::make_unique<PartitionedMatrixView<2,2, 3>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 2) &&
      (options.f_block_size == 4)) {
     return std::make_unique<PartitionedMatrixView<2,2, 4>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 2)) {
     return std::make_unique<PartitionedMatrixView<2,2, Eigen::Dynamic>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 3) &&
      (options.f_block_size == 3)) {
     return std::make_unique<PartitionedMatrixView<2,3, 3>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 3) &&
      (options.f_block_size == 4)) {
     return std::make_unique<PartitionedMatrixView<2,3, 4>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 3) &&
      (options.f_block_size == 6)) {
     return std::make_unique<PartitionedMatrixView<2,3, 6>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 3) &&
      (options.f_block_size == 9)) {
     return std::make_unique<PartitionedMatrixView<2,3, 9>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 3)) {
     return std::make_unique<PartitionedMatrixView<2,3, Eigen::Dynamic>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 4) &&
      (options.f_block_size == 3)) {
     return std::make_unique<PartitionedMatrixView<2,4, 3>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 4) &&
      (options.f_block_size == 4)) {
     return std::make_unique<PartitionedMatrixView<2,4, 4>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 4) &&
      (options.f_block_size == 6)) {
     return std::make_unique<PartitionedMatrixView<2,4, 6>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 4) &&
      (options.f_block_size == 8)) {
     return std::make_unique<PartitionedMatrixView<2,4, 8>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 4) &&
      (options.f_block_size == 9)) {
     return std::make_unique<PartitionedMatrixView<2,4, 9>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 2) &&
      (options.e_block_size == 4)) {
     return std::make_unique<PartitionedMatrixView<2,4, Eigen::Dynamic>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if (options.row_block_size == 2) {
     return std::make_unique<PartitionedMatrixView<2,Eigen::Dynamic, Eigen::Dynamic>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 3) &&
      (options.e_block_size == 3) &&
      (options.f_block_size == 3)) {
     return std::make_unique<PartitionedMatrixView<3,3, 3>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 4) &&
      (options.e_block_size == 4) &&
      (options.f_block_size == 2)) {
     return std::make_unique<PartitionedMatrixView<4,4, 2>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 4) &&
      (options.e_block_size == 4) &&
      (options.f_block_size == 3)) {
     return std::make_unique<PartitionedMatrixView<4,4, 3>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 4) &&
      (options.e_block_size == 4) &&
      (options.f_block_size == 4)) {
     return std::make_unique<PartitionedMatrixView<4,4, 4>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
   if ((options.row_block_size == 4) &&
      (options.e_block_size == 4)) {
     return std::make_unique<PartitionedMatrixView<4,4, Eigen::Dynamic>>(
-                   matrix, options.elimination_groups[0]);
+                   options, matrix);
   }
 
 #endif
@@ -179,7 +179,7 @@
   return std::make_unique<PartitionedMatrixView<Eigen::Dynamic,
                                                 Eigen::Dynamic,
                                                 Eigen::Dynamic>>(
-      matrix, options.elimination_groups[0]);
+      options, matrix);
 };
 
 }  // namespace ceres::internal
diff --git a/internal/ceres/partitioned_matrix_view.h b/internal/ceres/partitioned_matrix_view.h
index 533a0e5..f34fd4a 100644
--- a/internal/ceres/partitioned_matrix_view.h
+++ b/internal/ceres/partitioned_matrix_view.h
@@ -70,33 +70,25 @@
 
   // y += E'x
   virtual void LeftMultiplyAndAccumulateE(const double* x, double* y) const = 0;
-  virtual void LeftMultiplyAndAccumulateE(const double* x,
-                                          double* y,
-                                          ContextImpl* context,
-                                          int num_threads) const = 0;
+  virtual void LeftMultiplyAndAccumulateESingleThreaded(const double* x,
+                                                        double* y) const = 0;
+  virtual void LeftMultiplyAndAccumulateEMultiThreaded(const double* x,
+                                                       double* y) const = 0;
 
   // y += F'x
   virtual void LeftMultiplyAndAccumulateF(const double* x, double* y) const = 0;
-  virtual void LeftMultiplyAndAccumulateF(const double* x,
-                                          double* y,
-                                          ContextImpl* context,
-                                          int num_threads) const = 0;
+  virtual void LeftMultiplyAndAccumulateFSingleThreaded(const double* x,
+                                                        double* y) const = 0;
+  virtual void LeftMultiplyAndAccumulateFMultiThreaded(const double* x,
+                                                       double* y) const = 0;
 
   // y += Ex
   virtual void RightMultiplyAndAccumulateE(const double* x,
                                            double* y) const = 0;
-  virtual void RightMultiplyAndAccumulateE(const double* x,
-                                           double* y,
-                                           ContextImpl* context,
-                                           int num_threads) const = 0;
 
   // y += Fx
   virtual void RightMultiplyAndAccumulateF(const double* x,
                                            double* y) const = 0;
-  virtual void RightMultiplyAndAccumulateF(const double* x,
-                                           double* y,
-                                           ContextImpl* context,
-                                           int num_threads) const = 0;
 
   // Create and return the block diagonal of the matrix E'E.
   virtual std::unique_ptr<BlockSparseMatrix> CreateBlockDiagonalEtE() const = 0;
@@ -128,6 +120,8 @@
   virtual int num_cols_f()       const = 0;
   virtual int num_rows()         const = 0;
   virtual int num_cols()         const = 0;
+  virtual const std::vector<int>& e_cols_partition() const = 0;
+  virtual const std::vector<int>& f_cols_partition() const = 0;
   // clang-format on
 
   static std::unique_ptr<PartitionedMatrixViewBase> Create(
@@ -141,29 +135,34 @@
     : public PartitionedMatrixViewBase {
  public:
   // matrix = [E F], where the matrix E contains the first
-  // num_col_blocks_a column blocks.
-  PartitionedMatrixView(const BlockSparseMatrix& matrix, int num_col_blocks_e);
+  // options.elimination_groups[0] column blocks.
+  PartitionedMatrixView(const LinearSolver::Options& options,
+                        const BlockSparseMatrix& matrix);
 
-  void LeftMultiplyAndAccumulateE(const double* x, double* y) const final;
-  void LeftMultiplyAndAccumulateF(const double* x, double* y) const final;
-  void LeftMultiplyAndAccumulateE(const double* x,
-                                  double* y,
-                                  ContextImpl* context,
-                                  int num_threads) const final;
-  void LeftMultiplyAndAccumulateF(const double* x,
-                                  double* y,
-                                  ContextImpl* context,
-                                  int num_threads) const final;
-  void RightMultiplyAndAccumulateE(const double* x, double* y) const final;
-  void RightMultiplyAndAccumulateF(const double* x, double* y) const final;
-  void RightMultiplyAndAccumulateE(const double* x,
-                                   double* y,
-                                   ContextImpl* context,
-                                   int num_threads) const final;
-  void RightMultiplyAndAccumulateF(const double* x,
-                                   double* y,
-                                   ContextImpl* context,
-                                   int num_threads) const final;
+  // y += E'x
+  virtual void LeftMultiplyAndAccumulateE(const double* x,
+                                          double* y) const final;
+  virtual void LeftMultiplyAndAccumulateESingleThreaded(const double* x,
+                                                        double* y) const final;
+  virtual void LeftMultiplyAndAccumulateEMultiThreaded(const double* x,
+                                                       double* y) const final;
+
+  // y += F'x
+  virtual void LeftMultiplyAndAccumulateF(const double* x,
+                                          double* y) const final;
+  virtual void LeftMultiplyAndAccumulateFSingleThreaded(const double* x,
+                                                        double* y) const final;
+  virtual void LeftMultiplyAndAccumulateFMultiThreaded(const double* x,
+                                                       double* y) const final;
+
+  // y += Ex
+  virtual void RightMultiplyAndAccumulateE(const double* x,
+                                           double* y) const final;
+
+  // y += Fx
+  virtual void RightMultiplyAndAccumulateF(const double* x,
+                                           double* y) const final;
+
   std::unique_ptr<BlockSparseMatrix> CreateBlockDiagonalEtE() const final;
   std::unique_ptr<BlockSparseMatrix> CreateBlockDiagonalFtF() const final;
   void UpdateBlockDiagonalEtE(BlockSparseMatrix* block_diagonal) const final;
@@ -176,17 +175,26 @@
   int num_rows()         const final { return matrix_.num_rows(); }
   int num_cols()         const final { return matrix_.num_cols(); }
   // clang-format on
+  const std::vector<int>& e_cols_partition() const final {
+    return e_cols_partition_;
+  }
+  const std::vector<int>& f_cols_partition() const final {
+    return f_cols_partition_;
+  }
 
  private:
   std::unique_ptr<BlockSparseMatrix> CreateBlockDiagonalMatrixLayout(
       int start_col_block, int end_col_block) const;
 
+  const LinearSolver::Options options_;
   const BlockSparseMatrix& matrix_;
   int num_row_blocks_e_;
   int num_col_blocks_e_;
   int num_col_blocks_f_;
   int num_cols_e_;
   int num_cols_f_;
+  std::vector<int> e_cols_partition_;
+  std::vector<int> f_cols_partition_;
 };
 
 }  // namespace ceres::internal
diff --git a/internal/ceres/partitioned_matrix_view_benchmark.cc b/internal/ceres/partitioned_matrix_view_benchmark.cc
deleted file mode 100644
index 63d0e34..0000000
--- a/internal/ceres/partitioned_matrix_view_benchmark.cc
+++ /dev/null
@@ -1,147 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2022 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.
-
-#include <memory>
-#include <random>
-
-#include "benchmark/benchmark.h"
-#include "ceres/context_impl.h"
-#include "ceres/fake_bundle_adjustment_jacobian.h"
-#include "ceres/partitioned_matrix_view.h"
-
-constexpr int kNumCameras = 1000;
-constexpr int kNumPoints = 10000;
-constexpr int kCameraSize = 6;
-constexpr int kPointSize = 3;
-constexpr double kVisibility = 0.1;
-
-namespace ceres::internal {
-
-static void BM_PatitionedViewRightMultiplyAndAccumulateE_Static(
-    benchmark::State& state) {
-  const int num_threads = state.range(0);
-  std::mt19937 prng;
-  auto [partitioned_view, jacobian] =
-      CreateFakeBundleAdjustmentPartitionedJacobian<kPointSize, kCameraSize>(
-          kNumCameras, kNumPoints, kVisibility, prng);
-
-  ContextImpl context;
-  context.EnsureMinimumThreads(num_threads);
-
-  Vector x(partitioned_view->num_cols_e());
-  Vector y(partitioned_view->num_rows());
-  x.setRandom();
-  y.setRandom();
-  double sum = 0;
-  for (auto _ : state) {
-    partitioned_view->RightMultiplyAndAccumulateE(
-        x.data(), y.data(), &context, num_threads);
-    sum += y.norm();
-  }
-  CHECK_NE(sum, 0.0);
-}
-BENCHMARK(BM_PatitionedViewRightMultiplyAndAccumulateE_Static)
-    ->Arg(1)
-    ->Arg(2)
-    ->Arg(4)
-    ->Arg(8)
-    ->Arg(16);
-
-static void BM_PatitionedViewRightMultiplyAndAccumulateE_Dynamic(
-    benchmark::State& state) {
-  std::mt19937 prng;
-  auto [partitioned_view, jacobian] =
-      CreateFakeBundleAdjustmentPartitionedJacobian(
-          kNumCameras, kNumPoints, kCameraSize, kPointSize, kVisibility, prng);
-
-  Vector x(partitioned_view->num_cols_e());
-  Vector y(partitioned_view->num_rows());
-  x.setRandom();
-  y.setRandom();
-  double sum = 0;
-  for (auto _ : state) {
-    partitioned_view->RightMultiplyAndAccumulateE(x.data(), y.data());
-    sum += y.norm();
-  }
-  CHECK_NE(sum, 0.0);
-}
-BENCHMARK(BM_PatitionedViewRightMultiplyAndAccumulateE_Dynamic);
-
-static void BM_PatitionedViewRightMultiplyAndAccumulateF_Static(
-    benchmark::State& state) {
-  const int num_threads = state.range(0);
-  std::mt19937 prng;
-  auto [partitioned_view, jacobian] =
-      CreateFakeBundleAdjustmentPartitionedJacobian<kPointSize, kCameraSize>(
-          kNumCameras, kNumPoints, kVisibility, prng);
-
-  ContextImpl context;
-  context.EnsureMinimumThreads(num_threads);
-
-  Vector x(partitioned_view->num_cols_f());
-  Vector y(partitioned_view->num_rows());
-  x.setRandom();
-  y.setRandom();
-  double sum = 0;
-  for (auto _ : state) {
-    partitioned_view->RightMultiplyAndAccumulateF(
-        x.data(), y.data(), &context, num_threads);
-    sum += y.norm();
-  }
-  CHECK_NE(sum, 0.0);
-}
-BENCHMARK(BM_PatitionedViewRightMultiplyAndAccumulateF_Static)
-    ->Arg(1)
-    ->Arg(2)
-    ->Arg(4)
-    ->Arg(8)
-    ->Arg(16);
-
-static void BM_PatitionedViewRightMultiplyAndAccumulateF_Dynamic(
-    benchmark::State& state) {
-  std::mt19937 prng;
-  auto [partitioned_view, jacobian] =
-      CreateFakeBundleAdjustmentPartitionedJacobian(
-          kNumCameras, kNumPoints, kCameraSize, kPointSize, kVisibility, prng);
-
-  Vector x(partitioned_view->num_cols_f());
-  Vector y(partitioned_view->num_rows());
-  x.setRandom();
-  y.setRandom();
-  double sum = 0;
-  for (auto _ : state) {
-    partitioned_view->RightMultiplyAndAccumulateF(x.data(), y.data());
-    sum += y.norm();
-  }
-  CHECK_NE(sum, 0.0);
-}
-BENCHMARK(BM_PatitionedViewRightMultiplyAndAccumulateF_Dynamic);
-
-}  // namespace ceres::internal
-
-BENCHMARK_MAIN();
diff --git a/internal/ceres/partitioned_matrix_view_impl.h b/internal/ceres/partitioned_matrix_view_impl.h
index 689f09e..1ca2e14 100644
--- a/internal/ceres/partitioned_matrix_view_impl.h
+++ b/internal/ceres/partitioned_matrix_view_impl.h
@@ -45,11 +45,14 @@
 
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
-    PartitionedMatrixView(const BlockSparseMatrix& matrix, int num_col_blocks_e)
-    : matrix_(matrix), num_col_blocks_e_(num_col_blocks_e) {
+    PartitionedMatrixView(const LinearSolver::Options& options,
+                          const BlockSparseMatrix& matrix)
+
+    : options_(options), matrix_(matrix) {
   const CompressedRowBlockStructure* bs = matrix_.block_structure();
   CHECK(bs != nullptr);
 
+  num_col_blocks_e_ = options_.elimination_groups[0];
   num_col_blocks_f_ = bs->cols.size() - num_col_blocks_e_;
 
   // Compute the number of row blocks in E. The number of row blocks
@@ -79,6 +82,25 @@
   }
 
   CHECK_EQ(num_cols_e_ + num_cols_f_, matrix_.num_cols());
+
+  auto transpose_bs = matrix_.transpose_block_structure();
+  const int num_threads = options_.num_threads;
+  if (transpose_bs != nullptr && num_threads > 1) {
+    int kMaxPartitions = num_threads * 4;
+    e_cols_partition_ = parallel_for_details::ComputePartition(
+        0,
+        num_col_blocks_e_,
+        kMaxPartitions,
+        transpose_bs->rows.data(),
+        [](const CompressedRow& row) { return row.cumulative_nnz; });
+
+    f_cols_partition_ = parallel_for_details::ComputePartition(
+        num_col_blocks_e_,
+        num_col_blocks_e_ + num_col_blocks_f_,
+        kMaxPartitions,
+        transpose_bs->rows.data(),
+        [](const CompressedRow& row) { return row.cumulative_nnz; });
+  }
 }
 
 // The next four methods don't seem to be particularly cache
@@ -89,23 +111,14 @@
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
     RightMultiplyAndAccumulateE(const double* x, double* y) const {
-  RightMultiplyAndAccumulateE(x, y, nullptr, 1);
-}
-
-template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
-void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
-    RightMultiplyAndAccumulateE(const double* x,
-                                double* y,
-                                ContextImpl* context,
-                                int num_threads) const {
   // Iterate over the first num_row_blocks_e_ row blocks, and multiply
   // by the first cell in each row block.
   auto bs = matrix_.block_structure();
   const double* values = matrix_.values();
-  ParallelFor(context,
+  ParallelFor(options_.context,
               0,
               num_row_blocks_e_,
-              num_threads,
+              options_.num_threads,
               [values, bs, x, y](int row_block_id) {
                 const Cell& cell = bs->rows[row_block_id].cells[0];
                 const int row_block_pos = bs->rows[row_block_id].block.position;
@@ -125,15 +138,6 @@
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
     RightMultiplyAndAccumulateF(const double* x, double* y) const {
-  RightMultiplyAndAccumulateF(x, y, nullptr, 1);
-}
-
-template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
-void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
-    RightMultiplyAndAccumulateF(const double* x,
-                                double* y,
-                                ContextImpl* context,
-                                int num_threads) const {
   // Iterate over row blocks, and if the row block is in E, then
   // multiply by all the cells except the first one which is of type
   // E. If the row block is not in E (i.e its in the bottom
@@ -143,10 +147,10 @@
   const int num_row_blocks = bs->rows.size();
   const int num_cols_e = num_cols_e_;
   const double* values = matrix_.values();
-  ParallelFor(context,
+  ParallelFor(options_.context,
               0,
               num_row_blocks_e_,
-              num_threads,
+              options_.num_threads,
               [values, bs, num_cols_e, x, y](int row_block_id) {
                 const int row_block_pos = bs->rows[row_block_id].block.position;
                 const int row_block_size = bs->rows[row_block_id].block.size;
@@ -163,10 +167,10 @@
                   // clang-format on
                 }
               });
-  ParallelFor(context,
+  ParallelFor(options_.context,
               num_row_blocks_e_,
               num_row_blocks,
-              num_threads,
+              options_.num_threads,
               [values, bs, num_cols_e, x, y](int row_block_id) {
                 const int row_block_pos = bs->rows[row_block_id].block.position;
                 const int row_block_size = bs->rows[row_block_id].block.size;
@@ -190,7 +194,17 @@
     LeftMultiplyAndAccumulateE(const double* x, double* y) const {
   if (!num_col_blocks_e_) return;
   if (!num_row_blocks_e_) return;
+  if (options_.num_threads == 1) {
+    LeftMultiplyAndAccumulateESingleThreaded(x, y);
+  } else {
+    CHECK(options_.context != nullptr);
+    LeftMultiplyAndAccumulateEMultiThreaded(x, y);
+  }
+}
 
+template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
+void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
+    LeftMultiplyAndAccumulateESingleThreaded(const double* x, double* y) const {
   const CompressedRowBlockStructure* bs = matrix_.block_structure();
 
   // Iterate over the first num_row_blocks_e_ row blocks, and multiply
@@ -214,28 +228,19 @@
 
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
-    LeftMultiplyAndAccumulateE(const double* x,
-                               double* y,
-                               ContextImpl* context,
-                               int num_threads) const {
-  if (!num_col_blocks_e_) return;
-  if (!num_row_blocks_e_) return;
-
+    LeftMultiplyAndAccumulateEMultiThreaded(const double* x, double* y) const {
   auto transpose_bs = matrix_.transpose_block_structure();
-  if (transpose_bs == nullptr || num_threads == 1) {
-    LeftMultiplyAndAccumulateE(x, y);
-    return;
-  }
+  CHECK(transpose_bs != nullptr);
 
   // Local copies of class members in order to avoid capturing pointer to the
   // whole object in lambda function
   auto values = matrix_.values();
   const int num_row_blocks_e = num_row_blocks_e_;
   ParallelFor(
-      context,
+      options_.context,
       0,
       num_col_blocks_e_,
-      num_threads,
+      options_.num_threads,
       [values, transpose_bs, num_row_blocks_e, x, y](int row_block_id) {
         int row_block_pos = transpose_bs->rows[row_block_id].block.position;
         int row_block_size = transpose_bs->rows[row_block_id].block.size;
@@ -254,14 +259,24 @@
               y + row_block_pos);
         }
       },
-      transpose_bs->rows.data(),
-      [](const CompressedRow& row) { return row.cumulative_nnz; });
+      e_cols_partition());
 }
 
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
     LeftMultiplyAndAccumulateF(const double* x, double* y) const {
   if (!num_col_blocks_f_) return;
+  if (options_.num_threads == 1) {
+    LeftMultiplyAndAccumulateFSingleThreaded(x, y);
+  } else {
+    CHECK(options_.context != nullptr);
+    LeftMultiplyAndAccumulateFMultiThreaded(x, y);
+  }
+}
+
+template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
+void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
+    LeftMultiplyAndAccumulateFSingleThreaded(const double* x, double* y) const {
   const CompressedRowBlockStructure* bs = matrix_.block_structure();
 
   // Iterate over row blocks, and if the row block is in E, then
@@ -307,26 +322,19 @@
 
 template <int kRowBlockSize, int kEBlockSize, int kFBlockSize>
 void PartitionedMatrixView<kRowBlockSize, kEBlockSize, kFBlockSize>::
-    LeftMultiplyAndAccumulateF(const double* x,
-                               double* y,
-                               ContextImpl* context,
-                               int num_threads) const {
-  if (!num_col_blocks_f_) return;
+    LeftMultiplyAndAccumulateFMultiThreaded(const double* x, double* y) const {
   auto transpose_bs = matrix_.transpose_block_structure();
-  if (transpose_bs == nullptr || num_threads == 1) {
-    LeftMultiplyAndAccumulateF(x, y);
-    return;
-  }
+  CHECK(transpose_bs != nullptr);
   // Local copies of class members  in order to avoid capturing pointer to the
   // whole object in lambda function
   auto values = matrix_.values();
   const int num_row_blocks_e = num_row_blocks_e_;
   const int num_cols_e = num_cols_e_;
   ParallelFor(
-      context,
+      options_.context,
       num_col_blocks_e_,
       num_col_blocks_e_ + num_col_blocks_f_,
-      num_threads,
+      options_.num_threads,
       [values, transpose_bs, num_row_blocks_e, num_cols_e, x, y](
           int row_block_id) {
         int row_block_pos = transpose_bs->rows[row_block_id].block.position;
@@ -362,8 +370,7 @@
               y + row_block_pos - num_cols_e);
         }
       },
-      transpose_bs->rows.data(),
-      [](const CompressedRow& row) { return row.cumulative_nnz; });
+      f_cols_partition());
 }
 
 // Given a range of columns blocks of a matrix m, compute the block
diff --git a/internal/ceres/partitioned_matrix_view_template.py b/internal/ceres/partitioned_matrix_view_template.py
index 27f368f..8978d63 100644
--- a/internal/ceres/partitioned_matrix_view_template.py
+++ b/internal/ceres/partitioned_matrix_view_template.py
@@ -132,7 +132,7 @@
 #ifndef CERES_RESTRICT_SCHUR_SPECIALIZATION
 """
 FACTORY = """  return std::make_unique<PartitionedMatrixView<%s,%s, %s>>(
-                   matrix, options.elimination_groups[0]);"""
+                   options, matrix);"""
 
 FACTORY_FOOTER = """
 #endif
@@ -142,7 +142,7 @@
   return std::make_unique<PartitionedMatrixView<Eigen::Dynamic,
                                                 Eigen::Dynamic,
                                                 Eigen::Dynamic>>(
-      matrix, options.elimination_groups[0]);
+      options, matrix);
 };
 
 }  // namespace ceres::internal
diff --git a/internal/ceres/partitioned_matrix_view_test.cc b/internal/ceres/partitioned_matrix_view_test.cc
index 4653e06..fc1ebdd 100644
--- a/internal/ceres/partitioned_matrix_view_test.cc
+++ b/internal/ceres/partitioned_matrix_view_test.cc
@@ -76,76 +76,6 @@
       std::uniform_real_distribution<double>(0.0, 1.0);
 };
 
-TEST_F(PartitionedMatrixViewTest, DimensionsTest) {
-  EXPECT_EQ(pmv_->num_col_blocks_e(), num_eliminate_blocks_);
-  EXPECT_EQ(pmv_->num_col_blocks_f(), num_cols_ - num_eliminate_blocks_);
-  EXPECT_EQ(pmv_->num_cols_e(), num_eliminate_blocks_);
-  EXPECT_EQ(pmv_->num_cols_f(), num_cols_ - num_eliminate_blocks_);
-  EXPECT_EQ(pmv_->num_cols(), A_->num_cols());
-  EXPECT_EQ(pmv_->num_rows(), A_->num_rows());
-}
-
-TEST_F(PartitionedMatrixViewTest, RightMultiplyAndAccumulateE) {
-  Vector x1(pmv_->num_cols_e());
-  Vector x2(pmv_->num_cols());
-  x2.setZero();
-
-  for (int i = 0; i < pmv_->num_cols_e(); ++i) {
-    x1(i) = x2(i) = RandDouble();
-  }
-
-  Vector y1 = Vector::Zero(pmv_->num_rows());
-  pmv_->RightMultiplyAndAccumulateE(x1.data(), y1.data());
-
-  Vector y2 = Vector::Zero(pmv_->num_rows());
-  A_->RightMultiplyAndAccumulate(x2.data(), y2.data());
-
-  for (int i = 0; i < pmv_->num_rows(); ++i) {
-    EXPECT_NEAR(y1(i), y2(i), kEpsilon);
-  }
-}
-
-TEST_F(PartitionedMatrixViewTest, RightMultiplyAndAccumulateF) {
-  Vector x1(pmv_->num_cols_f());
-  Vector x2 = Vector::Zero(pmv_->num_cols());
-
-  for (int i = 0; i < pmv_->num_cols_f(); ++i) {
-    x1(i) = RandDouble();
-    x2(i + pmv_->num_cols_e()) = x1(i);
-  }
-
-  Vector y1 = Vector::Zero(pmv_->num_rows());
-  pmv_->RightMultiplyAndAccumulateF(x1.data(), y1.data());
-
-  Vector y2 = Vector::Zero(pmv_->num_rows());
-  A_->RightMultiplyAndAccumulate(x2.data(), y2.data());
-
-  for (int i = 0; i < pmv_->num_rows(); ++i) {
-    EXPECT_NEAR(y1(i), y2(i), kEpsilon);
-  }
-}
-
-TEST_F(PartitionedMatrixViewTest, LeftMultiplyAndAccumulate) {
-  Vector x = Vector::Zero(pmv_->num_rows());
-  for (int i = 0; i < pmv_->num_rows(); ++i) {
-    x(i) = RandDouble();
-  }
-
-  Vector y = Vector::Zero(pmv_->num_cols());
-  Vector y1 = Vector::Zero(pmv_->num_cols_e());
-  Vector y2 = Vector::Zero(pmv_->num_cols_f());
-
-  A_->LeftMultiplyAndAccumulate(x.data(), y.data());
-  pmv_->LeftMultiplyAndAccumulateE(x.data(), y1.data());
-  pmv_->LeftMultiplyAndAccumulateF(x.data(), y2.data());
-
-  for (int i = 0; i < pmv_->num_cols(); ++i) {
-    EXPECT_NEAR(y(i),
-                (i < pmv_->num_cols_e()) ? y1(i) : y2(i - pmv_->num_cols_e()),
-                kEpsilon);
-  }
-}
-
 TEST_F(PartitionedMatrixViewTest, BlockDiagonalEtE) {
   std::unique_ptr<BlockSparseMatrix> block_diagonal_ee(
       pmv_->CreateBlockDiagonalEtE());
@@ -174,134 +104,122 @@
   EXPECT_NEAR(block_diagonal_ff->values()[2], 37.0, kEpsilon);
 }
 
-const int kMaxNumThreads = 8;
-class PartitionedMatrixViewParallelTest : public ::testing::TestWithParam<int> {
+// Param = <problem_id, num_threads>
+using Param = ::testing::tuple<int, int>;
+
+static std::string ParamInfoToString(testing::TestParamInfo<Param> info) {
+  Param param = info.param;
+  std::stringstream ss;
+  ss << ::testing::get<0>(param) << "_" << ::testing::get<1>(param);
+  return ss.str();
+}
+
+class PartitionedMatrixViewSpMVTest : public ::testing::TestWithParam<Param> {
  protected:
-  static const int kNumProblems = 3;
   void SetUp() final {
-    int problem_ids[kNumProblems] = {2, 4, 6};
-    for (int i = 0; i < kNumProblems; ++i) {
-      auto problem = CreateLinearLeastSquaresProblemFromId(problem_ids[i]);
-      CHECK(problem != nullptr);
-      SetUpMatrix(i, problem.get());
-    }
-
-    context_.EnsureMinimumThreads(kMaxNumThreads);
-  }
-
-  void SetUpMatrix(int id, LinearLeastSquaresProblem* problem) {
-    A_[id] = std::move(problem->A);
-    auto& A = A_[id];
-    auto block_sparse = down_cast<BlockSparseMatrix*>(A.get());
+    const int problem_id = ::testing::get<0>(GetParam());
+    const int num_threads = ::testing::get<1>(GetParam());
+    auto problem = CreateLinearLeastSquaresProblemFromId(problem_id);
+    CHECK(problem != nullptr);
+    A_ = std::move(problem->A);
+    auto block_sparse = down_cast<BlockSparseMatrix*>(A_.get());
     block_sparse->AddTransposeBlockStructure();
 
-    num_cols_[id] = A->num_cols();
-    num_rows_[id] = A->num_rows();
-    num_eliminate_blocks_[id] = problem->num_eliminate_blocks;
-    LinearSolver::Options options;
-    options.elimination_groups.push_back(num_eliminate_blocks_[id]);
-    pmv_[id] = PartitionedMatrixViewBase::Create(options, *block_sparse);
+    options_.num_threads = num_threads;
+    options_.context = &context_;
+    options_.elimination_groups.push_back(problem->num_eliminate_blocks);
+    pmv_ = PartitionedMatrixViewBase::Create(options_, *block_sparse);
+
+    EXPECT_EQ(pmv_->num_col_blocks_e(), problem->num_eliminate_blocks);
+    EXPECT_EQ(pmv_->num_col_blocks_f(),
+              block_sparse->block_structure()->cols.size() -
+                  problem->num_eliminate_blocks);
+    EXPECT_EQ(pmv_->num_cols(), A_->num_cols());
+    EXPECT_EQ(pmv_->num_rows(), A_->num_rows());
   }
 
   double RandDouble() { return distribution_(prng_); }
 
+  LinearSolver::Options options_;
   ContextImpl context_;
-  int num_rows_[kNumProblems];
-  int num_cols_[kNumProblems];
-  int num_eliminate_blocks_[kNumProblems];
-  std::unique_ptr<SparseMatrix> A_[kNumProblems];
-  std::unique_ptr<PartitionedMatrixViewBase> pmv_[kNumProblems];
+  std::unique_ptr<LinearLeastSquaresProblem> problem_;
+  std::unique_ptr<SparseMatrix> A_;
+  std::unique_ptr<PartitionedMatrixViewBase> pmv_;
+  int num_cols_e;
   std::mt19937 prng_;
   std::uniform_real_distribution<double> distribution_ =
       std::uniform_real_distribution<double>(0.0, 1.0);
 };
 
-TEST_P(PartitionedMatrixViewParallelTest, RightMultiplyAndAccumulateEParallel) {
-  const int kNumThreads = GetParam();
-  for (int p = 0; p < kNumProblems; ++p) {
-    auto& pmv = pmv_[p];
-    auto& A = A_[p];
+TEST_P(PartitionedMatrixViewSpMVTest, RightMultiplyAndAccumulateE) {
+  Vector x1(pmv_->num_cols_e());
+  Vector x2(pmv_->num_cols());
+  x2.setZero();
 
-    Vector x1(pmv->num_cols_e());
-    Vector x2(pmv->num_cols());
-    x2.setZero();
+  for (int i = 0; i < pmv_->num_cols_e(); ++i) {
+    x1(i) = x2(i) = RandDouble();
+  }
 
-    for (int i = 0; i < pmv->num_cols_e(); ++i) {
-      x1(i) = x2(i) = RandDouble();
-    }
+  Vector expected = Vector::Zero(pmv_->num_rows());
+  A_->RightMultiplyAndAccumulate(x2.data(), expected.data());
 
-    Vector y1 = Vector::Zero(pmv->num_rows());
-    pmv->RightMultiplyAndAccumulateE(
-        x1.data(), y1.data(), &context_, kNumThreads);
+  Vector actual = Vector::Zero(pmv_->num_rows());
+  pmv_->RightMultiplyAndAccumulateE(x1.data(), actual.data());
 
-    Vector y2 = Vector::Zero(pmv->num_rows());
-    A->RightMultiplyAndAccumulate(x2.data(), y2.data());
-
-    for (int i = 0; i < pmv->num_rows(); ++i) {
-      EXPECT_NEAR(y1(i), y2(i), kEpsilon);
-    }
+  for (int i = 0; i < pmv_->num_rows(); ++i) {
+    EXPECT_NEAR(actual(i), expected(i), kEpsilon);
   }
 }
 
-TEST_P(PartitionedMatrixViewParallelTest, RightMultiplyAndAccumulateFParallel) {
-  const int kNumThreads = GetParam();
-  for (int p = 0; p < kNumProblems; ++p) {
-    auto& pmv = pmv_[p];
-    auto& A = A_[p];
-    Vector x1(pmv->num_cols_f());
-    Vector x2(pmv->num_cols());
-    x2.setZero();
+TEST_P(PartitionedMatrixViewSpMVTest, RightMultiplyAndAccumulateF) {
+  Vector x1(pmv_->num_cols_f());
+  Vector x2(pmv_->num_cols());
+  x2.setZero();
 
-    for (int i = 0; i < pmv->num_cols_f(); ++i) {
-      x1(i) = x2(i + pmv->num_cols_e()) = RandDouble();
-    }
+  for (int i = 0; i < pmv_->num_cols_f(); ++i) {
+    x1(i) = x2(i + pmv_->num_cols_e()) = RandDouble();
+  }
 
-    Vector y1 = Vector::Zero(pmv->num_rows());
-    pmv->RightMultiplyAndAccumulateF(
-        x1.data(), y1.data(), &context_, kNumThreads);
+  Vector actual = Vector::Zero(pmv_->num_rows());
+  pmv_->RightMultiplyAndAccumulateF(x1.data(), actual.data());
 
-    Vector y2 = Vector::Zero(pmv->num_rows());
-    A->RightMultiplyAndAccumulate(x2.data(), y2.data());
+  Vector expected = Vector::Zero(pmv_->num_rows());
+  A_->RightMultiplyAndAccumulate(x2.data(), expected.data());
 
-    for (int i = 0; i < pmv->num_rows(); ++i) {
-      EXPECT_NEAR(y1(i), y2(i), kEpsilon);
-    }
+  for (int i = 0; i < pmv_->num_rows(); ++i) {
+    EXPECT_NEAR(actual(i), expected(i), kEpsilon);
   }
 }
 
-TEST_P(PartitionedMatrixViewParallelTest, LeftMultiplyAndAccumulateParallel) {
-  const int kNumThreads = GetParam();
-  for (int p = 0; p < kNumProblems; ++p) {
-    auto& pmv = pmv_[p];
-    auto& A = A_[p];
-    Vector x = Vector::Zero(pmv->num_rows());
-    for (int i = 0; i < pmv->num_rows(); ++i) {
-      x(i) = RandDouble();
-    }
-    Vector x_pre = x;
+TEST_P(PartitionedMatrixViewSpMVTest, LeftMultiplyAndAccumulate) {
+  Vector x = Vector::Zero(pmv_->num_rows());
+  for (int i = 0; i < pmv_->num_rows(); ++i) {
+    x(i) = RandDouble();
+  }
+  Vector x_pre = x;
 
-    Vector y = Vector::Zero(pmv->num_cols());
-    Vector y1 = Vector::Zero(pmv->num_cols_e());
-    Vector y2 = Vector::Zero(pmv->num_cols_f());
+  Vector expected = Vector::Zero(pmv_->num_cols());
+  Vector e_actual = Vector::Zero(pmv_->num_cols_e());
+  Vector f_actual = Vector::Zero(pmv_->num_cols_f());
 
-    A->LeftMultiplyAndAccumulate(x.data(), y.data());
-    pmv->LeftMultiplyAndAccumulateE(
-        x.data(), y1.data(), &context_, kNumThreads);
-    pmv->LeftMultiplyAndAccumulateF(
-        x.data(), y2.data(), &context_, kNumThreads);
+  A_->LeftMultiplyAndAccumulate(x.data(), expected.data());
+  pmv_->LeftMultiplyAndAccumulateE(x.data(), e_actual.data());
+  pmv_->LeftMultiplyAndAccumulateF(x.data(), f_actual.data());
 
-    for (int i = 0; i < pmv->num_cols(); ++i) {
-      EXPECT_NEAR(y(i),
-                  (i < pmv->num_cols_e()) ? y1(i) : y2(i - pmv->num_cols_e()),
-                  kEpsilon);
-    }
+  for (int i = 0; i < pmv_->num_cols(); ++i) {
+    EXPECT_NEAR(expected(i),
+                (i < pmv_->num_cols_e()) ? e_actual(i)
+                                         : f_actual(i - pmv_->num_cols_e()),
+                kEpsilon);
   }
 }
 
-INSTANTIATE_TEST_SUITE_P(ParallelProducts,
-                         PartitionedMatrixViewParallelTest,
-                         ::testing::Values(1, 2, 4, 8),
-                         ::testing::PrintToStringParamName());
+INSTANTIATE_TEST_SUITE_P(
+    ParallelProducts,
+    PartitionedMatrixViewSpMVTest,
+    ::testing::Combine(::testing::Values(2, 4, 6),
+                       ::testing::Values(1, 2, 3, 4, 5, 6, 7, 8)),
+    ParamInfoToString);
 
 }  // namespace internal
 }  // namespace ceres